Thierry Poidras

Pierre Camps
email: camps@gm.univ-montp2.fr

Patrick Nicol

Controlled atmosphere Vibrating Thermo-Magnetometer (C at VTM) : A new device to optimize the absolute paleointensity determinations

The laboratory of paleomagnetism of Montpellier (France) developed a new one-axis Vibrating Thermal Magnetometer dedicated to the study of physical properties of natural rocks remanence. Among its key characteristics, this apparatus allows both to measure the magnetization moment on the interval from the room temperature to 700 • C with a precision of 2x10 -9 Am 2 , and to acquire a total or a partial TRM using a steady field from -100 up to 100 µT. Another point that is worth noting is that one can apply a controlled atmosphere by means of argon flux in order to prevent oxidation of the studied sample during heating. In the present paper, we report a technical description of this new instrument and review some specific applications in absolute paleointensity surveys.

Introduction

In 1994, the paleomagnetic laboratory in Montpellier bought a vibrating magnetometer with translation manufactured by Orion (Borok, Russia). The principal characteristic of this one-axis magnetometer is to measure continuously, during heating-cooling cycles between the room temperature and 700 • C, the remanent magnetization recorded in rock samples of ∼1 cm 3 size. A second specification is that continuous measurements are carried out either in a null magnetic field (demagnetization curve) or in presence of a weak-magnetic field between -100 and +100µT applied along one axis of the sample (magnetization curve). It is obvious that such apparatus is of great interest in paleomagnetic surveys since it allows to test directly the physical properties of the Natural Remanent Magnetizations (NRM). In particular, it opens new prospects in different fields such as:

1. The selection of well-suited samples for paleointensity determinations [START_REF] Camps | Paleomagnetic and geochronological study of a geomagnetic field reversal or excursion recorded in pliocene volcanic rocks from Georgia (Lesser Caucasus)[END_REF][START_REF] Carvallo | Mono Lake or Laschamp geomagnetic event recorded from lava flows in Amsterdam Island (southeastern Indian Ocean)[END_REF].

2. The development of new methodologies in paleointensity [START_REF] Biggin | The behaviour and detection of partial thermal thermoremanent magnetisation (ptrm) tails in Thellier palaeointensity experiments[END_REF][START_REF] Camps | Matuayma-Brunhes reversal on Tjornes Peninsula (northern Iceland) revisited[END_REF].

1
3. The experimental investigation of physical properties of thermo, thermoviscous, and chemical remanent magnetizations, according to the size of magnetic minerals [START_REF] Draeger | Single-domain chemical, thermochemical and thermal remanences in a basaltic rock[END_REF][START_REF] Biggin | First-order symmetry of weak-field partial thermoremanence in multi-domain ferromagnetic grains. 1. experimental evidence and physical implications[END_REF].

Unfortunately the former apparatus was failing in term of reliability, since we encountered multiple breakdowns of various origins, and unsatisfactory in term of sensitivity especially at high temperature, since some slightly magnetized basalts could not be analyzed. Consequently, during more than 10 years, we have worked on its development and upgrade with the objective to improve both its reliability and its sensitivity. Today, we have a new prototype called C at VTM in which only the detecting coils and the field application coils are of origin. This apparatus supplement the list of laboratory-made magnetometer, dedicated to rock magnetic studies, able to measure magnetic moments at elevated temperature [see 7, for a review].

In the first part of the present paper, we will report a technical description of Montpellier instrument by emphasizing its technical characteristics. Next, we will review some specific applications in absolute paleointensity surveys. Finally, we will conclude on what could be its future developments.

Magnetic Moment Measurement in a Vibrating Magnetometer

The principle of a vibrating magnetometer is to measure a voltage induced in a pick-up coil assembly by a changing magnetic flux. Usually, the changing flux is achieved by a periodically change of the sample position along the vibrating axis within a fixed detecting coil assembly. The motion of a magnetic sample produces a magnetic flux through the coil as:

Φ B = G • M
where Φ B is the magnetic flux in Webers, G is the vectorial geometry factor in Wb/Am 2 , and M is the magnetic moment of the sample in Am 2 . The voltage induced in the coils by the change in times of the magnetic flux is given by the Faraday's law of induction: v = -dΦ B dt where v is the voltage in volts, and t is the time. For vibrating magnetometer with a linear translation of M following z-axis, an oscillatory voltage is generated as:

v = -M
dG dz dz dt This equation helps us to identify three important parameters to achieve the best measurement performances. The first is the sample signal M, which is the sum of the remanent magnetic moment of the sample M r and an induced moment M i which is only present if a residual magnetic field exists in the measurement zone. Provided that M i is negligible, there is nothing we can do on sample signal except to choose a sample volume as large as possible. The second, and certainly the most important, is the variation of the vectorial geometry factor along the vibrating axis (dG/dz). This term characterizes the dependence of the signal induced in the coils upon the sample-coil system geometry. The third represents the sample motion (dz/dt). It depends on the quality of the mechanical system and on its control. In a vibrating magnetometer, the sample undergoes a sinusoidal vibration such as:

z = z 0 + a cos(ωt)
where z is the sample position along the vibration axis, a and ω are the amplitude and the angular frequency of the vibrations, respectively, and z 0 is the distance between the center of the motion and the center of the coil assembly. Then, we have:

dz dt = aω sin(ωt)
This equation clearly shows that, for a given sample magnetic moment, a larger vibration amplitude a as well as a higher frequency of vibration ω will increase the voltage induced in the coils, and thus will increase the magnetometer sensitivity.

C at VTM: Description and Solutions

The C at VTM magnetometer is a horizontal assembly composed by several main parts described hereafter and illustrated on Figure 1.

The sensor assembly

We kept the sensor assembly of origin as manufactured by the Orion company (Borok, Russia). This assembly is made of two coils in series-opposition centered in two concentric mumetal shields. Their physical characteristics are as following:

• Outer(Inner) coil diameter: 74.8(37.8) mm

• Coil width: 17 mm

• Outer(Inner) mumetal shield: 164(144) mm length and 100(77) mm diameter.

• Resistance of the two coils in series: 14955 Ω

• Wire diameter : 0.11 mm.

G(z) profile determination:

The knowledge of the vectorial geometry factor along z-axis is fundamental to fully describe the sensor assembly. The main reason is practical: the zone in which the sample is translated must be chosen such as dG/dz be as constant as possible over all the sample volume during its movement. We determined the G(z) profile by means of two different approaches. The first was to solve numerically the analytical solution for our sensor assembly geometry (see Appendix A). The second was to measure experimentally this profile. This is simply done by measuring the magnetic field generated along the vibrating z-axis when the sensor coil assembly is supplied with a known current I. G(z) profile is then estimated using the relation :

B(z) = G(z) • I
where B is the magnetic induction in Tesla, and I the current passing through the wires in Ampere. We designed and built a specific probe to achieve this measurement. This probe is described in section 3.2. We find, using a 0.01 mA DC current, a maximum field of 4.15 µT for z = -23.5 mm and a minimum field of -4.15 µT for z = 23.5 mm. The zone where (dG/dz) is nearly constant and offers a good range for sample motion is bounded between -13 mm and 13 mm. Over this amplitude range, the slope of (dG/dz) is 26.21 T/Am (Figure 2). Setting the value of the vibration amplitude a to 8.14 mm and the vibration frequency f to 8 Hz, with ω = 2.π.f , we can calculate a peak-value of 0.4092 m/s for (dz/dt). Then, the value of the magnetic moment measured by the lock-in amplifier (model SR830) and the preamplifier is:

M rms = - 1 K × v rms (dG/dz) (dz/dt) with, M rms = M √ 2
where M rms is the root mean square of the magnetic moment in Am 2 , K is the gain of the preamplifier, here K=6675, and M is the sample magnetic moment.

Offset and Noise:

The output of the lock-in amplifier v rms is given by:

v rms = vp of f set + vp noise +K -(M m ).(dG/dz).(a m ω m )/ √ 2 +vs noise
where vp of f set and vp noise are the constant offset and the constant noise of the preamplifier, respectively, K is the gain of the preamplifier, M n is the measured magnetic moment (M m = M + M noise ), a m is the amplitude of the vibration (a m = a + a noise ), ω m is the circular frequency (ω m = ω + ω noise ), and vs noise is the Johnson noise of the sensor coils. The offset and noise determinations were achieved experimentally using first a vibrating coil with no current (I=0 A) giving a null magnetic moment (M=0 Am 2 ). Setting I=0 A, we have:

v rms = vp of f set + vp noise + K.vs noise
The experiment lasted 20 minutes during which 776 data were acquired. We averaged the 20-mn-signal by calculating the arithmetic mean, and found 66±47 µV. The mean value corresponds to the offset of the preamplifier, while the standard deviation is the noise of the preamplifier plus the amplified noise of the sensor. A value of 47 µV corresponds to 0.6x10 -9 Am 2 rms or 3.1x10 -9 Am 2 peak-to-peak noise. Then, since the noise depends on the magnetic moment, we repeated the same experiment for different values of current through the calibration coil.

Johnson noise:

The detection coils produce a Johnson noise, generated by the thermal agitation of electrons, which happens regardless of any applied voltage. The rms value of the voltage across a resistor due to the Johnson noise, expressed in volts per root hertz, is given by: ν n = 4k b T R∆f where k B is the Boltzmann's constant in joules per kelvin, T is the coil temperature in kelvins, R is the resistor value in ohms, and ∆f is the bandwidth in Hertz over which the noise is measured. The amount of noise measured by the SR830 lock-in amplifier is given by:

v noise = 0.13 √ R W
where v noise is the rms voltage in nVolts, R is the resistance of the sensor coil system, and W is the equivalent noise bandwidth of the SR830 low-pass filter.

In our vibrating magnetometer, with R = 14955 Ω and W = 0.026, the minimum rms-noise we could achieve is theoretically 2.56 nV. This value corresponds to 0.34x10 -9 A.m 2 or 1.7x10 -9 A.m 2 peak-to-peak noise as notified in the SR830 manual (page 3.17).

Calibration:

As shown on Figure 2, harmonic distortions are suspected because (dG/dz) is not perfectly constant over the amplitude range. We quantified the ratedistortion D of the output waveform for two vibration amplitudes, 8.14 and 3.61 mm peak-values, keeping the vibration frequency at a constant value of 8 Hz. Results are presented in Table 1 for the first five harmonics. Note that a value of 8.14 mm is preferred for a general use, although the rate-distortion is higher than for a value of 3.61 mm, because the corresponding (dG/dz).(dz/dt) is higher (10.72 compared to 4.75). Then, one needs to calibrate the sensor in order to compensate the system non-linearities. One way to perform the calibration is to build a dedicated coil. Paperno et al. [START_REF] Paperno | Cylindrical induction coil to accurately imitate the ideal magnetic dipole[END_REF] described how to design a cylindrical induction coil to accurately simulate an ideal magnetic dipole. The calibration coil we built is 6.35 mm length and 9.5 mm diameter. This coil produces a magnetic moment of 1.488x10 -3 A.m 2 for one Ampere. An alternative solution is to perform a cross calibration between our 2G cryogenic magnetometer and the C at VTM by measuring on both systems a remanent magnetization imparted to a sample in the laboratory. This double measurement is repeated for several steps of AF-demagnetization until the limits of the C at VTM are reached.

Signal-to-noise ratio:

We use our calibration coil to determine the signal-to-noise ratio (STN). Each experiment lasts 20 minutes during which a constant current is applied. STN is given by: ST N = 20 log v signal v noise where v signal is the average of the output voltage, and v noise is the corresponding standard deviation. It is usually admitted that measurements with a signal-tonoise ratio lower than 10 dB are inexploitable. Results reported in Table 2 show that meaningful measurements can be performed from a sample magnetic moment as low as 2x10 -9 Am 2 (15dB). These results take into account the modulation "noise" induced by the variation of amplitude, which is estimated to be about ±50 µm from the equation in section 3.1.2, and by the variation of frequency, which is f n = ± 0.001 Hz.

Others sources of noise:

Additional noises on magnetic moment measurements are due to the resistive heater and to the solenoid applying the field on the sample. These noises are difficult to be estimated. During heatings, the physical and chemical properties of the magnetic minerals may suffer from changes. These variations can be considered as a noise in the way that we loose informations in the signal due to uncontrolled changes in the magnetic properties of the sample.

Magnetic shielding

The pick-up system is surrounded by two mumetal magnetic shields which provides a low reluctance path guiding the magnetic field around the measurement zone. An accurate estimation of the residual field is required. As the furnace inner diameter is only 17 mm, it is not obvious to measure the three components of the residual field in the measurement zone. We designed and built a specific 3-axes probe in order to achieve this measurement. We choose the Honeywell sensors HMC1001 and HMC1002. Those magnetoresistive sensors are able to detect a magnetic field as low as 3 nT for a frequency of 10 Hz. The sensors are mounted on a printed circuit board. Sensors control, signal acquisition and data processing are achieved with the MSC1211 component from Texas Instruments. The MSC1211 is a precision 24-bit analog-to-digital converter (ADC) with a 8051 microcontroller. This probe is a standalone system that can be used to measure the magnetic field for other needs in the laboratory.

The shielding factor or effectiveness S is defined as the ratio of the external field to the residual field within the shielded zone. Our system yields a good attenuation, S being measured between 175 and 250. Despite the mumetal magnetic shield, the field in the measurement zone is not totally cancelled. To obtain better results the system is positioned so that its axis be perpendicular to the ambient magnetic field. However there is still a little field measured between 100 and 200 nT which has to be cancelled.

Field application coils

A solenoid located between the pick-up coils and the sample enables the user to apply a field along the sample z-axis. This same solenoid is used to apply an offset to obtain a residual field as little as possible when a null field is required. The quality of the current source is a decisive factor. We use a National Instruments Digital to analog card followed by a voltage to current amplifier to drive the z-coil. It also appears that the transverse field in the measurement zone is not null and can be quite important depending on the quality of the mumetal shielding. We encountered problems with induced magnetization caused by this field. Two perpendicular Helmholtz coils assembly are wound between the magnetic shield and the pick-up coils to cancel this field. Two sources supply a small current to each coil to nullify the two components of the transverse field.

Sample Heater

The heating of samples is one part of the system to which we have to give attention. It is a resistive furnace made of a bifilar winding that is non-inductive in a theoretical point of view. We have to pay attention on how to wind in order to achieve the non-inductivity and to offer an homogeneous temperature profile in the measurement zone. The furnace is wound with Nikrothal 80 non-magnetic wire on a quartz tube of 17 mm inner diameter and 20 mm outer diameter. The total resistance is 75 Ohms and the maximum current is two amperes. A copper jacket and a thin layer of rockwool insulate the heater from the resistor coils. The jacket is designed so that the cooling liquid flow be regular on its whole length. The cooling circuit is closed-looped and the flow rate is regulated at 0.6 l/h. If for any reason the flow rate is to low, the security switches down the power supply. The power supply, designed and self-built in the laboratory, is a Pulse Width Modulation (PWM) supplied with a 150 V voltage-source. It is a full H-Bridge which enables a bi-directional current flow in the resistor wire. The specificity of this PWM power supply is that each current pulse is reversed, in order to minimize the effect of the imperfect non-inductive winding. The period of the PWM is 1 Khz. The temperature is measured with a thermocouple type E (-200 • C to 900 • C) or type S (0 • C to 1450 • C). The thermocouple is part of the sample holder and its exposed part is in direct contact with the sample. We use a Eurotherm controller to program the ramp rate and dwelling temperature needed to control the heating process.

The Sample Holder

The sample holder is made of three ceramics parts (Figure 3): the holder by itself is made in Macor. This machinable glass ceramic has a continuous use temperature of 800 • C and a peak temperature of 1000 • C. Its coefficient of thermal expansion readily matches most metals and sealing glasses. It exhibits zero porosity, and unlike ductile materials, wont deform. It is an excellent insulator at high voltages, various frequencies, and high temperatures. This holder minimizes the heat loss out of the heater. The bell that maintain the sample is made of Shapal. It is a machinable form of Aluminium Nitride ceramic with excellent mechanical strength and thermal conductivity. Shapal has zero porosity, good abilities to seal under vacuum, low thermal expansion coupled with a high heat resistance. Shapal also offers excellent machinability with conventional machine tools. The machined Shapal-bell (female) is screwed on the male screw thread of the Macor-holder. During the experiment they join together with a bond made of water and kaolin. We use some rockwool in the assembly to avoid direct contact between the sample and the bottom of the bell. Inside the Macor-holder, there is a ceramic rod (two holes) that contains the thermocouple.

A gas flow (argon) is sent through a hole at the base of the Macor-holder and get out of the Shapal-bell through a little hole. The flow is tuned at the flowmeter to be as the smallest rate as possible (∼ 0.05 l/mn). An electrical valve is driven by the computer to shut off the gas supply at the end of the experiment. A 70-cm-long fiber carbon rod is used to connect the sample holder to the motor.

Sample Motion Generator

As seen previously, the stability of both amplitude and frequency of the sample displacement is required to achieve a good measurement. This is why we used a linear motor. We chose the ThrustsTube Micro Motor model TB1104 from Copley Motion Systems. The mass of the system (forcer mass + sample holder + transmission rod) limits the speed performances of the assembly. The motion of 16.14 mm peak-to-peak at 8 Hz is a compromise we chose between sensitivity and long-time performance. Actually, some experiments can last more than 30 hours. The motor forcer moves along the axis of the thrust rod mass (Figure 4). It is heavier than the thrust rod mass but, as the thrust rod, is strongly magnetized and thus it has to be steady. The forcer is fixed to a low profile guide system. The carbon rod connected to the sample holder is fixed at the extremity of the system. The mechanical assembly has to be as simple as possible. We use to minimize the friction a miniature Drylin (linear bearing) system from Igus c . The two floating bearings used to guide the motor are very light and have a very low coefficient of friction. The mobile forcer is fixed to a plate to transmit the motion to the fiber carbon rod. A linear optical encoder module reads the position on a linear scale fixed on the mobile assembly. The resolution of the encoder assembly is 8.8 µm. We use a DC brushless digital servo-amplifier to drive the motor and program the sample displacement. The servo-amplifier is configured via a three-wire, full-duplex RS-232 port, by means of the CME2 software or with specific commands. The servo-amplifier provides a functiongenerator command to configure the amplitude and frequency of the motion. The stability of the vibration is ± 0.001 Hz for the frequency and ± 50 µm for the amplitude .

Signal Processing

As the pick-up coil assembly has a high resistance (15 KΩ) it behaves more than a current source than a voltage source. Thus, we used a current to voltage preamplifier followed by a low-pass filter. The output of the preamplifier is directly connected to a lock-in amplifier SR830 (Stanford Research System). The reference input comes from a photo-electric detector which is mounted on the motor assembly.

Data Processing

Labview is used to configure the experiment and to process the data. The computer drives the temperature controller, the applied field, the motor servoamplifier and the lock-in amplifier through the GPIB bus. Before to run the experiment, the user has to define a succession of temperature and applied field profiles. Each profile consists of a heating or cooling ramp rate (with field to apply) and the temperature to reach (with field to apply). There is no limitation in the number of profiles. During the experiment it is possible to remove or add a profile.

At the end of each profile the magnetization at the dwell temperature is computed. Each profile has its own datafile on which it is possible to work during the experiment. As soon as the experiment is finished, the results are available and can be plotted and analyzed.

4 Applications in Absolute Paleointensity Determination by the Thellier Method

The Paleointensity method and its limitations

The principle of absolute paleointensity determination with Thellier type experiments [START_REF] Thellier | Sur l'intensité du champ magnétique terrestre dans le passé historique et géologique[END_REF] rests on the modeling of the natural remanent magnetization (NRM) by an artificial thermo-remanent magnetization (TRM) acquired in the laboratory in presence of a known magnetic field H. Indeed, the methods are based on a direct comparison between the stepwise thermal demagnetization of the NRM and the acquisition of partial TRM (pTRM). Various aspect of the method applicability have been widely discussed [see 10, for a review and a description of the original Thellier method and its derivatives]. To be brief, igneous rocks used for paleointensity determinations must satisfy the following three conditions :

1. NRM must be a TRM not disturbed by significant secondary components.

2. The magnetic properties of the samples must be reasonably stable during the laboratory heating.

3. Thellier method is based on three assumptions concerning the properties of pTRM usually referred to as the Thellier laws [START_REF] Thellier | Sur l'aimantation des terres cuites et ses applications géophysiques[END_REF] which are valid only for fine remanence carriers, i.e. single-domain (SD) or small pseudo-singledomain (PSD) grains: (i) Additivity law: let two pTRMs imparted in the same magnetic field, pTRM(T 1 ,T 2 ) being acquired in the temperature interval T 1 -T 2 and pTRM(T 2 ,T 3 ) being acquired in the temperature interval T 2 -T 3 with T 1 > T 2 > T 3 then, pTRM(T 1 , T 2 ) + pTRM(T 2 , T 3 ) = pTRM(T 1 , T 3 ) (ii) Independence law: the blocking and unblocking temperatures of pTRMs must be equal. This means that a pTRM(T 1 ,T 2 ) is not thermally demagnetized when heated in zero field at a temperature lower than T 2 , while it is totally demagnetized when heated in zero field at a temperature equal or upper than T 1 .

(iii) The intensity of a pTRM(T 1 ,T 2 ) should not depend on the thermal history of the sample. This property can be checked for instance by comparing two different kinds of pTRM: a pTRM a (T 1 ,T 2 ) acquired when the upper temperature T 1 is reached by cooling from Curie Temperature (T c ), i.e. from above, and a pTRM b (T 1 ,T 2 ) when T 1 is reached by heating from T room , i.e. from below [START_REF] Vinogradov | Investigations in rock magnetism and palaeomagnetism, chapter On the effect of low temperature heating on the magnetic state of multidomain magnetite[END_REF][START_REF] Shcherbakova | Properties of partial thermoremanent magnetization in pseudosingle domain and multidomain magnetite grains[END_REF]. If the pTRMs are independent of the thermal history of the rock sample then,

pTRM a (T 1 , T 2 ) = pTRM b (T 1 , T 2 )
It is obvious that Thellier type experiments impose severe conditions on the samples that are rarely encountered in natural rocks. A direct consequence is that paleointensity determinations are very often characterized by a high failure rate, ranging between 70 and 90% [START_REF] Perrin | Paleointensity determination, magnetic domain structure and selection criteria[END_REF]. Hence, the number of determinations of absolute paleointensity available is extremely weak: To date approximately 3000 data, representative of instantaneous recordings of the paleomagnetic field, were published for the whole geological times. In addition to the low number of data available, their quality is unequal: 63% of the published data are eliminated with criteria of selection particularly not very severe and, in particular, almost all the determinations between 4 and 60 Ma [START_REF] Perrin | IAGA paleointensity database: distribution and quality of the data set[END_REF].

One of the specific purpose of using the C at VTM in Montpellier is to test directly the physical properties of NRM in order to provide a relevant preselection of samples that can be successfully used for a paleointensity determination. Examples of simple tests carried out to this end are reviewed hereafter.

Is the NRM a pure TRM ?

A positive answer to this question is the first prerequisite for obtaining a paleointensity. Usually, authors used directional arguments to detect the presence of significant secondary component in the remanence, as the grouping of NRM direction in the same cooling unit associated with a careful examination of the individual vector end-point diagrams obtained during the demagnetization in the paleodirection study. This question may be answered more directly by means of our C at VTM, provided that no noticeable chemical changes in the ferrimagnetic minerals occurred during the laboratory heatings. The experimental approach consists in comparing the shape of the continuous thermal demagnetization curve of the NRM with that of an artificial total TRM, imparted in a laboratory field H during the cooling from T c to T room on the same sample or on a sister sample from the same core. In this experiment, samples have to be drilled in the direction of the characteristic NRM due to the C at VTM limitation of measurement along a single axis. A similar shape of these two continuous demagnetization curves is a qualitative evidence that the NRM is actually a TRM (Figure 5) [START_REF] Camps | Paleomagnetic and geochronological study of a geomagnetic field reversal or excursion recorded in pliocene volcanic rocks from Georgia (Lesser Caucasus)[END_REF][START_REF] Carvallo | Mono Lake or Laschamp geomagnetic event recorded from lava flows in Amsterdam Island (southeastern Indian Ocean)[END_REF]. This involves as a corollary that a difference between the two curves suggests that the NRM is not a pure TRM, which is a sufficient condition to turn down samples from the selection. Moreover, the demagnetization curve of the NRM yields an estimate of the unblocking temperatures of NRM, which is a valuable information to choose the heating steps in the Thellier experiment.

It is obvious that a positive result in this experiment proves not only that the NRM is a TRM but also that the nature and/or the physical properties of the Fe-Ti oxides carrying the magnetization are not altered during heatings. Thus, condition 1 and condition 2 as stated above are checked simultaneously. In practice, the thermal stability is first verified independently by continuous measurements of magnetic susceptibility in weak field according to the temperature (KT curves). In Montpellier, we use the susceptibility meter Kappabridge KLY3 c associated to the furnace CS3, in which the heating-cooling cycles are performed under Argon flux in order to limit oxidation of the ferrimagnetic minerals. Because we believe important in the sample selection procedure to keep constant the experimental conditions, NRM and TRM demagnetization curves are measured in the C at VTM also under Argon flux.

Bolshakov and Shcherbakova's pTRM-tail test

The next step in sample selection for paleointensity determination is to assess the magnetic domain size in order to keep only sample with SD-like pTRM behavior. As stated in condition 3, the presence of multi-domain grains will invalidate the Thellier method. Usually, the size of the grains is estimated by the measurement of the hysteresis loop parameters for the whole rock. However this measurement often gives ambiguous answers between a behavior of PSD grain of that of a mixture of single-and multi-domain grains. Bolshakov et al. [START_REF] Bol | A thermomagnetic criterion for determining the domain structure of ferrimagnetics[END_REF] were the first to suggest that the domain structure can be inferred from a thermomagnetic criterion. This test, illustrated on Figure 6, consists of a direct verification of the Thellier law of independence of pTRM which implies that any pTRM(T 1 ,T 2 ), with T 1 > T 2 , must be entirely destroyed when heated en zero field strictly in the interval T 2 -T 1 . From an experimental study on synthetic powders of magnetite and hematite of various sizes, Shcherbakova et al. [START_REF] Shcherbakova | Properties of partial thermoremanent magnetization in pseudosingle domain and multidomain magnetite grains[END_REF] defined the A parameter by

A(T 1 , T 2 ) = tail[pTRM(T 1 , T 2 )] pTRM(T 1 , T 2 ) × 100%
as the relative intensity measured at room temperature of the pTRM remaining after heating to T 1 (pTRM-tail). According to this experimental study, this parameter can serve as a quantitative indicator of the domain structure of a sample. They found thresholds as following : the remanent carrier are predominantly SD grains for A(T 1 , T 2 ) < 4%, they behave as pseudo-single domain (PSD) grains for 4% < A(T 1 , T 2 ) < (15 -20)%, and they are predominantly MD grains for A(T 1 , T 2 ) > 20%. A major drawback of the Bolshakov [START_REF] Bol | A thermomagnetic criterion for determining the domain structure of ferrimagnetics[END_REF]'s pTRM-tail test is the initial heating to the Curie temperature necessary to demagnetized the NRM prior to the pTRM acquisition. Thus, this test can be applied only for samples thermally stable. The main application of determining the parameter A is usually to screen out samples ill-suited for paleointensity determinations [see 1, 14, for example]. Shcherbakov et al. [START_REF] Shcherbakov | Thermal stability of pTRMs created from different magnetic states[END_REF], arguing that this parameter is physically meaningful, even suggested to use it systematically before any paleointensity experiments. Further applications of A parameter have been proposed by Carvallo et al. [START_REF] Carvallo | Mono Lake or Laschamp geomagnetic event recorded from lava flows in Amsterdam Island (southeastern Indian Ocean)[END_REF] and Plenier et al. [START_REF] Plenier | Absolute palaeointensity of Oligocene (24-30 Ma) lava flows from the Kerguelen Archipelago (southern Indian Ocean)[END_REF] to help in the interpretation of NRM-TRM diagrams. In their study, they measured the coefficient A at increasing temperature intervals and used the results to validate their choice of the most suitable portion of the NRM-TRM diagrams when, in non-ideal cases, two slopes yielding a technically acceptable determination are present. Whatever the use of the A parameter, its determination when possible strengthens the reliability of a paleointensity determination.

Pseudo-Thellier Paleointensity

Perrin [START_REF] Perrin | Paleointensity determination, magnetic domain structure and selection criteria[END_REF] proposed a different approach to select suitable samples for paleointensity determination which do not required, as the pTRM-tail test, a total thermal demagnetization prior to the experiment. She proposed to perform using the C at VTM a Thellier-like experiment with the sole objective to estimate the shape of the NRM-TRM diagram and therefore the temperature range where a real Thellier experiment will be later performed. This method is called pseudo-Thellier because of the limitation of measuring the magnetization along one axis induces errors too much high to consider valid the paleointensity estimate. This conclusion remains true even if the small sample used in this experiment is drilled as close as possible along the characteristic NRM direction.

Shift of the Baseline in Thellier experiment

The idea is to perform a traditional Thellier [START_REF] Thellier | Sur l'intensité du champ magnétique terrestre dans le passé historique et géologique[END_REF] experiment except that the NRM left and the pTRM gained are measured at a temperature greater than the room temperature (for instance 200 • C). Remanence measurements are classically performed at room temperature just for purely technical reasons: paleointensity oven and magnetometer are usually two different apparatus. We report here an example of application of this original method on Icelandic lava flows from Tjornes peninsula, which are suspected to have recorded the Matuyama-Bruhnes transition [START_REF] Kristjansson | Brunhes-matuyama paleomagnetism in three lava sections in iceland[END_REF]. A very detailed rock magnetic study (low-and high-temperature susceptibility curves, Forc diagrams, thin sections analysis) revealed the joint presence of two titanomagnetite populations, the first, which has undergone a high temperature oxidation, shows a single-domain (SD) behavior and a mean curie temperature above 500 • C, while the second, not oxidized, presents a multidomain (MD) behavior and a mean curie temperature around 150 • C [START_REF] Camps | Matuayma-Brunhes reversal on Tjornes Peninsula (northern Iceland) revisited[END_REF]. The presence of MD grains does not enable us to obtain paleointensity estimates with a conventional method. However, by using the C at VTM and by choosing a temperature of measurement higher than the Curie temperature of the MD population, e.g., 200 • C, we obtained paleointensity diagrams of very good technical quality (Figure 7).

Conclusion

The C at VTM is an apparatus which should improve the success rate in absolute paleointensity determinations by investigating directly the physical properties of the NRM carried by natural rocks. Its use allows to efficiently preselect wellsuited samples for paleointensity experiments. Beyond this practical aspect, some experiments as the ones presented in this paper strengthen the reliability of the paleointensity estimates. This why we believe important to continue to work on further developments. To date, improvements can be made on each step of the measurement process. It is still possible to work on more sensitive coils system, and also on three-axis systems [START_REF] Goff | A new three-axis vibrating magnetometer for continuous high-temperature magnetization measurements[END_REF]. The use of a very low-noise preamplifier will improve the signal-to-noise ratio of the system by decreasing the noise at the source. As linear motors are also improving, it should be possible to speed up significantly the frequency of the sample motion. Another way to improve the signal-to-noise ratio is to increase the induced signal by going out of the low-frequency electronic noise area. The resistive heater is a important source of noise. It is an exciting challenge to find other ways to heat the sample by means of an external heat source as performed in the prototype of magnetometer developed by Wack et al. [START_REF] Wack | A new type of three-component spinner magnetometer to measure the remanence of rocks at elevated temperature[END_REF].

Acknowledgments

This work was initiated by Michel Prévot. He bought the original version of the vibrating magnetometer and did a lot of work to improve it. Without him the CatVTM would not have been built. His work is gratefully acknowledged. Reviews by Maxime LeGoff and Juan Morales helped to clarify several issues on the original manuscript and are much appreciated. We also benefited from discussions with Mireille Perrin and Brigitte Smith.

Appendix

A1: Analytical solution for G(z)

We present here the analytical solution used to calculate the G(z) profile in the CatVTM. In this vibrating magnetometer, the sensor assembly is made of two coils in series-opposition (e.g., Smith's arrangement). The vectorial geometry factor is then:

G(z) = G 1 (z) -G 2 (z)
where the subscripts 1 and 2 denote the first and the second coil, respectively. For the first coil, G 1 (z) is determined by means of the equation given by Durand [START_REF] Durand | Magnétostatique[END_REF]:

B(z) = J • µ 0 2 b a (cosγ 1 ± cos γ 2 )dρ
where B is the magnetic induction in Tesla, J is the current density in A/m 2 , a and b are the coil inner-and outer-radius in meters, respectively, and γ 1 and γ 2 two angles as illustrated on Figure 8. Then

G 1 (z) = B(z)
J is expressed in Tm 2 /A. To be directly comparable with the experimental measurement shown on Figure 2, we use the relation J = I S S where I S is the total intensity of current in amperes passing through the section S (A 1 B 1 B 2 A 2 on Figure 8), and S is the surface of the section in m 2 , S = 2c(ab), where 2c is the length of the coil. Thus, we obtain:

G 1 (z) = µ 0 2 • S b a (cosγ 1 ± cos γ 2 )dρ
where G 1 (z) is now in T/A, with G 1 (z) positive if the point P is inside the limits of the coil (-c,c) and negative if P is outside the coil. This equation can be developed using:

tan γ 1 = ρ z + c and, tan γ 2 = ρ |z -c| to obtain G 1 (z) = µ 0 2 • S b a cos(arctan( ρ z + c )) ± cos(arctan( ρ |z -c| )) • dρ G 2 (
z) is found using a translation between the centers of the two coils. total TRM acquired in a 50-µT field for a lava sample from Amsterdam Island [START_REF] Carvallo | Mono Lake or Laschamp geomagnetic event recorded from lava flows in Amsterdam Island (southeastern Indian Ocean)[END_REF]. In the present case, the similarity between the two curves suggests that the NRM is a TRM.

Figure 1 :Figure 2 :

 12 Figure1: Schematic diagram of the C a tVTM where a and b are the first and the second pick-up coil, respectively, connected in series opposition and perpendicular to the vibration axis, c is one of the X-Y field nulling coils, d is the 2-layers mumetal magnetic shield, e is the field solenoid, f is the woolrock insulation and water jacket, g is the sample holder with argon flow, h is the sample and the thermocouple, i is the heating coil, j is the carbon rod, k is the linear motor, and l is the linear bearing. Dashed lines are the connected parts of the same coils. The sketch is out of scale.

Figure 3 :

 3 Figure 3: Sketch of the sample holder. Sample size is 1 cm 3 .

Figure 4 :

 4 Figure 4: Sketch of the linear motor

Figure 5 :

 5 Figure5: Continuous thermal demagnetization curves of NRM and artificial total TRM acquired in a 50-µT field for a lava sample from Amsterdam Island[START_REF] Carvallo | Mono Lake or Laschamp geomagnetic event recorded from lava flows in Amsterdam Island (southeastern Indian Ocean)[END_REF]. In the present case, the similarity between the two curves suggests that the NRM is a TRM.

Figure 6 :Figure 7 : 1 Figure 8 :

 6718 Figure 6: pTRM tail-check: Continuous thermal demagnetization of pTRM(T 1 ,T 2 ) illustrating two different pTRM-tails. The low-temperature pTRM-tail corresponds to the part of pTRM(T 1 ,T 2 ) removed at T2. The hightemperature pTRM-tail corresponds to the part of pTRM(T 1 ,T 2 ) unremoved at T1. A(B) is the low-temperature pTRM-Tail measured at room temperature for PSD(MD) grains. C(D) is the high-temperature pTRM-Tail measured at room temperature for PSD (MD) grains. Figure redrawn from Plenier et al [18].

Table 1

 1

	: Measurements of rate-distortion Harmonics D % D %
	f = 8 Hz	a = 3.61 mm a = 8.14 mm
	1.f	100	100
	2.f	0.451	2.047
	3.f	1.128	1.528
	4.f	0.028	0.044
	5.f	0.085	0.026
	6.f	0.002	0.003
	Total	1.218 %	2.555 %
	Table 2: Signal-to-noise ratio
	Moment Signal-to-noise
	Am2	dB
	1.3e-9	8
	1.9e-9	15
	2.6e-9	21
	6.5e-9	39
	1.3e-8	52
	1.45e-7	97
	3.64e-7	99
	1.46e-6	109


Darchy Victor
email: victor.darchy-ext@enac.fr

Douvenot Rémi
email: remi.douvenot@enac.fr

Jamme Stéphane
email: stephane.jamme@isae-supaero.fr

Galiègue Hélène
email: helene.galiegue@enac.fr

Investigation on Vertical Resolution of Atmospheric Simulations for the Propagation in Realistic Turbulent Media

Keywords: parabolic equation, multiple phase screen, largeeddy simulation, marine atmospheric boundary layer, Kolmogorov spectrum

In split-step methods, the contribution of tropospheric scintillation on the propagating electric field is taken into account by modifying the phase of the signal. The extended-LES method has been proposed to generate phase screens that precisely model atmospheric turbulence. It consists in extending the 3D refractive index fields obtained for the large scales of the medium resolved from atmospheric large eddy simulations (LES) to subrange scales following the K41 theory. In this article, we investigate the impact of the vertical LES resolution on an X-band 2D propagation case. To do so, the extended-LES method is tested in a BOMEX Marine Boundary Layer numerically generated for 3 different vertical meshes. Finally, a log-amplitude variance comparison tends to show that propagation results are dependent on the vertical LES resolution at 10 GHz.

I. INTRODUCTION

The atmospheric boundary layer (ABL) is a dynamic medium where various atmospheric phenomena occur and can potentially degrade the performance of electromagnetic systems. In particular, tropospheric scintillation, that refers to the small scale fluctuations of the atmospheric refractive index, can cause variations in the amplitude and phase of the electromagnetic signal, inducing additionnal losses. An accurate modelling of turbulence is therefore crucial to quantify the impact of the atmosphere on the electromagnetic propagation.

The split-step wavelet (SSW) method is employed to model long-range propagation in complex media based on a 2D discrete "wide angle" parabolic wave equation (PWE) framework [START_REF] Zhou | Modeling the long-range wave propagation by a split-step wavelet method[END_REF] [START_REF] Kuttler | Theoretical description of the parabolic approximation/Fourier split-step method of representing electromagnetic propagation in the troposphere[END_REF]. In this iterative method, the contribution of the atmosphere is modelled through the implementation of a phase screen introducing a phase shift to the propagation signal. This 1D vertical phase screen is orthogonal to the propagation direction and is obtained by integrating the refractive index fields on all the propagation step.

When dealing with a turbulent atmosphere, the phase screens are usually randomly generated as a stochastic process following the classical K41 theory [START_REF] Tatarskii | The Effects of the Turbulent Atmosphere on Wave Propagation[END_REF]- [START_REF] Darchy | Discrete turbulent spectrum modelling for 2D split-step electromagnetic propagation schemes[END_REF]. However, this approach has two main drawbacks. First, it is limited to the modelling of the inertial range of scales of the turbulent flow. Secondly, it assumes a vertical homogeneity of turbulence rarely observed in the ABL.

To better represent the phenomenon at large scales, it has been proposed to generate phase screens from an ABL Large-Eddy Simulation (LES). Due to computational time constraints, the LES grid cannot be as fine as the propagation grid. The subgrid scales are therefore modelled based on the K41 theory. This technique is denoted as extended-LES [START_REF] Gilbert | Electromagnetic wave propagation through simulated atmospheric refractivity fields[END_REF]. In this paper, we test this approach for different vertical LES resolutions on an X-band propagation case.

This article is organized as follows. Section II focuses on the generation of the turbulent atmospheric media. It introduces the studied Barbados Oceanographic and Meteorological EXperiment (BOMEX) LES case and the extended-LES phase screens generation. Finally, Section III highlights the effects of the vertical LES resolution on an X-band 2D propagation case.

II. PHASE SCREENS GENERATION FROM A BOMEX LARGE-EDDY SIMULATION CASE

The objective of this section is to detail the 1D vertical phase screens generation method from LES outputs. The BOMEX LES case settings are initially outlined. The extended-LES generation method is then presented.

A. BOMEX case

The open-source atmospheric Large Eddy Simulation (LES) code MicroHH is used to simulate a BOMEX case of shallow cumulus [7] [8]. This classical case is initialized from data collected during the BOMEX measurement campaign conducted in the tropical Atlantic Ocean between May and July 1969. This simulation generates a convective cloud layer located between 500 and 1500 m of altitude, situated within a moist turbulent Marine Atmospheric Boundary Layer (MABL).

The time evolution of a 5 × 5 × 3 km LES domain is studied with a 50 × 50 m horizontal resolution over 6 hours for 3 different vertical resolutions d z : 50 m, 25 m and 10 m. The simulations become statistically steady after

∆n = n -n x,y , (1)
where • x,y indicates the mean operator over an horizontal plane. This figure highlights the vertical inhomogeneity of turbulence in the low troposphere. Details on boundary and initial conditions can be found in the comparative study led by Siebesma et al. [START_REF] Siebesma | A large eddy simulation intercomparison study of shallow cumulus convection[END_REF]. The main LES simulation parameters used for this study are summarized in Table I.

B. Phase screens computation

The purpose of this section is to generate 1D vertical phase screens from the 3D thermodynamic data obtained from the LES at each time step. To do so, the extended-LES phase screens generation method developed by Gilbert et al. [START_REF] Gilbert | Electromagnetic wave propagation through simulated atmospheric refractivity fields[END_REF] is used. It consists in considering the turbulent phase Φ at an arbitrary y location of the LES domain as the combination of two contributions so that

Φ = Φ r + Φ s . (2)
Φ r indicates the turbulence induced by the resolved scales of the turbulent flow corresponding, for 1D vertical phase screens, to scales down to the LES resolution d z . It is computed by directly integrating the fluctuations of the radiofrequency refractive index fields obtained from LES outputs on a propagation step ∆ x . Φ s denotes the effect of the subgrid scales and enables to extend the phase spectrum to finer scales than the resolution imposed by LES. Its generation follows the same principle as the classical stochastic method [START_REF] Fabbro | Comparison of 2D and 3D electromagnetic approaches to predict tropospheric turbulence effects in clear sky conditions[END_REF]. It is here modelled through the following 1D phase spectrum [3]

S 1D Φs (κ z ) = 2πk 2 0 ∆ x 0.055 C 2 n κ -8/3 z , if κ z ∈ [ 2π dz , 2π ∆z ], 0
otherwise.

(3) In this formula, k 0 is the electromagnetic wavenumber while κ z indicates the vertical component of the wavenumber of the turbulent flow. C 2 n corresponds to the structure constant of the refractive index that quantifies the scintillation intensity. Its value generally ranges from 10 -18 m -2/3 in a case of weak turbulence to 10 -12 m -2/3 under strong turbulence. It can be obtained at every point of the LES domain as proposed by Wilson and Fedorovich [START_REF] Wilson | Direct evaluation of refractive-index structure functions from large-eddy simulation output for atmospheric convective boundary layers[END_REF]. Here, C 2 n is fixed by computing the average value over the entire domain for the last 3 hours of simulation.

To illustrate the method, Figure 2 gives a schematic representation of the resulting total turbulent phase sprectrum. ∆ z corresponds to the vertical resolution of the propagation grid. Figure 3 shows the impact of the modelled subgrid scales on the resulting phase screen. The extended-LES turbulent phase presents finer fluctuations.

The effect of turbulence is cumulative with range. Its impact becomes significant only from several tens of kilometers for X-band frequencies. Since periodic boundary conditions are applied in both the x and y homogeneous directions, ergodicity hypothesis allows to successively use LES results recorded at uncorrelated times to virtually extend the computational domain [START_REF] Gilbert | Electromagnetic wave propagation through simulated atmospheric refractivity fields[END_REF]. Phase screens from each domain are computed following the presented method and stored. The propagation is then performed by randomly selecting the screens to avoid any correlation phenomena.

III. X-BAND PROPAGATION SIMULATION

The aim of this section is to test the presented extended-LES method for different vertical LES resolutions on an X-Band propagation case. The numerical results are analyzed through a log-amplitude variance comparison.

A. Propagation model and configuration

The propagation simulations are here performed with the SSW-2D software [START_REF] Douvenot | SSW-2D : Split-step wavelet 2D. Software[END_REF]. The propagation of a wave front at a 10 GHz frequency over a range of L x = 95 km is investigated. The computational domain has a vertical size L z of 3 km and a vertical discretization step ∆ z = 5λ/3 = 0.05 m. The value of C 2 n is computed from the atmospheric simulations as proposed in Section II giving the value 0.23×10 -12 m -2/3 . The propagation step ∆ x is fixed at 250 m which leads to the computation of N s = 380 phase screens following the generation method presented in Section II. A complex source point of size 6λ is located at 10 m from the first phase screen at an altitude z s = L z /2 = 1500 m. The main propagation parameters are summarized in Table II.

B. Numerical results for various vertical LES resolutions

To highlight the impact of the vertical LES resolution on the electromagnetic propagation, the log-amplitude is computed at every point of the propagation domain as

χ(x, z) = ln |E t (x, z)| |E 0 (x, z)| . (4)
E t indicates the electric field propagated in the modelled turbulent medium while E 0 is the field obtained in a turbulentfree atmosphere. The log-amplitude variance σ 2 χ quantifies the impact of turbulence on the electromagnetic propagation. It is numerically obtained at every range x by computing the variance of the vertical log-amplitude vector. σ 2 χ profiles obtained for different d z are plotted in Figure 4. The solid lines indicate a propagation through an atmosphere generated with the extended-LES method while the dashed lines denote the propagation through the non-extended atmospheric domain. These results are averaged over 500 Monte-Carlo runs.

Figure 4 first illustrates the cumulative effect of atmospheric turbulence with range. The log-amplitude variance computed using the extended-LES method reaches a value greater than 1 dB 2 after 95 km of propagation, which indicates a strong scattering [START_REF] Wheelon | Electromagnetic scintillation: volume 2, weak scattering[END_REF]. This figure also shows a significant difference between the propagation results in a non-extended LES atmosphere (dashed lines) according to the imposed resolution. These differences are due to the non-inclusion of scales that are not directly resolved by the atmospheric simulations. Moreover, these plots indicate that this gap is reduced by the use of the extended-LES method (solid lines). Indeed, simulations carried out for 3 different vertical resolutions seem to provide equivalent σ 2 χ profiles up to 30 km of range. However, these results diverge from the one obtained for the finest resolution (d z = 10 m) from 35 km for d z = 50 m and from 60 km for d z = 25 m. Thus, Figure 4 confirms the need to take into account the subgrid scales of the flow to accurately model propagation in turbulent media. These plots also point out that the accuracy of the extended-LES method depends on the vertical resolution used for generating the atmosphere. When using the extended-LES method, it appears that insufficient LES resolution can lead to an underestimation of the impact of turbulence. In this paper, the extended-LES method has been tested on an X-Band 95 km range propagation case through a shallow cumulus MABL. The impact of the vertical LES resolution on the resulting electromagnetic propagation has been quantified by computing the log-amplitude variance for 3 different vertical meshes. The results appear to be dependent on this parameter at 10 GHz. These initial conclusions suggest the need for further investigations at finer LES resolutions as well as at higher signal frequencies.

IV. CONCLUSION

Fig. 1 .

 1 Fig. 1. 2D Refractive index fluctuations field from a dz = 10 m BOMEX LES case.

Fig. 2 .

 2 Fig. 2. Extended-LES vertical phase spectrum representation.

Fig. 3 .

 3 Fig. 3. Comparison between a resolved scales phase screen (a) and the corresponding extended-LES phase screen (b) both obtained from a BOMEX LES case with dz = 10 m.

Frequency : f 10

 10 GHz Propagation domain size : Lx × Lz 95 × 3 km Propagation step : ∆x 250 m Number of screens : Ns 380 Propagation vertical resolution : ∆z 5λ/3 = 0.05 m C 2 n 0.23 × 10 -12 m -2/3

Fig. 4 .

 4 Fig. 4. Extended-LES (solid lines) vs LES (dashed lines) log-amplitude variance comparison for dz = 50, 25, 10 m.

ACKNOWLEDGMENT

The authors would like to thank the Centre National d' Études Spatiales (CNES) and the Fédération de recherche ENAC ISAE-SUPAERO ONERA (FONISEN) for the funding. This work was performed using HPC resources from CALMIP (Grant 2023-P22033).


domain: spi phys

text:
Gaëtan Antoine
email: gaetan.antoine@isae-supaero.fr

Bruno De Araujo

Romain Pascaud

Christophe Morlaas

Alexandre Chabory

Vincent Laquerbe

Gautier Mazingue

Comparison of Effective Permittivity Retrieval Methods of 3D-Printed Unit Cells for Dielectric Resonator Antenna Applications

Keywords: DRA, 3D printing, permittivity, unit cell

come

I. INTRODUCTION

Over the last 20 years, three-dimensional (3D) printing technologies have been developed to manufacture complex structures. Regarding antenna applications, horn [START_REF] Chieh | Development of a Ku-band corrugated conical horn using 3-D print technology[END_REF], lens [START_REF] Yi | 3-D printed millimeter-wave and terahertz lenses with fixed and frequency scanned beam[END_REF], reflector [START_REF] Nayeri | 3D printed dielectric reflectarrays: Lowcost high-gain antennas at sub-millimeter waves[END_REF], and more recently dielectric resonator antennas (DRAs) [START_REF] Xia | 3-D-printed wideband multi-ring dielectric resonator antenna[END_REF]- [START_REF] Xia | 3-D-printed wideband circularly polarized dielectric resonator antenna with two printing materials[END_REF], to cite a few, have already been 3D-printed.

DRAs have been studied since the 1980s [START_REF] Long | The resonant cylindrical dielectric cavity antenna[END_REF] and are still investigated today. This is mainly due to their attractive features such as compact size, high gain, and high radiation efficiency [START_REF] Petosa | Dielectric Resonator Antenna Handbook[END_REF]. 3D-printing of dielectric materials now opens up new possibilities in DRA design. Thus, the control of the shape [START_REF] Basile | Design and manufacturing of super-shaped dielectric resonator antennas for 5G applications using stereolithography[END_REF] [START_REF] Chietera | Dielectric resonators antennas potential unleashed by 3D printing technology: A practical application in the IoT framework[END_REF] and dielectric constant (structuring layers [START_REF] Xia | 3-D-printed wideband circularly polarized dielectric resonator antenna with two printing materials[END_REF] or unit cells [START_REF] Xia | 3-D-printed wideband multi-ring dielectric resonator antenna[END_REF]) of the DRA have been explored in the last 5 years. Most of the early works were using materials with dielectric constant up to 10. With new 3D-printing technologies being developed, materials with higher relative permittivity and low loss tangent can now be 3D-printed, as for example ceramics [START_REF] Mazingue | 3D printed ceramic lowprofile GNSS antenna for SmallSats[END_REF]. In particular, a DRA made up of periodic subwavelength unit cells mixing zirconia (ε r = 32.5 and tan δ = 1.9 × 10 -4 at 10 GHz) and air has been studied in [START_REF] Morales | 3D-printed ceramics with engineered anisotropy for dielectric resonator antenna applications[END_REF]. The filling ratio and the anisotropy of the unit cells allowed to control the resonance frequency and the polarization of the DRA.

In order to design a DRA made up of periodic unit cells, a full-wave analysis can be performed. Most of the time, it is costly in time and computer resources. As it is well known that the behavior of a DRA is mainly determined by its dimensions and effective permittivity, an alternative solution is to analyze a single unit cell of the DRA in order to retrieve its effective permittivity. From a transcendental equation or by comparing with a solid model, this permittivity can be linked to the resonance frequency of the 3D-printed DRA.

Methods to relate the filling ratio of a unit cell, i.e. the percentage of dielectric in the unit cell volume, to its effective permittivity have been studied in previous works using either Maxwell Garnett approximation (MGA) [START_REF] Mazingue | 3D printed ceramic antennas for space applications[END_REF][13], Sparameters method [START_REF] Mazingue | 3D printed ceramic antennas for space applications[END_REF][14], or plane wave expansion method (PWEM) [START_REF] Hehenberger | Effective permittivity measurement of 3D-printed dielectric crystals[END_REF]. However, the size of the unit cells has always been chosen so that they are electrically-small to ensure a homogeneous macroscopic behavior. For high frequency applications, mechanical limits of 3D printer may no longer allow to comply with this assumption. Therefore, it appears interesting to evaluate the effective permittivity obtained with these methods for larger unit cells.

In this work, the behavior of the effective permittivity for different unit cells is studied using the three methods mentionned before. In Section II, unit cell topologies are presented. In Section III, the methods to retrieve the effective permittivity from simulations are detailed. In Section IV, numerical analyses are performed and the computing time and accuracy of the methods are compared. Finally, Section V concludes the paper.

II. UNIT CELL TOPOLOGY

For 3D-printed antennas made up of periodic subwavelength unit cells, cubic cells are the most commonly used among other types of Bravais lattices (tetragonal, monoclinic, orthorhombic, rhombohedral, hexagonal, or triclinic). It can be explained by the feasibility of such structures, the cubic being the easiest to manufacture.

Here, three common topologies for cubic unit cell are compared: simple cubic (SC), body-centered cubic (BCC), and face-centered cubic (FCC). They can be seen in Fig. 1, with the 3D-printed dielectric in blue and its spherical inclusions of air. A complementary topology can exist, reversing dielectric and air, but it requires additional dielectric cylinders between the dielectric spheres in that case in order to be printable.

Finally, the topologies are defined as follows: • The SC lattice consists of one lattice point on each corner of the cube. Each atom at a lattice point is then shared equally between eight adjacent cubes, and the unit cell, therefore, contains a total of one atom (i.e., air inclusion). • The BCC lattice has one lattice point in the center of the unit cell in addition to the corner points. It contains a total of two atoms per unit cell. • The FCC lattice has lattice points on the faces of the cube in addition to the corner points. It contains a total of four atoms per unit cell.

III. PERMITTIVITY RETRIEVAL METHODS

In order to manufacture 3D-printed DRAs made up of periodic unit cells, dielectric behavior of the 3D-printed material has to be controlled. To find the resonance frequency associated with the fundamental mode of a dielectric resonator (DR), the most accurate way is to do an eigenmode analysis. For a DR made up of unit cells, it is also the most accurate method and it will be considered as our reference. However, for a DR with a large amount of unit cells, this method requires a lot of computing time and resources.

In order to reduce such requirements, a single unit cell of the DRA is often analyzed. The objective of this study is thus to compare three methods allowing to extract the effective permittivity of such unit cells. Note that all the simulations described in the following section are done with Ansys HFSS software.

A. Reference: Eigenmode analysis of a 3D-printed DR At first, a DR over an infinite ground plane is considered for eigenmode analysis. Its permittivity is directly assigned on Ansys HFSS and it is lower than the one of the 3Dprinted material. It is referred to as solid model and used as a reference as presented in Fig. 2a. More specifically, an eigenmode analysis of this solid model is performed to find the resonance frequency of its fundamental mode.

In a second step, an eigenmode analysis of the DR made up of periodic unit cells is performed. A DR made up of 64 BCC unit cells is presented in Fig. 2b as an example. The resonance frequency of its fundamental mode is identified by analyzing the electric and magnetic fields inside the DR. By varying the filling ratio of the unit cells, and thus, the effective permittivity of the DR, we set the computed resonance frequency of this mode at the same frequency as for the solid model. This method can be time-consuming. However, it allows to design with a good accuracy a DR made up of periodic unit cells that has the same fundamental mode resonance frequency as a solid DR.

B. Method 1: Maxwell Garnett approximation

The theory behind the Maxwell Garnett approximation (MGA) is described in detail in [START_REF] Choy | Effective Medium Theory: Principles and Applications[END_REF]. The MGA describes the effective relative permittivity ε ef f of a dielectric material with relative permittivity ε m filled with dielectric spheres with relative permittivity ε i as

ε ef f = ε m 2δ i (ε i -ε m ) + ε i + 2ε m 2ε m + ε i -δ i (ε i -ε m ) (1)
where δ i is the volume fraction of the inclusions. In our case, ε m refers to the ceramic dielectric constant, while ε i = 1 for the air inclusions. Note that MGA method is expected to be valid at low volume fractions δ i .

C. Method 2: S-Parameters method

This method relies on the Floquet simulation of a single cubic unit cell of side a uc . As shown in Fig. 3, the simulation setup consists of a single cubic unit cell and a vacuum box four times the size of the unit cell along the z-direction. The vacuum box has periodic boundary conditions on the four side walls along the xand y-axis and Floquet ports on the two faces along the z-axis. It provides a normal incident plane wave which polarization can be set.

The simulated S-parameters are used to retrieve the effective permittivity of the unit cell as a function of the frequency using the method described in [START_REF] Numan | Extraction of material parameters for metamaterials using a full-wave simulator [education column[END_REF]. In fact, the refractive index n and the impedance z of the unit cell are given by

z = (1 + S 11 ) 2 -S 2 21 (1 -S 11 ) 2 -S 2 21 , (2)
n = 1 k 0 a uc cos -1 1 2S 21 (1 -S 2 11 + S 2 21 ) . (3)
The effective relative permittivity is finally computed using the ratio between the refractive index and the impedance of the unit cell

ε ef f = n z . (4)
This method is fast and easy to implement but it relies on several assumptions that restrict its applicability. One of them being that the unit cell must remain electrically-small. Commonly, effective permittivity retrieval of subwavelength unit cells is done with this method [11][12].

D. Method 3: Plane wave expansion method

The theory behind the plane wave expansion method (PWEM) is described in [START_REF] Rumpf | Synthesis of spatially variant lattices[END_REF]. Only the eigenmode simulation setup and the exploitation of the data will be detailed hereafter.

Here, the simulation setup consists of a single unit cell with periodic boundary conditions along the x-, y-, and z-axis, as shown in Fig. 4.

The eigenmode simulation provides the frequencies associated with each eigenmode. Plotting the dispersion diagram, i.e., the frequency as a function of the wavenumber k, we can retrieve the effective relative permittivity with the following formula

ε ef f = kλ 0 2π 2 , ( 5
)
where λ 0 is the free-space wavelength.

In this work, simulations will be done for the Γ-X part of the Γ-X-M-Γ-R-X irreductible Brillouin zone (IBZ) used for cubic cells (see Fig. 4). This path is here considered because it is the same as for method 2. For future work, it can be interesting to study the anisotropy of the unit cell by analyzing other parts of the IBZ.

IV. RESULTS

A. Optimization of filling ratio of the unit cells

A cubic DR of side equal to 16 mm, placed on an infinite ground plane, and using a dielectric with ε r = 15 and tan δ = 1.9 × 10 -4 has first been designed. The resonance frequency f 0 of its first two modes, TE x δ11 and TE y 1δ1 , is found to be f 0 = 3.192 GHz, i.e., a wavelength λ in the dielectric equal to 24.27 mm.

Three DR of the same dimensions, but made up of different number of unit cells, are then simulated with the eigenmode analysis described in Sec.III.A. They are structured as shown in Table I and Fig. 5 (BCC unit cells are used in this figure), considering 1 to 64 unit cells. However, considering Fig. 5c, one may wonder whether a single 16 mm unit cell can still be analyzed as a cubic DR. Each DR has been optimized using SC, BCC, and FCC unit cells. As shown in Table II, the filling ratio of the unit cells has been chosen to obtain the same resonance frequency as for the reference solid model, namely 3.192 GHz.

One can observe that the filling ratio of the unit cells varies according to their topology, namely SC, BCC, or FCC. This is probably due to the spatial dispersion caused by their different degrees of isotropy. Moreover, as the size of a unit cell a uc increases, its filling ratio also varies whatever the type of unit cell.

Although the resonance frequency is fixed to the required value of 3.192 GHz, it should be noted that it is not sufficient to guarantee the same antenna behavior as the reference DRA at the end. For instance, another relevant parameter of the antenna performance of a DRA is its quality factor Q. As can be seen in Fig. 6, in the case of DRs made up of periodic unit cells, the value of Q varies with the size and topology of the unit cell. Therefore, the same electromagnetic properties as for the reference antenna cannot be assumed. Further full wave simulation of complete DRA are thus required.

B. Computing time

Simulations of DRs made up of unit cells have been performed on a PC with 256 Go of RAM and a processor Intel®Xeon®CPU E5-2690.

Considering for instance a DR made up of 64 BCC unit cells, i.e., a uc = 4 mm, the computing time necessary to obtain the resonance frequencies of the DR and to run method 2 and method 3 are given in Table III.

From Table III, it can be seen that the eigenmode simulation of a full DR made up of unit cells can be extremely timeconsuming. It shows why it can be useful to employ method 2 or method 3 in order to derive the permittivity of the DR through the permittivity of a single unit cell.

C. Comparison between methods

The results of the simulations are reported in Table IV. The effective permittivity ε ef f of the unit cell at the resonance frequency f 0 = 3.192 GHz of the DR, and the error in percentage from the reference value, namely ε r = 15, are given for each method, cell topology and size.

Values for the computed effective permittivity are very close between method 2 and method 3 for BCC and FCC unit cells when a uc is equal to 4 mm. The smallest SC unit cell indicates a trend: method 2 does not give accurate results for every unit cell topology. This trend is confirmed with 8 mm unit cells, where the permittivity of the SC topology at the resonance frequency of the DR is equal to 19.07 instead of being close to 15 as it is the case in method 3.

When a uc is too close to the wavelength, dispersion effects may cause the results of method 2 to be unreliable regardless of the topology of the unit cell. The results of method 2 are then noted N/A (not applicable). On the contrary, method 3 is still efficient for unit cell dimension around the wavelength.

To conclude, the PWEM remains the most accurate method whatever the unit cell size in our example.

V. CONCLUSION

In this paper, a comparison between three methods to extract the effective permittivity of periodic unit cells for 3D-printed DRA has been proposed.

When designing a DRA made up of periodic unit cells mixing dielectric and air, extracting its effective permittivity

Fig. 1 .

 1 Fig. 1. From left to right: SC, BCC, and FCC cell.

Fig. 2 .

 2 Fig. 2. Simulation setup for the eigenmode analysis for the (2a) solid and for instance, (2b) 64 BCC cells models.

Fig. 3 .

 3 Fig. 3. Floquet simulation setup of a single BCC cell.

Fig. 4 .

 4 Fig. 4. Eigenmode simulation setup of a single BCC cell and its irreductible Brillouin zone (IBZ)

Fig. 5 .

 5 Fig. 5. DR made up of BCC unit cells.

Fig. 6 .

 6 Fig. 6. Q of the DR as a function of auc for the different topologies of unit cells.

TABLE I UNIT

 I CELLS SPECIFICATIONS FOR EACH DR AT f 0 = 3.192 GHZ

	Number	auc	auc/λ
	64	4 mm	1/6
	8	8 mm	1/3
	1	16 mm	1/1.5

TABLE II UNIT

 II CELLS FILLING RATIO AT f 0 = 3.192 GHZ

			filling ratio
	auc	SC	BCC	FCC
	4 mm	47.11%	47.48%	47.54%
	8 mm	50.21%	48.51%	48.42%
	16 mm 48.81%	54.58%	52.37%

TABLE III COMPARISON

 III OF COMPUTING TIME BETWEEN THE DIFFERENT METHODS FOR A 4 MM BCC CELL

	Type of method	Computing time (hh:mm:ss)
	Reference	03:07:26
	Method 2	00:02:46
	Method 3	00:06:17

This paper's copyright is held by the author(s). It is published in these proceedings and included in any archive such as IEEE Xplore under the license granted by the "Agreement Granting EurAAP Rights Related to Publication of Scholarly Work." Authorized licensed use limited to: ISAE. Downloaded on June 16,2023 at 09:25:49 UTC from IEEE Xplore. Restrictions apply.

Authorized licensed use limited to: ISAE. Downloaded on June 16,2023 at 09:25:49 UTC from IEEE Xplore. Restrictions apply.

ACKNOWLEDGMENT

The authors would like to thank the French Agency for Innovation of Defense (Agence de l'Innovation de Défense, AID) and the French National Center of Space Studies (Centre Nationale d'Etudes Spatiales, CNES) for the funding.

can be very time-consuming depending on the number of cells. For this reason, three methods allowing this extraction using a single unit cell have been compared. The S-parameters method is commonly used in the literature but has assumptions about the size of the unit cell that restricts its applicability for high frequency applications. However, the PWEM remains accurate for larger unit cells, even close to the wavelength, and no matter the cell topology.

Future works include studying the results provided by these methods on the spatial and frequency dispersion behavior of cubic unit cells.