E Baudin

M E Hayden

G Tastevin

P J Nacher

Nonlinear NMR dynamics in hyperpolarized liquid 3 He

In a highly polarized liquid (laser-polarized 3 He-4 He mixtures in our experiment), dipolar magnetic interactions within the liquid introduce a significant nonlinear and nonlocal contribution to the Bloch equation that leads to instabilities during NMR evolution. We have launched a study of these instabilities using spin echo techniques. At high magnetizations, a simple 180 • rf pulse fails to refocus magnetization, so we use a standard solid-state NMR pulse sequence: the magic sandwich. We report an experimental and numerical investigation of the effect of this sequence on unstable NMR evolution. Using a series of repeated magic sandwich sequences, the transverse magnetization lifetime can be increased by up to three orders of magnitude.

Résumé

Dans un liquide fortement polarisé (ici, une solution d' 3 He hyperpolarisé dans l' 4 He liquide), les interactions dipolaires magnétiques apportent une importante contribution non linéaire et non locale à l'équation de Bloch qui rend l'évolution RMN instable. Nous avons entrepris une étude de cette instabilité par des techniques d'écho de spin. A forte aimantation, une impulsion rf de 180 • appliquée lors de l'évolution ne refocalise pas l'aimantation, aussi utilisons nous une technique standard de RMN du solide : le "magic sandwich". Nous décrivons les premiers résultats d'une étude expérimentale et numérique des effets de cette séquence sur une évolution RMN instable. Une augmentation de la durée de vie de l'aimantation transverse de 3 ordres de grandeur est obtenue en appliquant périodiquement des "magic sandwich".

INTRODUCTION

Conventional treatments of nuclear magnetic resonance (NMR) ignore collective effects for spin dynamics, and thus the evolution of the local magnetization of a single spin species is described by simple (linear) Bloch equations. This approach is usually justified in liquids, where short-range dipolar spin interactions (averaged out by fast atomic diffusion) only lead to relaxation. But it is actually limited to conventional implementations of NMR where equilibrium magnetizations are very small. The emergence of high-resolution NMR spectroscopy as a powerful tool for condensed matter physics, analytical chemistry, structural biology, and medicine, has motivated considerable efforts to improve signal-to-noise ratios through increased use of high static fields and of probes with high quality factors. These approaches enhance nonlinear contributions to spin dynamics that arise from radiation damping [1]. A more complex situation arises when additional nonlinear contributions of long-range dipolar interactions between the local magnetization and that of the remainder of the sample are considered. The influence of these dipolar interactions becomes significant at large magnetization densities and is thus expected to play an increasing role as technological advances lead to the use of higher and more uniform static fields. In particular, they are expected to lead to spin turbulence and spatiotemporal chaotic behavior [2,3].

Spectacular dynamical effects resulting from distant dipolar fields (DDF) have been observed in liquids obtained by condensation of hyperpolarized 3 He (i.e. having a high, out-of-equilibrium nuclear polarization produced by laser optical pumping). These effects range from spectral clustering (at small tipping angles, with the spontaneous appearance of long-lived geometry-dependent modes of coherent nuclear precession [4,5,6,7]), to precession instabilities (at large tipping angles, with exponential growth of spatially inhomogeneous magnetization patterns [8,9]). These effects have been accounted for, at least qualitatively, by theoretical approaches and numerical studies [10]. Experimental investigations of nonlinear DDF effects are challenging with thermally polarized liquids; the characteristic dipolar field for 1 H at 12 T is 43 nT. In contrast, dense laser-polarized systems such as liquefied 3 He are ideal model systems for the study of nonlinear NMR dynamics. 3 He has a high nuclear magnetic moment, and its liquid phase provides dipolar fields as large as 100 µT for pure 3 He polarized to 45%.

We have recently succeeded in generating spin echoes in highly polarized 3 He-4 He dilute mixtures (with dipolar fields of order 1 µT), and have demonstrated that it is possible to dynamically stabilize magnetization distributions in intrinsically unstable regimes [11]. Here, we report additional experimental results, provide a discussion of echo decays both at low and high DDF, and present a comparison to results of exploratory numerical simulations.

EXPERIMENTAL SETUP Preparation of hyperpolarized helium samples

The 3 He nuclei in our experiments are hyperpolarized by optical pumping, which amounts to the transfer of angular momentum from an intense circularly polarized laser beam to atoms by light absorption. 3 He gas is continuously injected at low pressure (mbar) into the room temperature part of a glass cell where it is optically pumped in a few seconds to 40% polarization by metastability exchange optical pumping [12], using a 2 W fibre laser at 1083 nm [13]. The polarized gas flows down a narrow tube into the low temperature (1 K) part of the apparatus (a 0.44 cm 3 spheroidal Pyrex cell) partially filled with superfluid liquid 4 He. Wall relaxation is inhibited through the use of Cs coating on the glass, and bulk dipole relaxation is of order several hours at the low 3 He molar fractions used in our experiments (typically 1-5%) [14]. We control the volume of liquid admitted to the cell (typically half full) and the nuclear polarization in the liquid phase (up to 40%). We can thus choose the initial magnetization M , and the 3 He spin diffusion coefficient D (in the 10 -3 -10 -2 cm 2 /s range, controlled both by the temperature T and the 3 He molar fraction X) for each experiment. At large tipping angles, NMR excitation induces a significant irreversible loss of this out-of-equilibrium magnetization, so that a new batch of polarized gas has to be condensed for each experiment.

NMR setup

NMR is performed at very low magnetic field (2.3 mT, shimmed to 20 ppm over the sample) using crossed coils. The frequency of the rf pulses (74.5 kHz) is set to lie within 1 Hz of the Larmor frequency corresponding to the average of the external field B ext over the sample. The magnetic field is produced over the entire apparatus by a series of seven coils. Slow fluctuations of B ext due to remote magnetic perturbations are measured with a fluxgate magnetometer, and compensated using additional bucking coils. Overall field fluctuations are thus reduced by more than a factor of 10, and the corresponding changes in Larmor frequency are found to be smaller than 0.5 Hz over the duration of experiments (several minutes).

The rf transmit coil is actively shielded [15] to prevent eddy currents from significantly heating the low temperature parts of the apparatus when an oscillating rf field B rf (≃ 30 µT) is applied. This also ensures that eddy currents do not distort the applied rf field, that is uniform to within 200 ppm over the volume of the sample cell. The receive coil is tuned to the Larmor frequency with a very low Q factor (≃ 2) so that radiation damping effects are negligible. In-house PC-controlled analog hardware allows flexible control of hard pulse sequences.

NMR EVOLUTION IN A MAGNETIZED LIQUID

We use the following modified Bloch equation to describe the NMR evolution of local magnetization M in the frame rotating at the rf angular frequency:

d M dt = γ( B eff + B DDF ) × M + D∆ M , (1)
where γ is the 3 He gyromagnetic ratio (γ/2π = 32.4 MHz/T) and B eff is the effective field resulting from B ext and B rf . B DDF is the distant dipolar field induced by the distribution of magnetization. The conditions of our experiments are such that contributions to the Bloch equation from relaxation and radiation damping can be neglected. Using a secular approximation in the rotating frame, the distant dipolar field can be expressed as:

B DDF ( r) = µ 0 8π d r ′ 2 M z ( r ′ ) -M ⊥ ( r ′ ) | r -r ′ | 3 (3 cos 2 θ -1), ( 2
)
where θ is the angle between r -r ′ and the average direction ẑ of B ext . B DDF brings a nonlinear nonlocal contribution to the Bloch equation (Eq. 1).

The main features of the dynamical evolution of magnetization can be discussed in terms of three parameters: γ∆B ext , the rate at which the external field inhomogeneities defocus magnetization (∆B ext being a measure of the typical variation of B ext over the sample); γB dip , a characteristic rate associated with the initial DDF (B dip =µ 0 M is the amplitude of the dipolar magnetic field created by homogeneous magnetization M before the first rf pulse of the sequence); D/r 2 , the diffusion-induced damping rate for a magnetization pattern of typical scale r. It is important to note that B DDF is the field created by the instantaneous distribution of local magnetization and evolves in both space and time, whereas B dip is fixed before each NMR experiment. 3 He-4 He mixtures. At low magnetization density (negligible dipolar interactions), we observe a characteristic FID half-life of 0.3 s, a value that is set by the homogeneity of the magnetic field over the cell (Fig. 1, trace a). In contrast, an experiment performed with B dip = 0.9µT reveals a dramatic collapse of the FID signal on a timescale of 70 ms (Fig. 1, trace b). The deviation of the NMR signal S from its initial value S 0 grows exponentially with time (inset to Fig. 1).

A model developed by J. Jeener [5] and numerical simulations described in section [3,8,16] have shed light on magnetization dynamics following a large tipping angle pulse: at large B dip , during transverse precession, unstable inhomogeneous magnetization patterns develop and grow exponentially with time. These patterns can arise from an initial seed of inhomogeneity following the tipping pulse, from stray or applied static field gradients, or from edge effects that depend on the sample shape. Jeener's model for infinite media predicts that the instability growth rate Γ S of the deviation S 0 -S is proportional to B dip . This analytical model and corresponding numerical simulations yield the value 2πΓ S /γB dip = 2 √ 2π/3 ≃ 3.0. Simulations of finite-sized samples and previous experiments on spheroidal samples agree with one another, but yield a significantly lower ratio 2πΓ S /γB dip ≃ 2.3 [8,16]. The present experiments confirm the latter result.

TIME-REVERSAL EXPERIMENTS USING MAGIC SANDWICH SEQUENCES

Our objective is to characterize the onset and growth of DDF-induced inhomogeneous magnetization patterns, beyond the mere measurement of a global growth rate. In particular, we are interested in the determination of spatial length scales associated with magnetization patterns that develop during unstable NMR precession. To this aim, we plan to use echo techniques and the established link between diffusion-induced NMR relaxation rates and dominant frequencies of spatial modulations of the magnetization density. The conventional Hahn echo sequence (90 • -τ -180 • ) effectively generates a spin echo in weakly magnetized systems, such as laser polarized gas or dilute liquid mixtures as long as the nuclear polarization is kept very low. In contrast, this sequence fails completely when used for experiments on highly magnetized mixtures [8,11]. This failure is due to the nature of the dipolar interactions; the 180 • pulse reverses both M and B DDF , leaving the quadratic term in Eq. ( 1) unchanged.

The Magic Sandwich sequence (MS) was introduced during the 1970s in the context of solid-state NMR to compensate for the strong effect of local dipolar fields [17]. In our experiment it has been adapted and used, for the first time, on a liquid system. It consists of continuously applied and appropriately phased rf pulses during a given time interval. It mixes longitudinal and transverse components of M , that contribute in a different way to the dipolar field (Eq. ( 2)), and causes a time-reversed evolution for the magnetization at half the pace of the free evolution while it is applied. Applying quasi-cw rf for up to several seconds without significantly heating the system is made possible through the use of shielded coils and by operating at very low NMR frequency.

Effect of rf detuning

The magic sandwich sequence introduced by Rhim, Pines and Waugh [17] can be summarized as:

90 • y , (180 • x , -180 • x ) n , -90 • y , (3)
where the 90 • and 180 • rotations are performed along orthogonal directions ŷ and x, and the (180 • x , -180 • x ) pulse pair is repeated n times. This scheme only works for perfectly resonant applied rf field. A detuning ∆ω between the Larmor and rf angular frequencies introduces a small angular tilt θ = ∆ω/γB rf between the desired rotation axis ( B rf ) and the effective rotation axis ( B eff ). The net result is that a (180 • ,-180 • ) pair will rotate a magnetization vector that is initially aligned with the ẑ-axis by an angle 4θ (Fig. 2). 2: Angular error associated with rf detuning. In the rotating frame, the effective rotation axes are tilted from the x-axis by an angle θ, proportional to the rf frequency offset. After the first 180 • pulse, the Bloch vector lies 2θ away from the targeted ẑ-axis, and the ( 180 • x , -180 • x ) pulse pair leads to a net 4θ angular tilt. Here, BDDF is assumed to have a negligible effect.

B ¡¡ ¢ B ¡¡ £ 0 180° (180°,-180°) 4θ 2θ θ θ FIG.
In our experiments, a typical MS sequence contains about 100 such rotation pairs and lasts 10-100 ms (rf intensity is such that a 90 • pulse lasts 0.5 ms). In order to avoid significant cumulated error (θ ≪ π/1000), one would need a small rf detuning (∆ω ≪ 10 rad/s). This cannot be achieved in our experiment due to the limited stability of the applied field B ext . The sequence becomes robust against rf detuning when integer number of turns are performed, e.g. using ( 360 • x , -360 • x ) pulse pairs, and the direct (1 st order) angular errors described above are eliminated. The remaining cumulative 4 th order angular error has a negligible effect for ∆ω = 10 rad/s. This modification of the basic MS sequence (3) is necessary to ensure that it correctly reverses the effect of dipolar fields while the rf is applied.

Experimental results

Figure 3 displays the data obtained using a 90 • -τ -MS sequence with MS duration 2τ , performed for various magnetization densities. Magnetization freely evolves for a time τ during which the FID signal is recorded. Then a MS sequence is applied to unwind the magnetization back to time t=0, so that the subsequent free evolution is expected to yield an NMR echo signal identical to the initial FID signal. This type of behavior is observed at the lowest magnetization shown in Fig. 3. At the highest magnetizations, echoes are still obtained long after the initial FID signals have vanished. However, the echo shapes progressively change, and their relative amplitudes decrease almost linearly with increasing B dip (inset to Fig. 3). This decrease hardly depends on the diffusion coefficient, which roughly scales with the inverse of the 3 He molar fraction over the investigated range. Imperfect refocussing can be attributed to sequence imperfections (in rf timing, phase or amplitude). Moreover, during the application of the MS sequence, while magnetization patterns that have developed during the free evolution period are driven back to their initial state, complementary families of patterns are expected to become unstable and grow [5].

The MS duration in the 90 • -τ -MS sequence can be set arbitrarily to achieve time reversal and bring the dipolarcoupled system back to various times. In particular, when the MS duration is set to 4τ (magnetization unwound back to t = -τ ) one observes an echo signal that peaks at time t = 6τ (i.e., a maximal refocussing of magnetization at time τ after the rf pulses are applied). As recently reported, echo trains can be obtained using a series of such MS sequences, periodically repeated [11]. In the following section, we focus on a dedicated echo sequence, analyze the conditions for optimal efficiency in our experiment, and discuss echo attenuation rates.

REPEATED MAGIC SANDWICHES

We use a Repeated Magic Sandwich (RMS) sequence:

90 • -τ -MS -2τ -MS -2τ -...
with MS cycles of duration 4τ . This is analogous to a conventional spin echo sequence in which the 180 • pulses have been replaced by MS sequences. We show below that the RMS sequence can be designed to refocus the phase dispersions caused by both linear and dipolar couplings, and compare spin echo decay rates for RMS sequences and Carr-Purcell-Meiboom-Gill (CPMG) sequences [18] in the limit of negligible dipolar interactions (Sect. ). RMS echo trains at high magnetizations are then described, and experimental and numerical results are compared (Sect. ).

Repeated magic sandwiches with negligible dipolar interactions

Expectations

The magic sandwiches used for the RMS sequence are of the form:

90 • y , (360 • x , -360 • x ) n , 90 • y .
Note that the sign of the last rotation has been changed with respect to the basic MS sequence (3) in order to remove the linear phase dispersion introduced by static field inhomogeneities (residual inhomogeneities or applied gradients).

To experimentally study echo attenuation on time scales comparable to those relevant in dipolar-coupled systems, a static field gradient is applied to the weakly polarized liquid mixtures to induce controlled spatial modulation of the magnetization. During each MS, the intense rf field has a dominant contribution in Eq. 1, and the evolution due to the gradient is frozen. The RMS sequence may thus be considered as being equivalent to a CPMG sequence in which the gradient is only applied between the rf pulses, i.e. for one third of the actual period of the RMS sequence (the period of the RMS sequence is T RMS = 6τ , while that of the conventional CPMG sequence is T CPMG = 2τ ).

Since spatial variations of M are damped by diffusion, echo amplitudes for both sequences are expected to decay exponentially in time with a rate:

1 T 2 = Dk 2 (γGT seq ) 2 12 , ( 4
)
where G is the amplitude of the gradient of the z component of B ext , T seq is the period of the sequence, and k a dimensionless factor [19]. This factor is given by:

k 2 = δ 2 (3T seq -2δ)/T seq 3 , ( 5
)
where δ is the effective duration of the applied gradient G. Hence, k 2 = 1 for the CPMG sequence (δ = T CPMG ) and k 2 = 7/27 for the RMS sequence (δ = T RMS /3). A linear dependence of the decay rate on (GT seq ) 2 is thus expected for both pulse sequences.

Experiments at very low magnetizations

A weakly polarized sample is easily prepared by addition of a small quantity of laser-polarized 3 He gas to a preliquified unpolarized mixture. The added 3 He quantity (∼ 0.5 % of the amount of 3 He in the sample) is sufficiently small to have negligible influence on the 3 He concentration, hence on the diffusion coefficient. Series of experiments can thus be performed on nearly identical samples.

Figure 4 shows echo trains obtained with RMS and CPMG sequences in the presence of an applied gradient. For both sequences, echo amplitudes are observed to decay exponentially with time. The decay rates 1/T 2 are deduced from exponential fits on the squared echo amplitudes. This self-weighting of data reduces a systematic bias in the analysis associated with Rician noise in amplitude data [20]. In our preliminary investigations, systematic measurements with both sequences have not been performed on identical samples.

Figure 5 shows a compilation of RMS and CPMG decay rates measured in liquid mixtures with 3 He molar fraction X = 1.9% (sample 1) and X = 4.5% (sample 2). Data are plotted as a function of the reduced parameter )

k 2 (γ G T seq ) 2 /12 (cm -2 )
FIG. 5: Comparison of RMS and CPMG decay rates under conditions where dipolar interactions are negligible. For sample 1 (X=1.9%, T =1.15 K), RMS measurements are performed with TRMS=145 ms and various gradient amplitudes (solid triangles), and a single CPMG measurement is performed with TCPMG=151 ms at G=1.1 µT/cm (open square). For sample 2 (X=4.5%, T =1.13 K), CPMG measurements are also performed with TCPMG=151 ms (solid squares). These data are used to infer limiting values for CPMG decay rates in sample 1 (dashed and dotted lines, see text). Solid lines a and b are linear fits to the measured RMS and CPMG decay rates, respectively.

k 2 (γGT seq ) 2 /12. Both sets of data exhibit the expected linear behavior in (GT seq ) 2 . The CPMG data (solid squares) reveal a small, but finite decay rate at zero applied gradient (imperfect shimming). The RMS decay rates (solid triangles) exhibit a larger intercept at G=0, and a more rapid increase with (GT RMS ) 2 . Much higher decay rates are indeed expected because of the lower 3 He concentration, as confirmed by the CPMG measurement performed on sample 1 (open square).

A second CPMG measurement on sample 1 would be required to quantitatively compare the attenuation induced by both CPMG and RMS sequences. The CPMG decay rate at G=0 can be tentatively inferred from that of sample 2, since conditions were comparable in both sets of experiments. The dashed line in Fig. 5 is obtained assuming these two rates to be equal (diffusion-independent decay rate due, e.g., to sequence imperfections). The dotted line is obtained by appropriate scaling assuming, instead, that decay rates are solely due to diffusion-induced attenuation (in residual inhomogeneities). Actual CPMG decay rates are expected to lie between those two boundaries. In the following subsection, they are compared to RMS decay rates measured using sample 1.

Discussion

The slope of each line in Fig. 5 directly measures the diffusion coefficient for the corresponding helium mixture (cf. Eq. 4). As expected, the linear fit to RMS data and the lines demarcating CPMG decay rate limits in sample 1 have very similar slopes (they only differ by 1.5% for the dotted line and 8.5% for the dashed line). Direct comparison of RMS and CPMG decay rates at fixed 3 He molar fraction is in progress. The RMS decay rates in sample 1 yield D=2.7 × 10 -3 cm 2 /s for X=1.9%. The diffusion coefficient measured with the CPMG sequence in sample 2 is smaller (D=1.3 × 10 -3 cm 2 /s for X=4.5%). Both values agree with published data [21].

At zero gradient, the echo decay rate is higher for the RMS sequence than for the CPMG sequence. The physical origin of this additional damping has not yet been identified. Evaporation of 3 He atoms due to rf heating (of order 1 mK/s during RMS experiments) has a negligible influence on signal decay. Additional signal losses incurred during RMS sequences, that are observed to be B dip -independent, may be caused by minor rf amplitude and timing imperfections that will soon be addressed.

Repeated magic sandwiches in the presence of dipolar interactions

The upper panel of Fig. 6 displays two echo trains recorded for RMS experiments performed at zero gradient on highly polarized samples. It illustrates the spectacular change in time behavior of the signal amplitudes when the MS repetition time is modified. For the short T RMS , the echo amplitude slowly and steadily decreases. For the long T RMS , the echo amplitude rapidly decreases during the first seconds, then decays as slowly as for the short T RMS . Dynamic stabilization of transverse magnetization by rf-driven time reversal yields a tremendous increase in precession lifetime: 40 s is measured for the slow signal decays, which is three orders of magnitude larger than the 70 ms obtained for a FID at the same B dip (Fig. 1,trace b). This stabilization results from fast switching between free evolution and rf-driven evolution, i.e., from the frequent sign inversion of the evolutionary effect of dipolar fields. The initial rapid decrease observed at long T RMS is reminiscent of the strong DDF-induced attenuation observed at high magnetizations in the single-shot experiments with MS duration 2τ (Fig. 3). Finally, when the RMS pulse sequence is stopped (about 9 s after the initial 90 • pulse), signal amplitudes abruptly decrease as they do for simple FIDs following 90 • tipping angle pulses (cf. Fig. 1).

The lower panel of Fig. 6 displays RMS echo trains obtained by numerical simulation. We use a general-purpose model to compute NMR dynamics based on an exact calculation of the time evolution of coupled magnetic moments on a cubic lattice [8,16]. All processes that contribute to the evolution of M are taken into account, including diffusion, dipolar interactions, rf field, and coupling to the detection coil (radiation damping). However, this coarse-grained description of a continuous fluid has obvious limitations when spatial variations at the scale of the lattice constant are involved. Computations are performed in the rf rotating frame using a standard secular approximation. The magnetic field induced at each site by the remainder of the sample is efficiently computed by toggling between real and Fourier spaces [22]. The time evolution is computed by integrating the nonlinear Bloch equation using a standard Runge-Kutta technique. Free evolution usually involves slow local changes, and computations are quite fast on a standard PC computer. When intense rf fields are applied, the induced precession is fast and the computing load is significantly increased (e.g., by two orders of magnitude for a magic sandwich). Results presented in Fig. 6 have been obtained using 16 × 16 × 16 sites, with computation times of order 1000 s per second of physical evolution time. For these exploratory simulations, the diffusion coefficient D has been deliberately increased fivefold with respect to the experimental value to efficiently damp spatial variations of M at scales smaller than the lattice constant, as required to obtain meaningful results. To obtain traces c and f, input parameters B dip and T RMS have been set to their experimental values for traces a and b, respectively. T RMS has been set to intermediate values for traces d and e (see caption). The simulated echo trains exhibit the same qualitative features as the experimental data: an initial rapid decrease of echo amplitudes (for all but the shortest T RMS ), a plateau, and an abrupt signal decrease following the last magic sandwich.

The plateau regime corresponds to the dynamic stabilization of transverse magnetization obtained in the experiments, but the observed slow decay of echo amplitudes is not reproduced. A flat plateau is obtained for all simulations. The computed magnetization maps show that M is almost uniform, which explains why no noticeable attenuation occurs despite the use of a large diffusion coefficient. The experimental slow decays must thus result from other effects (cf. Sect. ). The amplitude of the plateau increases when T RMS is decreased (lower panel), as experimentally observed (upper panel and ref. [11]). It also increases when D is increased (data not shown). This indicates that a higher transverse magnetization can be dynamically stabilized when the development of unstable magnetization patterns is more efficiently prevented. A reliable computation of the evolution of M during the first seconds will require using a larger number of sites to more accurately describe the development of inhomogeneities at small length scales.

FIG. 1 :Figure 1

 11 FIG.1: FID response of laser-polarized3 He in liquid4 He to a 45 • tipping pulse for B dip =51 nT (negligible dipolar interactions, trace a) and a 90 • tipping pulse for B dip =0.9 µT (trace b). The inset shows the deviation of trace b from its initial value S0, in semi-log scale. Experiments are performed on the same sample at T = 1.14 K, with X ≃ 2%

FIG. 3 :

 3 FIG.3: Spin echo formation using a 90 • -τ -MS sequence with MS duration 2τ =140 ms. Data are shown for increasing values of B dip , from 0.3 µT to 1.5 µT (from back to front). FID signals have been normalized to facilitate comparison of data. Increasing B dip decreases the echo amplitude (inset), and significantly alters its shape. Experiments are performed at T ≃1.14 K for the 3 He molar fractions X indicated in the legend.

FIG. 4 :

 4 FIG.4: Echo trains obtained under conditions where dipolar interactions are negligible using a CPMG sequence (TCPMG = 151 ms, trace a) and a RMS sequence (TRMS = 145 ms, trace b). Both experiments are performed on the same sample (X=1.9%) at T =1.15 K with G=1.1 µT/cm. Decay rates 1/T2 of echo amplitudes are 0.28 s -1 (trace a) and 0.12 s -1 (trace b).

FIG. 6 :

 6 FIG.6: Experimental and simulated RMS echo trains at zero gradient. Upper panel: Experiments (traces a and b) are performed with TRMS=97 ms and 241 ms, respectively, for identical conditions: B dip =0.9 µT, X=4.5%, and T =1.12 K. Lower panel : Numerical simulations (traces c to f) are performed using the experimental value of B dip , but a larger diffusion coefficient: D=6.4 × 10 -3 cm 2 /s instead of D=1.3 × 10 -3 cm 2 /s (see text). The MS repetition time increases from trace c to trace f, with TRMS= 97, 180, 217, and 241 ms for these computations.


Christophe Chesneau
email: chesneau@math.unicaen.fr

A tail bound for sums of independent random variables : application to the symmetric Pareto distribution

Keywords: AMS 2000 subject classifications: 60E15 Tail bound, symmetric Pareto distribution imsart-generic ver Sum-Pareto.tex date: November 20, 2007

published or not. The documents may come

A tail bound for sums of independent random variables : application to the symmetric Pareto distribution

MOTIVATION

Let (Y i ) i∈N * be independent random variables. For any n ∈ N * , we wish to determine the smallest sequence of functions p n (t) such that

P n i=1 Y i ≥ t ≤ p n (t), t ∈ [0, ∞[.
This problem is well-known; numerous results exist. The most famous of them is the Markov inequality. Under mild assumptions on the moments of the X i 's, it gives a polynomial bound p n (t). In many cases, this bound can be improved. For instance, if the X i 's are almost surely absolutely bounded, or admit finite moments of all orders (and these moments satisfy some inequalities), the Bernstein inequalities provide better results. The obtained bounds p n (t) are exponential. See [START_REF] Petrov | Limit Theorems of Probability Theory[END_REF] and [START_REF] Pollard | Convergence of Stochastic Processes[END_REF] for further details and complete bibliography.

In this note, we present a new inequality which provides a bound p n (t) of the form p n (t) = v n (t) + w n (t), where v n (t) is polynomial, and w n (t) is exponential. It can be applied under mild assumptions on the X i 's; as for the Markov inequality, only knowledge of the order of a finite moment is required. The main interest of our inequality is that it can be applied when the 'Bernstein conditions' are not satisfied, and can give better results than the Markov inequality. In order to illustrate this, we investigate the bound of the tail probability for a sum of n weighted i.i.d. random variables having the symmetric Pareto distribution. This is particularly interesting because the exact expression of the distribution of such a sum is really difficult to identify. See, for instance, [START_REF] Ramsay | The distribution of sums of certain i.i.d. Pareto variates[END_REF]. Moreover, there are some applications in economics, actuarial science, survival analysis and queuing networks.

The note is organized as follows. Section 2 presents the main result. In Section 3 we illustrate the use of this result by considering the symmetric Pareto distribution. The technical proofs are postponed to Section 4.

MAIN RESULT

Theorem 2.1 below presents a bound of the tail probability for a sum of n independent random variables. As mentioned in Section 1, it requires knowledge only of the order of a finite moment.

Theorem 2.1. Let (Y i ) i∈N * be independent random variables. We suppose that

• for any n ∈ N * , and any i ∈ {1, ..., n}, we have, w.l.o.g., E(Y i ) = 0, • there exists a real number p ≥ 2 such that, for any n ∈ N * , and any i ∈ {1, ..., n}, we have

E(|Y i | p ) < ∞.
Then, for any t > 0, and any n ∈ N * , we have

P n i=1 Y i ≥ t ≤ C p t -p max r n,p (t), (r n,2 (t)) p/2 + exp - t 2 16b n , (2.1)
where, for any u ∈ {2, p}, r n,u (t) =

n i=1 E |Y i | u 1 {|Yi|≥ 3bn t } , b n = n i=1 E Y 2 i and C p = 2 2p+1 max p p , p p/2+1 e p ∞ 0 x p/2-1 (1 -x) -p dx .
The proof of Theorem 2.1 uses truncation technics, the Rosenthal inequality and one of the Bernstein inequalities. See [START_REF] Rosenthal | On the subspaces of L p (p ≥ 2) spanned by sequences of independent random variables[END_REF] and [START_REF] Petrov | Limit Theorems of Probability Theory[END_REF].

Clearly, Theorem 2.1 can be applied for a wide class of random variables. However, if the variables are almost surely absolutely bounded, or have finite moments of all orders, the Bernstein inequalities can give more optimal results than (2.1). But, when these conditions are not satisfied, Theorem 2.1 becomes of interest. This fact is illustrated in Section 3 below for the symmetric Pareto distribution. Other examples can be studied in a similar fashion. Proposition 3.1. Let s > 2 and (X i ) i∈N * be i.i.d. random variables with the probability density function f (x) = 2 -1 s|x| -s-1 1 {|x|≥1} . Let (a i ) i∈N * be a sequence of nonzeros real numbers such that n i=1 |a i | s < ∞. Then, for any n ∈ N * , any t ∈ (0, 3bn ρn ), where ρ n = (

APPLICATION: SYMMETRIC PARETO DISTRIBUTION

n i=1 |a i | s )
1/s , and any p ∈ (2, s),

we have P n i=1 a i X i ≥ t ≤ K p t -2p+s b p-s n n i=1 |a i | s + exp - t 2 16b n , (3.1) where b n = s s-2 n i=1 a 2 i , K p = 3 p-s max s s-p , s s-2 p/2
C p , and C p = 2 2p+1 max p p , p p/2+1 e p ∞ 0 x p/2-1 (1x) -p dx . Notice that, since the distribution of the variables is symmetric, the constant C p (associated to the Rosenthal inequality) can be improved. For its optimal form, we refer to [START_REF] Ibragimov | On an exact constant for the Rosenthal inequality[END_REF].

In the literature, there exist several results on the approximation of the tail probability of a sum of n i.i.d. random variables having the symmetric Pareto distribution. But, to our knowledge, contrary to Proposition 3.1, these results are asymptotic (i.e. t → ∞). See, for instance, [START_REF] Goovaerts | The Tail Probability of Discounted Sums of Pareto-like Losses in Insurance[END_REF].

Illustration. Here, we consider a simple example to compare the precision between (3.1) and the bound obtained via the Markov inequality.

Let s > 2 and (X i ) i∈N * be i.i.d. random variables with the probability density function f (x) = 2 -1 s|x| -s-1 1 {|x|≥1} . For any integer n such that n 1/2-1/s (log n) -1/2 > 2 3/2 3 s 1/2 , and any p ∈ max( s 2 , 2), s , if we take t = t n = 2 3/2 (sn log n) 1/2 , then we can balance the two terms of the bound in (3.1); there exist two constants, Q 1 > 0 and Q 2 > 0, such that

P n i=1 X i ≥ t n ≤ Q 1 n 1-s/2 (log n) -p+s/2 + n 1-s/2 ≤ Q 2 n 1-s/2 . (3.2)
Under the same framework, for any p < s, the Markov inequality combined with the Rosenthal inequality (see Lemma 4.1 below) implies the existence of two constants, R 1 > 0 and R 2 > 0, such that

P n i=1 X i ≥ t n ≤ t -p n E n i=1 X i p ≤ R 1 t -p n n p/2 ≤ R 2 (log n) -p/2 . (3.3)
Therefore, for n large enough, the rate of convergence in (3.2) is really faster than those in (3.3). In this case, (3.1) gives a better result than the Markov inequality.

imsart-generic ver. 2007/09/18 file: Sum-Pareto.tex date: November 20, 2007

PROOFS

. Proof of Theorem 2.1. Let n ∈ N * . For any t > 0, we have

P n i=1 Y i ≥ t = P n i=1 (Y i -E (Y i )) ≥ t ≤ U + V,
where

U = P n i=1 Y i 1 {|Yi|≥ 3bn t } -E Y i 1 {|Yi|≥ 3bn t } ≥ t 2
and

V = P n i=1 Y i 1 {|Yi|< 3bn t } -E Y i 1 {|Yi|< 3bn t } ≥ t 2 .
Let us bound U and V , in turn.

The upper bound for U . The Markov inequality yields

U ≤ 2 p t -p E n i=1 Y i 1 {|Yi|≥ 3bn t } -E Y i 1 {|Yi|≥ 3bn t } p . (4.1)
Now, let us introduce the Rosenthal inequality. See [START_REF] Rosenthal | On the subspaces of L p (p ≥ 2) spanned by sequences of independent random variables[END_REF].

Lemma 4.1 (Rosenthal's inequality). Let p ≥ 2 and (X i ) i∈N * be independent random variables such that, for any n ∈ N * , and any i ∈ {1, ..., n}, we have E(X i ) = 0 and E(|X i | p ) < ∞. Then we have

E n i=1 X i p ≤ c p max   n i=1 E (|X i | p ) , n i=1 E X 2 i p/2   ,
where c p = 2 max p p , p p/2+1 e p ∞ 0 x p/2-1 (1x) -p dx .

For any i ∈ {1, ..., n}, set

Z i = Y i 1 {|Yi|≥ 3bn t } -E Y i 1 {|Yi|≥ 3bn t } . Since E(Z i ) = 0 and E (|Z i | p ) ≤ 2 p E |Y i | p 1 {|Yi|≥ 3bn t } ≤ 2 p E (|Y i | p ) < ∞, Lemma 4.1 applied to the independent variables (Z i ) i∈N * gives E n i=1 Z i p ≤ c p max   n i=1 E (|Z i | p ) , n i=1 E Z 2 i p/2   , (4.2)
where c p = 2 max p p , p p/2+1 e p ∞ 0 x p/2-1 (1x) -p dx .

It follows from (4.1) and (4.2) that

U ≤ 2 p t -p c p max n i=1 E (|Z i | p ) , n i=1 E Z 2 i ≤ 2 2p t -p c p max   n i=1 E |Y i | p 1 {|Yi|≥ 3bn t } , n i=1 E Y 2 i 1 {|Yi|≥ 3bn t } p/2   = C p t -p max r n,p (t), (r n,2 (t)) p/2 , (4.3)
where

C p = 2 2p c p .
The upper bound for V . Let us present one of the Bernstein inequalities. See, for instance, [START_REF] Petrov | Limit Theorems of Probability Theory[END_REF].

Lemma 4.2 (Bernstein's inequality). Let (X i ) i∈N * be independent random variables such that, for any n ∈ N * and any i ∈ {1, ..., n}, we have E(X i ) = 0 and |X i | ≤ M < ∞. Then, for any λ > 0, and any n ∈ N * , we have

P n i=1 X i ≥ λ ≤ exp - λ 2 2(d 2 n + λM 3 )
,

where d 2 n = n i=1 E(X 2 i ).
For any i ∈ {1, ..., n}, set

Z i = Y i 1 {|Yi|< 3bn t } -E Y i 1 {|Yi|< 3bn t } . Since E(Z i ) = 0 and |Z i | ≤ |Y i |1 {|Yi|< 3bn t } + E |Y i |1 {|Yi|< 3bn t } ≤ 6bn
t , Lemma 4.2 applied with the independent variables (Z i ) i∈N * , and the parameters λ = t 2 and M = 6bn t , gives

V ≤ exp   - t 2 8 n i=1 V Y i 1 {|Yi|< 3bn t } + t 6 6bn t   . Since n i=1 V Y i 1 {|Yi|< 3bn t } ≤ n i=1 E Y 2 i = b n , it comes V ≤ exp - t 2 16b n . (4.4)
Putting (4.3) and (4.4) together, we obtain the inequality

P n i=1 Y i ≥ t ≤ U + V ≤ C p t -p max r n,p (t), (r n,2 (t)) p/2 + exp - t 2 16b n .
Theorem 2.1 is proved.

Y i ) = a i E(X i ) = 0 and n i=1 E(Y 2 i ) = s s-2 n i=1 a 2 i < ∞.
In order to apply Theorem 2.1, let us bound the term r n,u

(t) = n i=1 E |Y i | u 1 {|Yi|≥ 3bn t } = n i=1 |a i | u E |X i | u 1 |Xi|≥ 3bn |a i |t
for any u ∈ {2, p}, and any p ∈ max( s 2 , 2), s .

Recall that ρ n = ( C p , and C p = 2 2p+1 max p p , p p/2+1 e p ∞ 0 x p/2-1 (1x) -p dx . Proposition 3.1 is proved.

  Proposition 3.1 below investigates the bound of the tail probability for a sum of n weighted i.i.d. random variables having the symmetric Pareto distribution. imsart-generic ver. 2007/09/18 file: Sum-Pareto.tex date: November 20, 2007

  Proof of Proposition 3.1. Let n ∈ N * . Set, for any i ∈ {1, ..., n}, Y i = a i X i . Clearly, (Y i ) i∈N are independent random variables such that E(

imsart-generic ver. 2007/09/18 file: Sum-Pareto.tex date: November 20, 2007

.

  = sup i=1,...,n |a i |, we have E |X i | u 1 |Xi|≥ 3bn

					n i=1 |a i | s )	1/s . Since t ∈ 0, 3bn ρn	⊆ 0, 3bn σn , where
	σ n |a i |t s s-u 3bn |ai|t u-s . Hence,	= s	∞ 3bn |a i |t	x u-s-1 dx =
					r n,u (t) =	s s -u		3b n t	u-s n i=1	|a i | s .
	Therefore,
	max r n,p (t), (r n,2 (t)) p/2 ≤ R p	3b n t	p	max	  3b n tρ n	-s	,	3b n tρ n	-s p/2	  ,
	where R p = max	s s-p , s s-2	p/2	.
	Since t ∈ 0, 3bn ρn	and p > 2, we have max	3bn tρn	-s	,	3bn tρn	-s p/2	=
	3bn tρn	-s	. Hence,
			max r n,p (t), (r n,2 (t)) p/2 ≤ R p		3b n t	p-s n i=1	|a i | s .	(4.5)
	Putting (4.5) in Theorem 2.1, we obtain
		P	n i=1	a i X i ≥ t ≤ K p t -2p+s b p-s n	n i=1	|a i | s + exp -	t 2 16b n	,
	where b n =	s s-2	n i=1 a 2 i , K p = 3 p-s max	s s-p , s s-2	p/2

imsart-generic ver. 2007/09/18 file: Sum-Pareto.tex date: November 20, 2007


Jun Shen

Nabil Zemiti

Jean-Louis Dillenseger

Philippe Poignet

Fast and Simple Automatic 3D Ultrasound Probe Calibration Based on 3D Printed Phantom and an Untracked Marker

Tracking the pose of an ultrasound (US) probe is essential for an intraoperative US-based navigation system. The tracking requires mounting a marker on the US probe and calibrating the probe. The goal of the US probe calibration is to determine the rigid transformation between the coordinate system (CS) of the image and the CS of the marker mounted on the probe. We present a fast and automatic calibration method based on a 3D printed phantom and an untracked marker for three-dimensional (3D) US probe calibration. To simplify the conventional calibration procedures using and tracking at least two markers, we used only one marker and did not track it in the whole calibration process. Our automatic calibration method is fast, simple and does not require any experience from the user. The performance of our calibration method was evaluated by point reconstruction tests. The root mean square (RMS) of the point reconstruction errors was 1.39 mm.

I. INTRODUCTION

Ultrasound (US) imaging is a safe, non-ionizing, low cost and portable medical imaging modality. Three-dimensional (3D) imaging can generate a 3D US volume of the region of interest (ROI) which provides more information about the anatomical structures than 2D US imaging. There are two kinds of 3D US imaging techniques commonly being used by clinicians. The first one is called freehand 3D US system [START_REF] Gee | Engineering a freehand 3D ultrasound system[END_REF] and uses a marker mounted on a 2D US probe. When the probe sweeps over the ROI, the acquired slices are stacked to reconstruct a 3D volume according to their corresponding locations in the space. However, the disadvantages of using freehand 3D US system are that it is time-consuming and does not provide real-time US volume acquisition, which is required for many intraoperative applications. In order to solve these problems, the second imaging technique generates a 3D US volume by using a dedicated 3D probe. For instance, the 2D matrix array transducer is able to perform volumetric imaging in realtime and the swept motor transducer acquires a 3D volume in 1 to 4 seconds.

Tracking the pose of the ultrasound probe is required for many intraoperative applications. However, tracking systems **This work was supported in part by the French ANR within the Investissements d'Avenir Program (Labex CAMI, ANR-11-LABX0004, Labex NUMEV, ANR-10-LABX-20, and the Equipex ROBOTEX, ANR-10-EQPX-44-01), the ARC Foundation and the Région Bretagne, France.

1 Jun Shen, Nabil Zemiti and Philippe Poignet are with the Laboratory of LIRMM, University of Montpellier, CNRS, Montpellier, France. Email: Jun.Shen@lirmm.fr, Nabil.Zemiti@lirmm.fr, poignet@lirmm.fr 2 Jun Shen and Jean-Louis Dillenseger are with INSERM UMR 1099, Rennes, France; Universiy of Rennes 1, LTSI, Rennes, France. Email:

jun.shen@univ-rennes1.fr, jean-louis.dillenseger@univ-rennes1.fr generally only record the pose of a marker (rigidly mounted on the US probe) with respect to a fixed coordinate system (CS). Therefore, it is essential to calibrate the pose of the US image with respect to this marker. The process of determining the rigid transformation between the CS of the US image and the CS of the marker (mounted on the probe) is called US probe calibration [START_REF] Mercier | A review of calibration techniques for freehand 3-D ultrasound systems[END_REF]. Common calibration approaches are based on calibration phantoms with known geometry. These phantoms can be points [3][4], wires [START_REF] Carr | Surface reconstruction in 3D medical imaging[END_REF] [START_REF] Prager | Rapid calibration for 3-D freehand ultrasound[END_REF] or planes [6][7]. The calibration problem of freehand 3D US system has been studied extensively, referring to the reviews of many published techniques [2][8]. However, the calibration procedure for calibrating freehand 3D US systems is different from calibrating 3D US probes. This is due to the fact that freehand 3D US system calibration uses points or 2D features extracted from US slices, while 3D US probe calibration uses volumetric features extracted from a 3D volume. In this study, we focus on 3D US probe calibration techniques.

Bergmeir et al. presented in [START_REF] Bergmeir | Comparing calibration approaches for 3D ultrasound probes[END_REF] a calibration method based on a tracked phantom. This method achieved a higher accuracy than another method performing US imaging around an untracked fixed phantom and registering the US volumes to each other, as in hand-eye calibration. Some works directly used the tip of a tracked stylus as an imaging feature for the probe calibration [START_REF] Poon | Comparison of calibration methods for spatial tracking of a 3-D ultrasound probe[END_REF]. However, the tracked stylus had to be calibrated before using it, and the stylus calibration introduced non-negligible errors, as reported in [START_REF] Hsu | Accurate fiducial location for freehand 3D ultrasound calibration[END_REF]. Some other works used planar phantoms and the phantoms' positions had to be determined in a fixed CS by a tracked calibrated stylus [START_REF] Lange | Automatic calibration of 3D ultrasound probes[END_REF]- [START_REF] Abeysekera | Calibration for position tracking of swept motor 3-d ultrasound[END_REF]. Their methods suffer not only from the stylus calibration errors, but also from phantom calibration errors. The use of tracked wire phantoms, as proposed in [START_REF] Bergmeir | Comparing calibration approaches for 3D ultrasound probes[END_REF] and [START_REF] Poon | Comparison of calibration methods for spatial tracking of a 3-D ultrasound probe[END_REF], includes challenges in manufacturing the wires precisely, ensuring they are straight, and correctly positioning them in water. Unlike the above calibration techniques using 3D US volumes, a study presented in [START_REF] Abeysekera | Calibration for position tracking of swept motor 3-d ultrasound[END_REF] performed calibration on multiple slices before the 3D scan conversion process. Then, they found the best-fit path through the multiple slices calibrations and used it as the final calibration solution for a swept motor 3D probe. Although, they achieved a high calibration accuracy, their calibration procedure seems timeconsuming for multiple US slices. In addition, it seems difficult to implement this method for calibrating 2D matrix array transducers which perform volumetric imaging by electronic scanning.

To the best of our knowledge, all the 3D US probe calibration methods proposed in the literature need tracking systems to track at least two markers: the first one is rigidly mounted on the US probe, and the second one is fixed to the calibration phantom. Indeed, some of the methods need to track a third marker (attached on a stylus) for the phantom calibration to find the relationship between the features and the marker of the phantom. Since using and tracking more markers may introduce more uncertainties, we propose to use only one marker and not to track it during the calibration procedure, so that decreases the source of measurement uncertainties. This also avoids errors due to stylus calibration and phantom calibration steps.

In this study, a fast and automatic calibration method based on a 3D printed phantom and an untracked marker is proposed for 3D US probe calibration. The calibration solution is easily found by a rigid registration between the phantom's computer-aided design (CAD) model and the phantom's US image. This is because the marker's CAD model was used during the design of the phantom's CAD model, so that the position of the phantom is known in the marker's CS. The proposed automatic method is easy to implement and does not require any experience from the user. In addition, our method not only works for calibrating swept motor 3D probes, but also for calibrating 2D matrix array transducers.

II. METHOD

Fig. 1 shows all the CSs and transformations used in a conventional US calibration procedure. As we can see, a tracking system has two mobile parts (marker n1 and n2) and a stationary counterpart (an optical pose-tracking system). Usually, the stationary part of the tracking system is defined as world coordinate system (WCS) w, and the markers' CSs and CAD models are provided by the manufacturer. The first marker n1 is rigidly fixed on the probe and its CS is represented by m. The second marker n2 is mounted on the phantom and the phantom calibration is required to relate the CS of the phantom to the CS of the marker n2. As classically presented in the literature, the CS p refers to the CS of the marker n2 after phantom calibration.

The goal of the calibration is to determine the rigid transformation m T i between the CS of the 3D US image i and the CS of the marker n1 m. In Fig. 1, the conventional mathematical calibration framework can be presented as

m T i = w T -1 m w T p p T i (1)
where b T a is the transformation from the CS a to CS b, and i, m, p and w represent the CS of the 3D US image, the marker n1 mounted on the probe, the marker n2 attached to the phantom, the world, respectively. Therefore, to obtain the calibration solution m T i , the main issues below have to be solved: 1) Phantom calibration: using a calibrated stylus to determine the phantom position in the CS of the marker n2 and thus in WCS, so that w T p is measured. 2) Registration of the phantom's US image and the phantom's geometries, so that p T i is obtained.

3) Tracking of the markers n1 and n2 computing the transformations w T -1 m and w T p . In order to simplify this calibration procedure, we propose a method based on a custom-made phantom, in which solving only problem 2) is enough to find the calibration solution. As shown in Fig. 2, a marker is assembled and secured on the probe. There is no tracking of the marker needed in the whole calibration process, so we call this marker 'untracked marker'. The phantom is mounted on the untracked marker and its position is known in this marker's CS. This is due to that the marker's CAD model was used during the design of the phantom, and their CSs m and p was coincided with each other. Thus, the conventional mathematical calibration framework of Eq.1 can be simplified to

m T i = p T i (2)
Therefore, solving a series of problems in the conventional calibration procedure becomes solving only the registration problem between the mesh model generated from the US image and the phantom's CAD model. There is no tracking of the marker needed in our method. Coordinate systems and transformations used in conventional calibration method: tracking system tracks the marker n1 and n2.

A. Calibration Phantom

The calibration phantom was created as shown in Fig. 3(a). At one end (right part of the Fig. 3 3(a)) was designed to be attached to the marker. As shown in Fig. 4, the hole in the middle of the phantom perfectly fits the protrusion of the marker, and the protrusion of p3 can be perfectly embedded into p3 . p1 and p1 , p2 and p2 having the same diameter were fixed together by two screws. After attaching the phantom to the marker, the CS p in Fig. 4 is the same as the CS m.

A 3D modeling software SketchUp (Trimble Inc.) was used to create the CAD model of our calibration phantom, and the designed model was 3D printed by the Stratasys Fortus 400mc prototyping machine with Polycarbonates material. The printing resolution was 0.178 mm. The CAD model of the phantom can be shared upon request. The marker used in this study was the Boomerang active marker from Atracsys with RMS calibration error of 0.04 mm. The phantom was designed to adapt to this marker, but it can adapt to other markers after few modifications according to other markers' CAD model.

B. Data Acquisition

Before performing 3D US imaging on the calibration phantom, the phantom, the untracked marker and the probe must be assembled together, as shown in Fig. 2. The marker was rigidly mounted on the probe by a 3D printed mounting frame, then, the phantom was mounted on this marker, as explained in Fig. 4. The US data acquisition was simply done by placing the calibration setup into water (Fig. 5(a)).

In this study, the 3D US image was acquired by a Sonix-Touch Q+ ultrasound machine (BK Ultrasound, powered by Analogic Corporation), equipped with the 4DL14-5/38 Linear 4D US probe. The ultrasound machine performed an US acquisition on the phantom at a transmit frequency of 10 MHz and a depth of 75 mm. These parameters were chosen to be compatible with our clinical application for tongue base surgery, but they can be adjusted for other applications.

C. US to Phantom's CAD model Registration

The goal of the calibration process is to determine the rigid transformation between the CS of the 3D US image and the CS of the marker which is mounted on the probe. As shown in Eq.2, our calibration solution is found by solving only the rigid registration problem between two mesh models from the US image and the phantom. As presented in Fig. 6, the olive green model was the CAD model of the phantom, and the orange mesh model was generated from an automatic segmentation of the US volume. The registration between these two 3D mesh models was accomplished by using 3D Slicer [START_REF] Fedorov | 3D Slicer as an image computing platform for the quantitative imaging network[END_REF] surface registration. It took few seconds to compute this transformation.

The automatic segmentation was achieved by computing the directional gradients of the 3D US image and then applying Standard Hough Transform (SHT) [START_REF] Gerig | Fast contour identification through efficient Hough transform and simplified interpretation strategy[END_REF] on the gradients image. Thus, 5 horizontal lines and 5 vertical lines on the plane having 4 circles were detected. The US image was automatically segmented based on the position information of these lines segments in the US image and the known geometry of the phantom.

After completing the calibration process, the phantom was removed by unscrewing it from the marker. Meanwhile, there was no relative movement between the marker and the probe, because the marker was mounted and secured on the probe by a mounting frame.

III. RESULTS

A. Feature Exaction Accuracy

In order to evaluate the registration performed in subsection II-C, we used the target registration error (TRE) [START_REF] Fitzpatrick | Predicting error in rigid-body point-based registration[END_REF] which computes the distance, after registration, between corresponding points. First, landmarks were manually placed on the inner corners of the box and the centers of the circles in the 3D US image. Then, these landmarks and the corresponding features in the phantom's CAD model were used for the TRE computation. The root mean square of the target registration error (RmsTRE) was computed to evaluate the registration performance.

We used 10 US volumes from different acquisitions and automatically segmented them. Then, the registration was performed on each of the mesh models from the automatic segmentation and the phantom's CAD model generating 10 transformations ( m T i ) 1...10 . These transformations were used to compute the RmsTRE 1...10 . The mean of the RmsTREs was 0.361 mm and the standard deviation was 0.156 mm. The calibration procedure (including assembling calibration setup, US acquisition, automatic segmentation and mesh model registration) took 3 to 5 minutes. The user does not need to be trained.

We found that our automatic method and the manual method gave similar RmsTRE 1...10 values. A non-specialist was asked to manually segment the US volumes over several days. Then the registration was performed between each of the mesh models from the manual segmentation and the phantom's CAD model. Then, 10 transformations ( m T i ) 1...10 were obtained and used to compute the RmsTRE 1...10 . The mean of the RmsTREs was 0.344 mm and the standard deviation was 0.079 mm, which was slightly better than the values obtained by our automatic method. However, manual segmentation was time consuming and user-dependent. As shown in Fig. 7, the RmsTRE values from automatic method have bigger variance than the values from manual method. This is because the automatic method is more sensitive to the image quality than the manual method. The US image quality can be improved by adjusting imaging parameters during the data acquisition. RmsTRE values obtained by manual method and automatic method during registration evaluation. The smaller RmsTRE means better registration performance.

B. Point Reconstruction Accuracy

Calibration accuracy refers to how close the transformation m T i obtained from the calibration method is to the real transformation between the probe and the US volume. We used point reconstruction tests, as presented in [START_REF] Abeysekera | Calibration for position tracking of swept motor 3-d ultrasound[END_REF], to evaluate our calibration method. The tests calculate the distance between a point (stylus tip) with a known location in the WCS w and this point transformed from US volume to the WCS w using our calibration solution. A tracked stylus was then used for the point reconstruction tests. The stylus is custom-made and consists of a marker at one end and a sharp tip at the other end. The point position of the sharp tip was obtained by a stylus calibration [START_REF] Leotta | Performance of a miniature magnetic position sensor for three-dimensional ultrasound imaging[END_REF] with an error of 0.83 mm.

We placed the stylus tip into water and the calibrated 3D probe was used to image the tip. 20 US volumes were acquired. Both the stylus and the probe were moved between each acquisition, meanwhile, their positions were recorded by the tracking system (WCS). The tip in US volume was manually located. The point reconstruction error was calculated by

e = w s -w T m m T i i s (3)
where w s represents the coordinate of stylus tip in the WCS, i s represents the coordinate of the stylus tip in the US image, m T i represents our calibration solution and w T m is the transformation between the marker CS m and the stationary part of the tracking system (WCS).

In the point reconstruction tests, we used the Atracsys easyTrack 500 tracking system (0.2 mm RMS error at 1 m distance) with the Boomerang active marker (0.04 mm RMS calibration error). The stylus was 3D printed and its calibration error was 0.83 mm. The root mean square (RMS) of point reconstruction errors was 1.39 mm, mean of the point reconstruction errors was 1.26 mm and standard deviation was 0.62 mm.

IV. DISCUSSION

As a part of an intraoperative application, a 3D US probe must be calibrated and localized in a fast and simple way with high accuracy. Our calibration procedure takes less than 5 minutes, which is much faster than the method (taking 20 minutes) presented in the study of [START_REF] Baumann | 3-D ultrasound probe calibration for computer-guided diagnosis and therapy[END_REF].

The stylus used in the point reconstruction tests provides points with a known localization for validating our calibration method. However, the stylus is prone to errors and its pointer calibration normally has RMS errors in the range of [0.6 mm, 0.9 mm], and as hight as 1.5 mm [START_REF] Hsu | Accurate fiducial location for freehand 3D ultrasound calibration[END_REF]. In our study, the stylus calibration had an RMS error of 0.83 mm. This might introduce some errors into the point reconstruction accuracy, which may be improved with a better calibrated stylus.

The comparison of the calibration performance between our method and previous studies is difficult, due to the differences of experimental setup, protocol and environment. Furthermore, different studies use different definitions of point reconstruction accuracy to evaluate their methods. Some well known studies reported their results using RMS or mean of point reconstruction errors. For instance, Bergmeir et al. used a cross-wire phantom and compared hand-eye calibration (HE) and typical tracked phantom method (TP). They reported mean errors of 3.5 mm (HE) and 3.3 mm (TP) [START_REF] Bergmeir | Comparing calibration approaches for 3D ultrasound probes[END_REF]. Poon and Rohling compared methods based on an IXI-shaped wires phantom, a cube phantom and a stylus, and found RMS point reconstruction errors of 2.15 mm, 4.91 mm and 2.36 mm, respectively [START_REF] Poon | Comparison of calibration methods for spatial tracking of a 3-D ultrasound probe[END_REF]. Compared to above methods, we achieved smaller RMS and mean of point reconstruction errors of 1.39 mm and 1.26 mm, respectively. In [START_REF] Abeysekera | Calibration for position tracking of swept motor 3-d ultrasound[END_REF], they proposed performing calibration on multiple 2D slices, instead of a 3D volume. Then, they found the best-fit path through the multiple slices calibrations and used it as the calibration solution for a swept motor 3D probe. They achieved an RMS point reconstruction errors of 0.93 mm. The stylus used in their point reconstruction tests had a high accuracy which was computed by averaging the standard deviations of stylus location at 30 points achieving a value of 0.11 mm. However, our stylus was made in our laboratory and calibrated by sphere fit method with an accuracy of 0.83 mm. Our RMS point reconstruction accuracy will likely be improved by using a more accurate stylus in point reconstruction tests. Baumann et al. reported their method achieving an RMS point reconstruction error of 0.9 mm, however, instead of using a stylus, they used a bead phantom made of crossed wires in their point reconstruction tests [START_REF] Baumann | 3-D ultrasound probe calibration for computer-guided diagnosis and therapy[END_REF]. Therefore, as explained, it is very difficult to compare the performance among different methods, since they used different evaluation methods.

Besides the calibration precision discussed above, conventional calibration approaches use and track two markers: one is mounted on the probe and the other is attached to the calibration phantom. In addition, some techniques have introduced a third tracked marker (mounted on a stylus) for the phantom calibration [10][12][13] [START_REF] Abeysekera | Calibration for position tracking of swept motor 3-d ultrasound[END_REF]. Since using and tracking more markers may introduce more uncertainties into the calibration precision, we propose not to track any markers during the calibration procedure. Indeed, in conventional calibration approaches, multiple steps are prepared for the probe calibration (phantom calibration, stylus calibration, ...), which increase the calibration process duration. In our study, only one untracked marker was needed and mounted on the probe. The idea of assembling the phantom on this untracked marker, so that they share the same CS, makes the calibration procedure simple, fast and easy to implement. The fully automatic calibration method can be performed by a non-specialist. The whole calibration procedure (including assembling calibration setup, US acquisition, automatic segmentation and mesh model registration) takes less than 5 minutes.

V. CONCLUSION

A simple and efficient approach is presented for calibrating 3D US probes including 2D matrix array transducers or swept motor 3D probes. The idea of assembling the phan-tom, the untracked marker and the probe together greatly simplifies the conventional calibration procedures. The fully automatic calibration method does not require any experience or skills from the user. Future work will improve the stylus used in the validation process.

Fig. 1 .

 1 Fig. 1.Coordinate systems and transformations used in conventional calibration method: tracking system tracks the marker n1 and n2.

Fig. 2 .

 2 Fig. 2. Coordinate systems used in our calibration method: the calibration phantom CS p is the same as the untracked marker CS m.

  (a)), it has geometric features for US imaging. Fig.3(b) shows the axial view of this part which consists of a box shape (25 × 20 × 10 mm 3 ) and 4 holes with dimensions of ∅ = 8 mm, ∅ = 6 mm, ∅ = 4 mm and ∅ = 4 mm, respectively. These features were experimentally chosen in consideration of the US probe's field of view and the quality of the acquired image. The other end of the phantom (left part of the Fig.

Fig. 3 .

 3 Fig. 3. (a) CAD model of the calibration phantom; (b) Axial view of the phantom features used for US imaging.

Fig. 4 .

 4 Fig. 4. Phantom's CAD model designed to be attached on the marker. CS p = m in the design.

Fig. 5 (Fig. 5 .

 55 Fig. 5. (a) Performing US acquisition on the calibration phantom; (b) Sagittal (left) and axial (right) view of the acquired 3D US image.

Fig. 6 .

 6 Fig. 6. CAD mesh model of the phantom (olive green) and mesh model from US segmentation (orange) before (left) and after (right) registration.

Fig. 7 .

 7 Fig. 7.RmsTRE values obtained by manual method and automatic method during registration evaluation. The smaller RmsTRE means better registration performance.


est destinée au dépôt

Systematic study of the dynamics about and between the Libration Points of the Sun-Earth-Moon System

Bastien Le Bihan, Josep Masdemont, Gérard Gómez, Stéphanie Lizy-Destrez

To cite this version:

Bastien Le Bihan, Josep Masdemont, Gérard Gómez, Stéphanie Lizy-Destrez. Systematic study of the dynamics about and between the Libration Points of the Sun-Earth-Moon System. International Symposium on Space Flight Dynamics (ISSFD), Jun 2017, Matsuyama, Japan. pp.1-10. hal-01925297