Habib Ammari
email: ammari@cmapx.polytechnique.fr

Yves Capdeboscq
email: yves.capdeboscq@math.uvsq.fr

Hyeonbae Kang
email: hkang@math.snu.ac.kr

Eunjoo Kim
email: eunjoo@cmapx.polytechnique.fr

Mikyoung Lim
email: mklim@cmapx.polytechnique.fr

ATTAINABILITY BY SIMPLY CONNECTED DOMAINS OF OPTIMAL BOUNDS FOR THE POLARIZATION TENSOR

The notion of polarization tensor is employed for the derivation of the leading-order boundary perturbations in the steady-state voltage potentials that are due to the presence of conductivity inclusions of small diameter. Recently, Capdeboscq and Vogelius obtained optimal bounds of Hashin-Shtrikman type for the trace of the polarization tensor, showing that every pair satisfying these optimal bounds arises as the eigenvalues of a polarization tensor associated with a coated ellipse. In this paper, we give numerical evidence of the fact that the set of possible polarization tensor eigenvalue pairs can also be obtained using simply connected domains. Our numerical computations are based on a boundary integral method.

Introduction

This paper is concerned with the notion of polarization tensor (PT) associated with a bounded Lipschitz domain and an isotropic constant conductivity. The notion of PT appeared in problems of potential theory related to certain questions arising in hydrodynamics, in electrostatics, and in low-frequency scattering; see [START_REF] Shiffer | Virtual mass and polarization[END_REF][START_REF] Pólya | Isoperimetric Inequalities in Mathematical Physics[END_REF][START_REF] Dassios | Low Frequency Scattering[END_REF]. The PT is a key mathematical concept in efficiently imaging small conductivity inclusions from boundary measurements and also in calculating effective electrical properties of composite materials consisting of inclusions of one material of known shape embedded homogeneously into a continuous matrix of another having electrical properties different from its own. It is now known that the leading-order term in the boundary perturbations due to the presence of an inclusion inside a conductor as well as in the asymptotic expansion of the effective conductivity of a dilute composite material in terms of the volume fraction of the inclusions can be expressed by means of the PT of the inclusions shape; see [START_REF] Cedio-Fengya | Identification of conductivity imperfections of small diameter by boundary measurements: Continuous dependence and computational reconstruction[END_REF][START_REF] Brühl | A direct impedance tomography algorithm for locating small inhomogeneities[END_REF][START_REF]Optimal asymptotic estimates for the volume of internal inhomogeneities in terms of multiple boundary measurements[END_REF][START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF][START_REF] Ammari | Anisotropic polarization tensors and effective properties of anisotropic composite materials[END_REF][START_REF] Ammari | Boundary layer techniques for deriving the effective properties of composite materials[END_REF][START_REF] Jikov | Homogenization of Differential Operators and Integral Functionals[END_REF][START_REF] Milton | The Theory of Composites[END_REF][START_REF] Kleinman | Rayleigh scattering in Low and High Frequency Asymptotics[END_REF]. The asymptotic expansion of the effective conductivity is motivated by the practically important inverse problem of determining the volume fraction of a suspension of complicated shaped particles from boundary measurements of voltage potentials. Therefore, in both applications, it is important from an imaging point of view to precisely characterize the class of PTs and extract from the (detected) PT some important information, such as the size and the orientation of the inclusion. Indeed, it is important to obtain bounds for the matrix-elements of the PT since a better bound for the PT yields a better estimate for the size of the inclusion; see [START_REF]Optimal asymptotic estimates for the volume of internal inhomogeneities in terms of multiple boundary measurements[END_REF][START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF].

Recently, based on variational techniques originally described in [START_REF] Kohn | On bounding the effective conductivity of anisotropic composites, in Homogenization and Effective Moduli of Materials and Media[END_REF], Capdeboscq and Vogelius [START_REF]Optimal asymptotic estimates for the volume of internal inhomogeneities in terms of multiple boundary measurements[END_REF] obtained geometry independent optimal bounds of Hashin-Shtrikman (HS) type for the trace of the PT. These bounds had already been derived, by an alternative approach, by Lipton [START_REF] Lipton | Inequalities for electric and elastic polarization tensors with applications to random composites[END_REF]. Throughout this paper, these bounds are called the HS bounds. To mathematically state the main estimates provided in [START_REF]Optimal asymptotic estimates for the volume of internal inhomogeneities in terms of multiple boundary measurements[END_REF], let M denote the PT associated with the bounded domain D ⊂ R d whose volume |D| = 1 and the constant conductivity 0 < k = 1 < +∞. Then, for d = 2, 3, these bounds are as follows:

(1. [START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF] Trace(M )

≤ (k -1)(d -1 + 1 k ),
and

(1.2) Trace(M -1 ) ≤ d -1 + k k -1 .
In particular, in the two-dimensional case, if λ 1 and λ 2 are two eigenvalues of M , then

(1.3) λ 1 + λ 2 ≤ (k -1)(k + 1) k ,
and

(1.4) 1 λ 1 + 1 λ 2 ≤ k + 1 k -1 .
Figure 1 shows these bounds graphically, where the region A is the upper part of the HS bounds, namely the set of all (λ 1 , λ 2 ) satisfying

λ 1 + λ 2 < (k -1)(k + 1) k , 1 λ 1 + 1 λ 2 ≤ k + 1 k -1 , λ 1 ≤ λ 2 ,
and the region B is the set of all (λ 1 , λ 2 ) satisfying

λ 1 + λ 2 < (k -1)(k + 1) k , 1 λ 1 + 1 λ 2 ≤ k + 1 k -1 , λ 1 ≥ λ 2 .
The question we deal with in this paper is whether each point inside the HS bounds is a pair of eigenvalues of a PT associated with a domain of unit area in R 2 . This question was positively answered by Capdeboscq and Vogelius in [START_REF]A review of some recent work on impedance imaging for inhomogneieties of low volume fraction[END_REF]. These authors showed that each point inside the HS bounds is attained as a PT associated with a coated ellipse, or a washer of elliptic shape. In fact, every point on the lower bound 1/λ 1 + 1/λ 2 = (k + 1)/(k -1) corresponds to an ellipse, and as ellipses get thinner, corresponding points on the lower bound move to the upper or lower corner. If we start from an ellipse corresponding to a point on the lower bound, and make confocal washers of elliptic shape, then corresponding points move toward the upper bound following a certain curve as the washers get thinner and larger. These curves make foliations and cover all regions inside the bounds except the upper bound.

Our goal in this paper is to give numerical evidence of the fact that this optimal set can also be obtained using simply connected domains. Our approach is to compute numerically the polarization tensors corresponding to a radial interpolation

(k-1)/k k-1 k-1 (k-1)/k 1/x+1/y=(k+1)/(k-1) x+y=(k-1)(k+1)/k A B 0 Figure 1.
The optimal bounds for the polarization tensor in R 2 . between ellipses of unit area and crosses with right angles, equal legs, and smooth corners. It turns out that if we start from the unit disk and vary the domain to make a thin and long cross as shown in Figure 2, the corresponding eigenvalues move from the intersection point of the lower hyperbola and the line λ 1 = λ 2 toward the intersection point of the upper bound and the line λ 1 = λ 2 following the line λ 1 = λ 2 . Note that the intersection point of the lower hyperbola and the line λ 1 = λ 2 is the pair of eigenvalues of the PT associated with the unit disk. We also note that the cross-shaped domain in Figure 2 is invariant under rotation by π/2, and hence the corresponding PT is of the form λI for some λ where I is the 2 × 2 identity matrix. Thus, by interpolating a cross-shape domain and an ellipse, we can obtain foliation of the bounds (1.3) and (1.4). This is the basic idea of our work. Based on this idea, we will show, by numerical computations, that each point inside the HS bounds is attained by cross-shaped domains whose sides are different (see section 3 for precise parameterizations of domains). We will also present a method to find a domain whose PT has a given pair of eigenvalues. We emphasize that the result of this paper is numerical and analytic computations of the PTs associated with the domains we construct in this work, for arbitrary aspect ratios, seem unlikely. However, we will show that it is possible to compute the eigenvalues of mixtures of ellipses and cross-shape domains in the high aspect ratio limit.

. . . . . The paper is organized as follows. In section 2, we review the definition of the PT and its representation in terms of boundary integrals. This boundary integral representation is used for our numerical computations throughout this paper. In sections 3 and 4, we parameterize the domains we consider using two parameters, one being the ellipticity of the ellipses (or the length over width ratio of the crosses' legs) and the second being the proportion of ellipses and crosses, to show numerical evidence of the fact that the corresponding PTs make foliation of the HS bounds. In section 5 we first provide the asymptotic behavior of the PT associated with crosses as crosses get thinner and longer. We then show that this results validates the numerical results presented in section 4 for thin an long inclusions -which are the most difficult to accurately simulate numerically.

In section 6, we present a method to identify a parameterized domain whose PT has a given pair of eigenvalues. We derive an approximate formula to extrapolate from the computed eigenvalue pairs data the shape of the domain. We use a combination of least square fitting with Lagrange interpolation. The robustness of this extrapolation, compared to the direct simulation, is accurately tested.

In concluding this introduction, we briefly mention a conjecture of Pólya-Szegö. Observe that among the points in HS-bounds, the intersection point of the lower bound and the line λ 1 = λ 2 , (2(k -1)/(k + 1), 2(k -1)/(k + 1), has the minimal trace, and this is the eigenvalue of the PT associated with the disk of unit area. Pólya and Szegö [START_REF] Pólya | Isoperimetric Inequalities in Mathematical Physics[END_REF] conjectured that the disk (or the sphere for three dimensions) is a unique domain with minimal trace. This conjecture has not yet been proved.

Preliminaries

Let us recall the definition of PT. Let D be a bounded Lipschitz domain in R 2 . The conductivity of R 2 \ D and D are assumed to be 1 and 0 < k = 1 < +∞, respectively. Let W 1,2 be the collection of all functions ψ such that R 2 |∇ψ| 2 < +∞ and ψ(x) = O(|x| -1 ). The polarization tensor M = (M ij ) i,j=1,2 associated with the domain D is defined to be the 2 × 2 matrix given by (2.1)

M ij = (k -1) δ ij |D| + ∂D x j ∂ψ i ∂ν - ds , i, j = 1, 2,
where ψ i , i = 1, 2, is the unique solution in W 1,2 to the following transmission problem:

(2.2)

∆ψ i = (1 -k)∇ • (χ(D) (∇ψ i + e i )) , or equivalently,                  ∆ψ i = 0 in (R 2 \ D) ∪ D, ψ i | + -ψ i | -= 0 on ∂D, ∂ψ i ∂ν + -k ∂ψ i ∂ν - = (k -1)ν i on ∂D, ψ i (x) = O(|x| -1 ) as |x| → ∞.
Here, χ(D) is the characteristic function of D, e i is the standard basis of R 2 , the subscripts + and -denote the limits from outside and inside D, respectively, and ν = (ν 1 , ν 2 ) is the unit outward normal to ∂D. Then one can show (see [START_REF]A review of some recent work on impedance imaging for inhomogneieties of low volume fraction[END_REF]) that for any ξ ∈ R 2

(2.3) M ξ • ξ = R 2 γ D ∇w + k -1 k χ(D)ξ 2 + k -1 k |ξ| 2 |D|,
where w is the solution to

∆w = (1 -k)∇ • (χ(D) (∇w + ξ)) .
Here and throughout this paper γ D := 1 + (k -1)χ(D). Moreover, the solution w is the minimizer of the functional

R 2 γ D ∇φ + k -1 k χ(D)ξ 2 + k -1 k |ξ| 2 |D| over φ ∈ W 1,2 . It then follows that (2.4) R n γ D ∇w + χ(D) k -1 k ξ 2 = - R 2 γ D |∇w| 2 + (k -1) 2 k |ξ||D|.
See [START_REF]A review of some recent work on impedance imaging for inhomogneieties of low volume fraction[END_REF] for details. The definition (2.1) of PT is proved to be equivalent to the following representation using layer potential techniques [START_REF]High-order terms in the asymptotic expansions of the steady-state voltage potentials in the presence of conductivity inhomogeneities of small diameter[END_REF][START_REF] Ammari | Reconstruction of Small Inhomogeneities from Boundary Measurements[END_REF]. For φ ∈ L 2 (∂D), let K * D φ be defined by

K * D φ(x) = 1 2π p.v. ∂D x -y, ν x |x -y| 2 φ(y)dσ(y)
, where p.v. stands for the Cauchy principal value. Then M can be represented as follows:

(2.5)

M ij = ∂D y j (λI -K * D ) -1 (ν i )(y)dσ(y), λ := k + 1 2(k -1)
.

This boundary integral representation of PT is particularly useful for computational purpose and we will use this representation for all the computations provided in this paper.

Parameterizations of domains and attainability

In this section, we show numerically that every point inside the bounds (1.3) and (1.4) is a pair of eigenvalues of a PT associated with a simply connected domain of unit area. The domains varies from ellipses to crosses as shown in the introduction. The crosses to be considered in this section will be slightly different from those considered in the previous section in order to avoid the corners.

Let us briefly explain how we numerically compute the PT M . The domain D to be considered is symmetric with respect to x and y axes, and centered at the origin (0, 0). In order to compute M using its definition given in (2.5), we compute (λI -K * D ) -1 (ν i )(y) in the following way. We first choose 2m points x 1 , x 2 , . . . , x 2m on ∂D, and find φ i , i = 1, 2, satisfying

(λI -K * D )φ i (x k ) = ν i (x k ), k = 1, 2, • • • , 2m
. This equation is solved by making use of the collocation method and trapezoid rule (see, for instance, [START_REF] Kress | Linear Integral Equations[END_REF]). Because of the symmetry of D, the discretized operator

(λI -K * D ) is a matrix of the form A B B A
where A and B are m × m matrices, and the vector

ν i = ν i (x k ) for k = 1, . . . , 2m, takes the form a i -a i . Therefore, the solution φ i = (φ i (x k )) 2m k=1 should be of the form b i -b i , where the vector b i ∈ R m is the solution of the algebraic equation (3.1) (A -B)b i = a i , i = 1, 2.
In order to solve (3.1), we use the decomposition P (A -B) = LU where the matrix L, U and P are a lower triangular, an upper triangular, and a permutation matrix, respectively. Once (3.1) is solved, we can in a standard way use (2.5) to finally calculate M . As it will be shown in the following, the domain D can be a thin long cross. As the cross gets thinner and longer, we need more points on ∂D, and hence the matrix A -B becomes larger, causing more difficulties for solving (3.1). But, thanks to the symmetry of D, A -B is a sparse matrix and then we can decompose it into several symmetric matrices of small sizes.

Note that, because of an obvious symmetry reason, it suffices to consider the region λ 2 ≥ λ 1 . In fact, we recall that if two domains D and D are related by D = R(D) where R is a unitary transformation, then their corresponding PTs, M and M , satisfy the relation:

(3.2) M = RM R T .
For t ≥ 1, let E(t) be the ellipse of the form

x 2 a 2 + y 2 b 2 = 1, a ≤ b
D(t) := (x, y), (x, -y), (-x, y), (-x, -y)|(x, y) ∈ D ,
where |D(t)| = 1, and for x ≥ 0 and y ≥ 0, t := r 1 /r 2 , and the set D is given by

D =                                          (r 1 , y) for 0 ≤ y ≤ r2 2 , (r 1 - r 2 2 + r 2 2 x x 2 + y 2 , r 2 2 + r 2 2 y x 2 + y 2 ) for 0 ≤ tan -1 y x ≤ π 2 , (x, r 2 ) for 3r2 2 ≤ x ≤ r 1 -r2 2 , ( 3r 2 2 + r 2 2 x x 2 + y 2 , 3r 2 2 + r 2 2 y x 2 + y 2 ) for π ≤ tan -1 y x ≤ 3π 2 , (r 2 , y) for 3r2 2 ≤ y ≤ r 1 -r2 2 , ( r 2 2 + r 2 2 x x 2 + y 2 , r 1 - r 2 2 + r 2 2 y x 2 + y 2 ) for 0 ≤ tan -1 y x ≤ π 2 , (x, r 1 ) for 0 ≤ x ≤ r 2 2 .
The domain D(t) is a cross whose corners are polished to be arcs. We make the corners smooth to make the computations easier. Figure 3 shows the shape of D(t = 3). The classes of domains E(t) and D(t) are the basic building blocks for representing the optimal shapes constructed in this paper.

It is known that the eigenvalue pair of the PT of the ellipse E(t) lies on the lower bound (hyperbola) of the region A. In fact, according to [START_REF] Brühl | A direct impedance tomography algorithm for locating small inhomogeneities[END_REF], the PT for the ellipse E(t) is given by

(3.3) M = |E| (k-1)(a+b) a+kb 0 0 (k-1)(a+b) ka+b .
Thus the eigenvalue pair of the PT for On the other hand, since the domain D(t) is invariant under the rotation by π/2, by the relation (3.2) the eigenvalue pair of its PT lies on the line

E(t) is (k -1)(1 + t) 1 + kt , (k -1)(1 + t) k + t . 0 r 1 r 1 -r 1 -r 1
λ 1 = λ 2 . Let U (s, 2) = (1 -s)D(2) + sE(2), 0 ≤ s ≤ 1.
As shown in Figure 4, the eigenvalue pair of the PT associated with the domain U (s, 2) divides the region A into two parts A 1 and A 2 . Therefore, we will use different parameterizations of domains to attain points in A 1 and A 2 . This is solely for the convenience of the computations.

Attainability of the region

A 1 . Let U (s, t) = (1 -s)D(t) + sE(t), 0 ≤ s ≤ 1, 2 ≤ t.
By this notation, we do not mean the convex combination of the two shapes. Notice for every t, that D(t) and E(t) are both star-shaped domains around the origin, i.e., there exist two functions 2 dϕ is not a priori equal to one. We therefore rescale the polarization tensor by a factor 1/|U (s, t)|.

d t ∈ C 1 ([0, 2π]) and e t ∈ C ∞ ([0, 2π]) such that ∂D(t) = {(ϕ, d t (ϕ)), ϕ ∈ [0, 2π]} and ∂E(t) = {(ϕ, e t (ϕ)), ϕ ∈ [0, 2π]}. The precise definition of U (s, t) for a given s ∈ [0, 1] is the bounded domain limited by the C 1 curve ∂U (s, t) = {(ϕ, (1 -s)d t (ϕ) + se t (ϕ)), ϕ ∈ [0, 2π]}. Note that |U (s, t)| = 1 2 2π 0 ((1 -s)d t (ϕ) + se t (ϕ))
Our goal in this section is to show that the points inside A 1 are attained by the domain U (s, t). Note that U (1, t) = E(t) is an ellipse and U (0, t) = D(t) is a cross. As t gets larger, the ellipse U (1, t) becomes thinner and the cross U (0, t) becomes thinner and longer. For a fixed t, U (s, t) changes its shape from a cross to an ellipse as s varies from 0 to 1. The left-hand side of Figure 5 shows the variation of the domain U (s, t). Numerical computations show that as t → ∞, the eigenvalue pair of PT associated with U (0, t) moves toward the point (λ 1 , λ 2 ) = ((k -1)(k + 1)/(2k), (k -1)(k + 1)/(2k)) following the line λ 1 = λ 2 while that of U (1, t) moves toward the point ((k -1)/k, k -1) following the hyperbola. Both of these facts will be rigorously established in section 5. For each fixed t, the eigenvalue pair for U (s, t), 0 ≤ s ≤ 1, makes a curve in the region A 1 , and these curves make foliation of the region A 1 as t varies. Figure 5 shows variations of U (s, t) and the result of numerical computations.

Figure 5. The variation of U (s, t) and the corresponding eigenvalue pairs. As t → ∞, the cross and ellipse become thinner and longer, and corresponding eigenvalues move toward the upper bound. As s → 1, the cross becomes an ellipse and corresponding eigenvalues move toward the lower hyperbola.

Attainability of the region

A 2 . Let L(s, t) := (1 -s)E(t) + sD(2), 0 ≤ s ≤ 1 and 1 ≤ t < 2.
As above, this notation does not mean the convex combination of E(t) and D(2) but the bounded domain limited by the curve {(ϕ, (1 -s)e t (ϕ) + sd 2 (ϕ)), ϕ ∈ [0, 2π]}. Note that D(2) is fixed. For each fixed t, E(t) corresponds to a point on the lower hyperbola. So, L(s, t) interpolates the point on the hyperbola and the fixed point on the line λ 1 = λ 2 corresponding to D [START_REF]High-order terms in the asymptotic expansions of the steady-state voltage potentials in the presence of conductivity inhomogeneities of small diameter[END_REF]. Numerical computations give the evidence of the fact that the points in A 2 arise as the eigenvalue pairs of the PT corresponding to the domain L(s, t). Figure 6 shows the variation of L(s, t) and the corresponding eigenvalues of its PT. Figure 7 shows attainability of both the regions A 1 and A 2 .

An approximate formula for the foliation

In the previous section, we parameterized domains to show that as the parameter s varies from 0 to 1 for a fixed t, the eigenvalue pairs make a curve. Then, by varying t these curves, a foliation of the region, say A 1 , can be obtained.

(2) Figure 6. The variations of L(s, t) and the corresponding eigenvalue pairs of its PT. In this section, we derive an approximate formula for this foliation. Since the formula for the region A 2 can be obtained in almost the same way, we omit its derivation.

Choose finite sets of parameters, S = {s 0 = 0, s 1 , . . . , s m = 1} and T = {t 0 = 2, t 1 , . . . , t n }. Then numbers m, n are at our disposal. We may take large ones if we need very accurate approximations. We then compute the eigenvalues (λ 1 (s i , t j ), λ 2 (s i , t j )) of the PT associated with the domain U (s i , t j ) for i = 0, . . . , m and j = 0, . . . , n. Our objective is to find a formula to relate λ 1 (s i , t j ) with λ 2 (s i , t j ).

Let us introduce a transformation which makes the computations easier. Let V be given by

V (x, y) := R π 4 x - (k -1)(k + 1) 2k , y - (k -1)(k + 1) 2k ,
where R θ is the rotation by θ. Then define

(ρ 1 (s, t), ρ 2 (s, t)) := V (λ 1 (s, t), λ 2 (s, t)).
We seek to find a function P (x, t) of the form

P (x, t) = m l=0 a l (t)x l , such that ρ 2 (s i , t j ) = P (ρ 1 (s i , t j ), t j ) = m l=0 a l (t j )ρ 1 (s i , t j ) l .
We can obtain the desired formula by completing the following steps.

An approximate formula for the foliation S1. For each s i ∈ S (i = m) and l = 1, 2, we look for a function P l i on (1, 3/2] satisfying

P l i 1 + t j t j = ρ k (s i , t j ), t j ∈ T,
or minimizing

tj ∈T P l i ( 1 + t j t j ) -ρ k (s i , t j ) .
We obtain such a polynomial P l i using the least square curve fitting. If i = m or s i = 1, then U (s m , t) is an ellipse, and the corresponding eigenvalues are given explicitly by (3.3). We use these explicit formulas for our P l m , l = 1, 2. S2. We seek a function P (x, t), which is the unique m-th degree polynomial in x satisfying P (ρ 1 (s i , t j ), t j ) = ρ 2 (s i , t j ), s i ∈ S and t j ∈ T. Such a polynomial is given by

P (x, t) = m i=0 a i (t)x i = m i=0 l =i l=0,...,m ρ 2 (s i , t) x -ρ 1 (s l , t) ρ 1 (s i , t) -ρ 1 (s l , t) = m i=0 l =i l=0,...,m P 2 i 1 + t t x -P 1 l 1+t t P 1 i 1+t t -P 1 l 1+t t
, where P 1 i , P 2 i , i = 0, . . . , m, are the polynomials obtained in step S1.

Notice that for the sake of a simpler parameterization we have used two different ways of interpolation for the lower and upper parts.

Figure 8 shows results of numerical tests using Matlab program. For these computations we assume the conductivity k = 3. The boundary integral is discretized at 16, 000 boundary points. This large number of points will be used throughout the paper. After doing experiments on the accuracy of the method, we have chosen this number to make the numerical computations precise, since as the cross gets thinner, the domain variations affect very little the PT. The dots in the left figure of Figure 8 are the computed eigenvalue pairs of the PT associated with U (s, t) for s = 0, 0.2, • • • , 1 and t = 2, 3, • • • , 18. The solid-curves are graphs of P (x, t) fitting dots, which are computed following the procedure described in steps S1 and S2. We use S := {0, 0.4, 1} and T := {2, 3, • • • , 18} for the computation of P (x, t). The right hand side figure is for L(s, t), s = 1, 1.1, . . . , 1.9, t = 2, 3, • • • , 18. Numerical values for the eigenvalues are recorded in Tables 1 and2. Figure 9 shows whole foliation by P (x, t) of the HS bounds. Table 2. The eigenvalue pairs (λ L 1 (s, t), λ L 2 (s, t)) associated with L(s, t).

Asymptotics of PT for thin crosses, and analytic accuracy of the foliation

On Figure 5 we see that the upper trace bound is reached by the PT corresponding to E(t) as t tends to infinity. It also seems that as the cross becomes thinner and longer, the corresponding PT approaches the upper trace bound.

The first result of this section show that this is indeed what is predicted by the asymptotics.

V -1

M D(t) = (k -1)(k + 1) 2k
1 0 0 1 .

On Figure 9, it seems that as the "thinness" parameter t becomes large, the curve described by (λ U 1 (s, t), λ U 2 (s, t)) as s varies between 0 and 1 appears to be a straight line. This can be measured on table 1 as well. The cross-shape was designed to be "simple" from a computational point of view. The drawback is that obtaining precise analytical results for this shape is unreasonably technical. We shall show that the affine nature of the curve (λ U 1 (s, t), λ U 2 (s, t)) is asymptotically correct for a similar (but different) cross for which analytical computations are simpler.

The alternate cross we consider is given as follows. Let F t and S t be a " flat" elliptic beam and a "standing" elliptic beam given by

F t : x 2 b 2 + y 2 a 2 = 1, S t : x 2 a 2 + y 2 b 2 = 1, a = 1 √ 2πt , b = ta.
Define X t := F t ∪ S t , and following the approach introduced in section 3, define for all s ∈ [0, 1] and t ∈ [1, ∞),

U (s, t) = sE(t) + (1 -s)X(t),
the domain of unit area obtained by radial interpolation between X(t) and E(t).

The second result of this section shows that as s varies between 0 and 1, the curve described by (λ U 1 (s, t), λ U 2 (s, t)) is asymptotically a line. Theorem 5.2. The PT M U (s,t) associated to U (s, t) satisfies

M U (s,t) = √ 2s + 1 -s 2 √ 2s + 1 -s 2 + (1 -s) 2 M U (1,t) + (1 -s) 2 √ 2s + 1 -s 2 + (1 -s) 2 M U (0,t) + (t),
where (t) tends to zero as t tends to infinity.

We shall start by proving Theorem 5.1. Note that D(t) can be approximated by Dt = H t ∪ V t , where

H t = - t 8 , t 8 × - 1 √ 8t , 1 √ 8t and V t = - 1 √ 8t , 1 √ 8t × - t 8 , t 8 .
The difference between the two domains are the smooth corners of D(t), and the surface of H t ∪ V t which is not exactly one. It is easy to see that, for every t > 0 there exists t 1 and t 2 and λ such that

|t 1 -t 2 | + |λ -1| ≤ C t ,
where C is a constant independent of t, and Dt 1 ⊂ D(t) ⊂ λ Dt 2 . Since the polarization tensors are monotonous with respect to inclusion, we deduce that

M D(t) = M D(t) + O( 1 t
).

In the sequel, we will identify D and D(t). The asymptotic limit of the polarization tensor of a single thin beam is known [START_REF] Beretta | Asymptotic formulas for steady state voltage potentials in the presence of thin inhomogeneities. A rigorous error analysis[END_REF]. We have

lim t→∞ M H(t) = 1 2 k -1 0 0 k -1 k , and lim t→∞ M V (t) = 1 2 k -1 k 0 0 k -1 . Given ξ ∈ R 2 , let w ξ t ∈ W 1,2 be the solution to ∆w t = (1 -k)∇ • (χ(D(t)) (∇w t + ξ)) ,
and let h ξ t , and v ξ t be the solutions to the same equations with D(t) replaced with H t and V t respectively.

The following lemmas hold.

Lemma 5.3. For any unit vector ξ = (ξ 1 , ξ 2 ), the functions h ξ t and v ξ t satisfy (5.2)

∇h ξ t + χ(H t ) 1 -k k ξ 2 e 2 L 2 (R 2 )
= (t), and

∇v ξ t + χ(V t ) 1 -k k ξ 1 e 1 L 2 (R 2 ) = (t),
where lim t→∞ (t) = 0.

Proof. Write M (H t ) = (M i,j t ) i,j=1,2 . Since M

= (t).

Since h ξ t = ξ 1 h e 1 t + ξ 2 h e 2 t , we obtain the first limit in (5.2). The second estimate can be proved in a similar fashion. Lemma 5.4. For each unit vector ξ, we have

(5.3) ∇w ξ t -∇h ξ t + ∇v ξ t L 2 (R 2 ) = (t),
where lim t→∞ (t) = 0.

Proof. Let us, for simplicity, drop superscripts ξ. Since

∇h t L 2 (H t ∩V t ) ≤ ∇h t + χ(H t ) 1 -k k ξ 2 e 2 L 2 (H t ∩V t ) + ξ 2 |H t ∩ V t | = (t),
it follows from (5.2) that (5.4) ∇h t L 2 (Vt) = (t).

Likewise we get (5.5)

∇v t L 2 (H t ) = (t).
By simple algebra, we obtain

∇ • γ D(t) ∇ (w t -h t -v t ) = (k -1) ∇ • (χ(V t )∇h t + χ(H t )∇v t ) -(k -1) ∇ • (χ(H t ∩ V t ) (∇h t + ∇v t + ξ)) .
If we multiply the above equation by a test function φ and integrate by parts, we get

R 2 γ Xt ∇ (w t -h t -v t ) • ∇φ = (k -1) (A + B + C) with |A| = H t ∩V t (∇h t + ∇v t + ξ) • ∇φ ≤ ∇h t + ∇v t + ξ L 2 (Ht∩Vt) ∇φ L 2 (R 2 ) , |B| = V t ∇h t • ∇φ ≤ ∇h t L 2 (V t ) ∇φ L 2 (R 2 ) , |C| = Ht ∇v t • ∇φ ≤ ∇v t L 2 (Ht) ∇φ L 2 (R 2 ) .
Thanks to (5.4) and (5.5), we obtain

R 2 γ D(t) ∇ (w t -h t -v t ) • ∇φ ≤ ∇φ L 2 (R 2 ) (t),
for any test function φ, and hence (5.3) follows.

We are now ready to prove Theorem 5.1.

Proof of Theorem 5.1. By (2.1) we have

M D(t) ξ • ξ = (k -1) 1 + (k -1) D(t) ξ • ∇w ξ t dx , for any ξ ∈ R 2 . It then follows from (5.3) that (5.6) M D(t) ξ • ξ = (k -1) 1 + (k -1) D(t) ξ • (∇h ξ a + ∇v ξ d )dx + (t).
From (5.2) we obtain that

D(t) ξ • (∇h ξ t + ∇v ξ t )dx = H t ξ • ∇h ξ t + V t ξ • ∇v ξ t + D(t)\H t ξ • ∇h ξ t + D(t)\V t ξ • ∇v ξ t = H t ξ • ∇h ξ t + V t ξ • ∇v ξ t + (t).
Therefore, (5.6) yields

M D(t) ξ • ξ = M V (t) ξ • ξ + M H(t) ξ • ξ + (t).
Now, (5.1) immediately follows from [START_REF] Ammari | Reconstruction of closely spaced small inclusions[END_REF].

Proof of Theorem 5.2. We have introduced

U (s, t) = sE(t) + (1 -s)X(t),
that is, a domain with vertical and horizontal symmetries defined in the first quadrant as the domain inside the curve given, in polar coordinates, by

for 0 ≤ ϕ ≤ π 4 , ρ U (ϕ) = t 2π *   √ 2s t 2 cos 2 (ϕ) + sin 2 (ϕ) + 1 -s t 2 sin 2 (ϕ) + cos 2 (ϕ)   for π 4 ≤ ϕ ≤ π 2 , ρ U (ϕ) = t 2π * √ 2s + 1 -s t 2 cos 2 (ϕ) + sin 2 (ϕ)
.

With this definition, we can compute

|U (s, t)| = π 4 0 1 2 ρ 2 U (ϕ)dϕ + π 2 π 4 1 2 ρ 2 U (ϕ)dϕ = (1 -s) 2 2 + √ 2s + 1 -s 2 2 + (t).
Furthermore, we see that, up to a small volume, U (s, t) = Fst ∪ Sst + (t), where Fst :

x 2 ã2 + y 2 b2 = 1, Sst : x 2 c2 + y 2 d2 = 1, a = √ 2s + 1 -s √ 2πt , b = ta and c = td, d = (1 -s) √ 2πt .
Arguing as before, we see that the PT corresponding to U (s, t) is equal to the PT corresponding to Fst ∪ Sst , up to an error (t) which tends to zero as t tends to infinity. The polarization tensors of Fst and Sst are known [START_REF] Brühl | A direct impedance tomography algorithm for locating small inhomogeneities[END_REF] and given by (3.3). Thus we have

lim t→∞ M F st = F st k -1 0 0 k -1 k , and lim t→∞ M S(t) = S st k -1 k 0 0 k -1 .
Using the same arguments that in the proof of Theorem 5.1, with D(t), H(t) and V (t) replaced by U (s, t), Fst and Sst respectively, formula (5) becomes

U (s, t) M U (s,t) ξ • ξ = M F st ξ • ξ + M S st ξ • ξ + (t),
which, together with (5), concludes the proof of Theorem 5.2.

In concluding this section, we would like to refer to the recent paper by Capdeboscq and Kang [START_REF] Capdeboscq | Improved Hashin-Shtrikman bound for thick domains[END_REF] where it is shown that if the domain contains certain set of volume and is not thin than its PT stays away from the upper HS-bound.

It is also worth noticing that the asymptotic result in this section holds in the three-dimensional case as well. Indeed, taking combinations of ellipsoids, we can show by exactly the same arguments as those presented above the optimality of the upper bounds (1.1) 3. The numerical results for the identification algorithm. For a given eigenvalue pair (λ 1 , λ 2 ), we compute the parameters t * , s * of the domain and find the eigenvalues (λ 1 (s * , t * ), λ 2 (s * , t * )). The last column is the relative error.

Identification of domains

In this section, given a pair (λ 1 , λ 2 ) inside the HS bounds, we present a method to find the parameters s and t, such that this pair arises as the eigenvalue pair of the PT associated with U (s, t) or L(s, t).

Our procedure is described in the following.

Identification algorithm

S'1. Given (λ 1 , λ 2 ) ∈ A 1 , transform it into (ρ 1 , ρ 2 ) := V (λ 1 , λ 2 ). S'2. Find the zero t * ≥ 2 satisfying P (ρ 1 , t * ) = ρ 2 .

Here P (x, t) is the polynomial computed in the previous section. S'3. Compute (ρ 1 (s i , t * ), ρ 2 (s i , t * )) for s i ∈ S. From these data, we derive the relation between s and ρ 1 (s, t * ) as follows:

x 1 (s) = m i=0 l =i ρ 1 (s i , t * ) s -s l s i -s l , 0 ≤ s ≤ 1.
S'4. Find s * so that

x 1 (s * ) = ρ 1 .
We can check the validity of the above algorithm by computing the relative error. If (s * , t * ) is the computed pair of parameters, then we calculate the eigenvalues λ 1 (s * , t * ) and λ 2 (s * , t * ) of its associated PT. We can see the order of the relative error by computing

|(λ 1 , λ 2 ) -(λ 1 (s * , t * ), λ 2 (s * , t * ))| |(λ 1 , λ 2 )| .
As shown in Table 3, the above algorithm identifies the parameterized domain pretty well.

Conclusion

In this paper we proved that the PT associated with crosses approaches to the upper bound as cross gets thinner and longer. We then show, by numerical computations based on a boundary integral method, that every point inside the Hashin-Shtrikman bounds for polarization tensors can be attained by a simply connected domain which is given explicitly in terms of two parameters. We also present a method to find the parameters corresponding to a given pair of eigenvalues inside the Hashin-Shtrikman bounds.

Figure 2 .

 2 Figure 2. The variation of cross-domains starting from the unit disk.

Figure 3 .Figure 4 .

 34 Figure 3. The graph of D(t) for t = r 1 /r 2 = 3.

Figure 8 .

 8 Figure 8. Numerical results for U and L. The dots are the computed eigenvalue pairs of the PT of the domains and solid curves are curves fitting dots.

Figure 9 .

 9 Figure 9. The numerical results for the foliation.

  , where t := b/a and the area |E(t)| = 1. The cross-shaped domain D(t), t ≥ 2, is defined by

Table 1 .

 1 The eigenvalue pairs (λ U

	t s	0	0.2	0.4	0.8	1
	2	(1.0111,1.0111) (0.9736,1.0463) (0.9376,1.0847) (0.8778,1.1640) (0.8571,1.2000)
	3	(1.0443,1.0443) (0.9812,1.1021) (0.9194,1.1656) (0.8227,1.2886) (0.8000,1.3333)
	4	(1.0737,1.0737) (0.9945,1.1480) (0.9160,1.2292) (0.7938,1.3803) (0.7692,1.4286)
	5	(1.0981,1.0981) (1.0074,1.1852) (0.9171,1.2798) (0.7767,1.4498) (0.7500,1.5000)
	6	(1.1182,1.1182) (1.0187,1.2158) (0.9195,1.3208) (0.7656,1.5042) (0.7368,1.5556)
	7	(1.1351,1.1351) (1.0285,1.2413) (0.9220,1.3548) (0.7578,1.5482) (0.7273,1.6000)
	8	(1.1491,1.1491) (1.0366,1.2627) (0.9243,1.3833) (0.7520,1.5844) (0.7200,1.6364)
	9	(1.1622,1.1622) (1.0443,1.2815) (0.9268,1.4078) (0.7475,1.6147) (0.7143,1.6667)
	10	(1.1717,1.1717) (1.0497,1.2970) (0.9282,1.4286) (0.7438,1.6407) (0.7097,1.6923)
	11	(1.1811,1.1811) (1.0552,1.3111) (0.9299,1.4470) (0.7410,1.6629) (0.7059,1.7143)
	12	(1.1881,1.1881) (1.0591,1.3229) (0.9308,1.4629) (0.7382,1.6824) (0.7027,1.7333)
	13	(1.1964,1.1964) (1.0640,1.3344) (0.9325,1.4773) (0.7361,1.6995) (0.7000,1.7500)
	14	(1.2025,1.2025) (1.0675,1.3441) (0.9335,1.4901) (0.7342,1.7147) (0.6977,1.7647)
	15	(1.2077,1.2077) (1.0703,1.3527) (0.9341,1.5013) (0.7325,1.7281) (0.6957,1.7778)
	16	(1.2134,1.2134) (1.0736,1.3609) (0.9351,1.5117) (0.7309,1.7403) (0.6939,1.7895)
	17	(1.2205,1.2205) (1.0779,1.3690) (0.9368,1.5213) (0.7299,1.7511) (0.6923,1.8000)
	18	(1.2272,1.2272) (1.0819,1.3764) (0.9384,1.5301) (0.7282,1.7613) (0.6909,1.8095)

1 (s, t), λ U 2 (s, t)) associated with U (s, t).

  and (1.2).(λ 1 , λ 2 ) t * s * (λ 1 (s * , t * ), λ 2 (s * , t * )) |(λ1,λ2)-(λ1(s * ,t * ),λ2(s * ,t * ))| |(λ 1 ,λ 2 )|

	(1.0470,1.2893) 9.5607	0.1992	(1.0480,1.2902)	8.0838E-4
	(0.8277,1.4614) 6.8182	0.6081	(0.8248,1.4586)	2.3709E-3
	(1.2060,1.2060) 14.4493 4.7730E-15	(1.2062,1.2062)	2.1222E-4
	(0.9387,1.5378) 19.0739	0.4008	(0.9403,1.5394)	1.2567E-3
	(0.7274,1.7702) 18.4077	0.8089	(0.7253,1.7681)	1.5114E-3
	(1.1015,1.4132) 26.7781	0.1981	(1.1065,1.4182)	3.9419E-3
	Table

Acknowledgements

The question we consider in this paper was raised to us by Michael Vogelius. We acknowledge it with gratitude. H. Kang and E. Kim are partly supported by grant R02-2003-000-10012-0 from the Korea Science and Engineering Foundation. M. Lim is partly supported by the post-doctoral fellowship program of the Korea Science and Engineering Foundation.


A Landragin

F Pereira

Dos Santos

Accelerometer using atomic waves for space applications

The techniques of laser cooling combined with atom interferometry make possible the realization of very sensitive and accurate inertial sensors like gyroscopes or accelerometers. Besides earth-based developments, the use of these techniques in space should provide extremely high sensitivity for research in fundamental physics, Earth's observation and exploration of the solar system.

-Introduction

Inertial sensors are useful devices in both science and industry. Higher precision sensors could find scientific applications in the areas of general relativity [1], geodesy and geology. There are also important applications of such devices in the field of navigation, surveying and analysis of earth structures. Matter-wave interferometry was envisaged for its potential to be an extremely sensitive probe for inertial forces [2]. In 1991, atom interference techniques have been used in proof-of-principle work to measure rotations [3] and accelerations [4]. In the following years, many theoretical and experimental works have been performed to investigate this new kind of inertial sensors [START_REF] Berman | Atom Interferometry[END_REF]. Some of the recent works have shown very promising results leading to a sensitivity comparable to c Società Italiana di Fisica 1 other kinds of sensors, as well as for rotation [START_REF] Gustavson | [END_REF] as for acceleration [7] and possibility of realizing a full inertial base within the same device [8]. The most developed atominterferometer inertial sensors are today atomic state interferometers [9] which use twophoton velocity selective Raman transitions [10,11] to manipulate atoms while keeping them in long-lived ground states.

Atom interferometry is nowadays one of the most promising candidates for ultraprecise and ultra-accurate measurement of gravito-inertial forces on ground or for space [START_REF]Special Issue: Quantum Mechanics for Space Application: From Quantum Optics to Atom Optics and General Relativity[END_REF]. From performances on ground, one can expect unprecedented sensitivity in space, leading to many mission proposals since 2000 [START_REF]Hyper-Precision Cold Atom Interferometry in Space (HYPER)[END_REF]. This technology is now mature enough that several groups are developing instruments for practical experiments: in the field navigation [START_REF] Kasevich | [END_REF], or fundamental physics (gradiometer for the measurement of G [START_REF] Bertoldi | [END_REF], gravimeter for the watt balance experiment [16], interferometer for the measurement of fine structure constant thanks to h/m [17]). Moreover, the realization of Bose-Einstein condensation (BEC) of a dilute gas of trapped atoms in a single quantum state [18,19,20] has produced the matter-wave analog of a laser in optics [21,22,23,24] and open new possibilities. Alike the revolution brought by lasers in optical interferometry [1,25,26], it is expected that the use of Bose-Einstein condensed atoms will bring the science of atom optics, and in particular atom interferometry, to an unprecedented level of accuracy [27,28]. In addition, BEC-based coherent atom interferometry would reach its full potential in space-based applications where micro-gravity will allow the atomic interferometers to reach their best performance. Applications of accelerometer in space concern fundamental physics, like testing the equivalence principle by comparing the free fall of two different atomic species [START_REF]Proposal to Cosmic vision[END_REF], and more generally testing all aspects of gravity as in [START_REF] Wolf | Submitted to publication to Experimental Astronomy[END_REF], which aims at testing the laws of gravity at large scale, as well for fundamental physics as for exploration of the solar system.

In the following part, we investigate the sources of noise and systematic effects for such atom interferometers, thanks to the results obtained with the gravimeter under development in our laboratory. This setup is based on same concepts than previous works [7]: atomic source realized with a magneto-optical traps and manipulation of the atomic wavepackets with Raman transitions. In our experiment, we have studied in details the influence of any perturbations on the sensitivity of the sensor. Reducing there impacts on the noise of the instrument, we finally reach an excellent sensitivity of 1.4 × 10 -8 m.s -2 .Hz -1/2 , despite a rather short interrogation time (100 ms only). We think this experiment is a good benchmarck to oversee the performances of best space accelerometers based on same technologies. Performances in space environment will be derived, taking into account the specifications of the environment and the much longer interrogation time.

-Experimental setup

The principle of our gravimeter is based on the coherent splitting of matter-waves by the use of two-photon Raman transitions [START_REF] Kasevich | [END_REF]. These transitions couple the two hyperfine levels F = 1 and F = 2 of the 5 S 1/2 ground state of the 87 Rb atom. An intense beam of slow atoms is first produced by a 2D-MOT. Out of this beam 10 7 atoms are loaded within 50 ms into a 3D-MOT and subsequently cooled in a far detuned (-25 Γ) optical molasses. The lasers are then switched off adiabatically to release the atoms into free fall at a final temperature of 2.5 µK. Both lasers used for cooling and repumping are then detuned from the atomic transitions by about 1 GHz to generate the two off-resonant Raman beams. For this we have developed a compact and agile laser system that allows us to rapidly change the operating frequencies of these lasers, as described in [32]. Before entering the interferometer, atoms are selected in a narrow vertical velocity distribution (σ v ≤ v r = 5.9 mm/s) in the |F = 1, m F = 0 state, using a combination of microwave and optical Raman pulses.

The interferometer is created by using a sequence of three pulses (π/2 -π -π/2), which split, redirect and recombine the atomic wave packets. Thanks to the relationship between external and internal state [9], the interferometer phase shift can easily be deduced from a fluorescence measurement of the populations of each of the two states. Indeed, at the output of the interferometer, the transition probability P from one hyperfine state to the other is given by the well-known relation for a two wave interferometer:

P = 1 2 (1 + C cos ∆Φ)
, where C is the interferometer contrast, and ∆Φ the difference of the atomic phases accumulated along the two paths. The difference in the phases accumulated along the two paths depends on the acceleration a experienced by the atoms. It can be written as [START_REF] Kasevich | [END_REF] 2 , where φ(0, T, 2T ) is the difference of the phases of the lasers, at the location of the center of the atomic wavepackets, for each of the three pulses [33]. Here k ef f = k 1 -k 2 is the effective wave vector (with | k ef f | = k 1 + k 2 for counter-propagating beams), and T is the time interval between two consecutive pulses.

∆Φ = φ(0) -2φ(T ) + φ(2T ) = -k ef f • aT
The Raman light sources are two extended cavity diode lasers based on the design of [34], which are amplified by two independent tapered amplifiers. Their frequency difference is phase locked onto a low phase noise microwave reference source. The two overlapped beams are injected in a polarization maintaining fiber, and guided towards the vacuum chamber. We obtain counter-propagating beams by placing a mirror and a quarterwave plate at the bottom of the experiment. Four beams are actually sent onto the atoms, out of which only two will drive the counter-propagating Raman transitions, due to conservation of angular momentum and the Doppler shift induced by the free fall of the atoms.

Experimental setups, based on the same principle, can be realized for space experiments. They will benefit from the technical developments realized in the frame of the PHARAO atomic clock [35] for the space ACES project. -The sensitivity function characterizes the influence of fluctuations in the Raman lasers phase difference φ onto the transition probability [START_REF] Dick | Proc. Nineteenth Annual Precise Time and Time Interval[END_REF], and thus on the interferometer phase. This function is defined by :

g(t) = 2 lim δφ→0 δP (δφ, t) δφ . ( 1
)
where δφ is a jump on the Raman phase difference φ, which occurs at time t during the interferometer sequence, and induces a change of δP (δφ, t) in the transition probability.

The expression of the sensitivity function can easily be derived when considering that the Raman pulses are infinitesimally short. In that case, the interferometer phase Φ is given by [4]: Φ = φ 1 -2φ 2 + φ 3 , where φ 1 ,φ 2 ,φ 3 are the the Raman laser phase differences at the three laser interactions, taken at the position of the center of the atomic wavepacket [33]. Usually, the interferometer is operated at mid fringe (Φ = π/2), in order to maximize the sensitivity to interferometer phase fluctuations. If the phase step δφ occurs for instance between the first and the second pulses, the interferometric phase changes by δΦ = -δφ, and the transition probability by δP = -cos(π/2+δΦ)/2 ∼ -δφ/2 in the limit of an infinitesimal phase step. Thus, in between the first two pulses, the sensitivity function is -1. The same way, one finds for the sensitivity function between the last two pulses : +1.

In the more general case of finite duration Raman laser pulses, the sensitivity function will depend on the time evolution of the atomic state during the pulses. This function is calculated in [START_REF] Cheinet | [END_REF] considering laser waves as plane waves and quantizing atomic motion in the direction parallel to the laser beams, in the case of a constant Rabi frequency (square pulses) and resonance condition fulfilled. We calculated the change in the transition probability for a infinitesimally small phase jump at any time t during the interferometer, and deduce g(t).

The sensitivity function is an odd function, whose expression is given here for t > 0:

g(t) =    sin(Ω R t) 0 < t < τ R 1 τ R < t < T + τ R -sin(Ω R (T -t)) T + τ R < t < T + 2τ R (2)
where we choose the time origin at the middle of the second Raman pulse and where Ω R /2π is the Rabi frequency.

Using this function, we can now evaluate the fluctuations of the interferometric phase Φ for an arbitrary Raman phase noise φ(t) on the lasers

δΦ = +∞ -∞ g(t)dφ(t) = +∞ -∞ g(t) dφ(t) dt dt. (3)

3

.

2. Influence of the phase noise onto the sensitivity of the interferometer . -The sensitivity of the interferometer is characterized by the Allan variance of the interferometer phase fluctuations, σ 2 Φ (τ ), defined as

σ 2 Φ (τ ) = 1 2 ( δΦ k+1 -δΦ k ) 2 (4) = 1 2 lim n→∞ 1 n n k=1 ( δΦ k+1 -δΦ k ) 2 . ( 5
)
where δΦ k is the average value of δΦ over the interval [t k , t k+1 ] of duration τ . The Allan variance is equal, within a factor of two, to the variance of the differences in the successive average values δΦ k of the interferometric phase. As the interferometer is operated sequentially at a rate

f c = 1/T c , τ is a multiple of T c : τ = mT c .
In order to evaluate correctly the stability of the interferometer phase Φ, it is necessary to take into account that the measurement is pulsed. The sensitivity of the interferometer is limited by an aliasing phenomenon, similar to the Dick effect for atomic clocks [START_REF] Dick | Proc. Nineteenth Annual Precise Time and Time Interval[END_REF]: only the phase noise at multiples of the cycling frequency contribute to the Allan variance, weighted by the Fourier components of the transfer function. For large averaging times τ , the Allan variance of the interferometric phase is given by

σ 2 Φ (τ ) = 1 τ ∞ n=1 |H(2πnf c )| 2 S φ (2πnf c ) (6)
where S φ (ω) is the power spectral density of the Raman phase, and H(ω) is the transfer function, given by H(ω) = ωG(ω), where G is the Fourier transform of the sensitivity function.

G(ω) = +∞ -∞
e -iωt g(t)dt (7) At low frequency (ω << Ω R ), G can be approximated by

G(ω) = - 4i ω sin 2 [ω(T + 2τ R )/2] (8)
The transfer function |H 2 | has two important features. First, it presents oscillations at a frequency given by 1/(T + 2τ R ), leading to zeros at frequencies given by f k = k T +2τR . The second is a low pass first order filtering due to the finite duration of the Raman pulses.

For white phase noise S 0 φ , and to first order in τ R /T , the phase stability is given by:

σ 2 Φ (τ ) = πΩ 2 S 0 φ T c τ (9)
where Ω is the Rabi frequency.

This illustrates the filtering of the transfer function: the shorter the pulse duration τ R , the greater the interferometer noise.

3

.

3. The 100 MHz source oscillator . -To reduce the noise of the interferometer, the frequency difference between the Raman beams needs to be locked to a very stable microwave oscillator, whose frequency is close to the hyperfine transition frequency (ν mw = 6.834 GHz for 87 Rb). The reference frequency will be delivered by a frequency chain, which transposes in the microwave domain the low phase noise of an RF oscillator, typically a quartz oscillator. When this transposition is performed without degradation, the phase noise power spectral density of the RF oscillator, of frequency ν rf , is multiplied by (ν mw /ν rf ) 2 .

No single quartz oscillator fulfills the requirements of ultra low phase noise over a sufficiently large frequency range. We present in figure 1 the specifications for different high stability quartz : a Premium 10 MHz-SC from Wenzel, a BVA OCXO 8607-L from Oscilloquartz, and a Premium 100 MHz-SC quartz from Wenzel. The phase noise spectral density is displayed at the frequency of 100 MHz, so that one can compare the different quartz oscillators. A 100 MHz source for a space interferometer could be realized by combining two quartz oscillators. A 100 MHz quartz would be locked onto one of the above mentioned high stability 10 MHz reference oscillators. The bandwidth of this lock will correspond to the frequency below which the phase noise of the reference quartz is lower than the noise of the 100 MHz quartz.

The phase noise properties of such a combined source can be seen in figure 1, where we display as a solid line the phase noise spectral density of the 100 MHz source developed by THALES for the PHARAO space clock project. This combined source has been optimized for mimimal phase noise at low frequency, where it reaches a level of noise lower than any commercially available quartz. An atomic clock is indeed mostly limited by the lower frequency part of the frequency spectrum, so the requirements on the level of phase noise at higher frequency (f > 1kHz) is less stringent than for an atom interferometer. A medium performance 100 MHz oscillator is thus sufficient.

Using a simple model for the phase lock loop, we calculated the phase noise spectral density of the different combined sources we can realize by locking the Premium 100 MHz-SC either on the Premium 10 MHz-SC (Source 1), or on the BVA (source 2), or even on the Pharao source (source 3). We then estimated the impact on the interferometer of the phase noise of the 100 MHz source, assuming we are able to transpose the performance of the source at 6.8 GHz without degradation. We calculated using 6 the Allan standard deviation of the interferometric phase fluctuations for the different configurations and for various interferometer parameters. The results are presented in table I Table I. -Contribution of the 100 MHz source phase noise to the interferometrer phase fluctuations. The calculation has been performed for a Rb interferometer, for each of the three different sources and for various parameters of the interferometer. First case corresponds to the usual parameters of our gravimeter, and second (resp. third) case corresponds to typical values for space experiments with a MOT atomic source (resp. ultra-cold source).

For short interrogation times (such as 2T = 100 ms, which is the maximum interrogation time of our cold gravimeter), Source 1 behaves better, whereas for large interrogation times, where the dominant contribution to the noise comes from the lower decades (0.1-10 Hz), Source 2 and 3 are better. We assumed here that for any of these sources, the phase noise below 1 Hz is well described as flicker noise, for which the spectral density scales as S φ (f ) = S φ (1Hz)/f 3 . If the phase noise would behave as pure flicker noise over the whole frequency spectrum, the Allan standard deviation of the interferometer phase would scale as T . This is roughly the behavior we notice in the table for long interrogation times.

The sensitivity of the interferometer for acceleration scales as T 2 , so that the sensitivity to acceleration gets better when the interrogation time gets larger. For example, for 2T = 10s and T c = 15 s (resp. 2T = 2 s and T c = 3 s), the phase noise of source 3 would limit the sensitivity to acceleration of the interferometer to 1.8 × 10 -10 m.s -2 (resp. 4.7 × 10 -10 m.s -2 ) at 1 s for 87 Rb.

3

. 4. The frequency chain. -The microwave signal is generated by multiplication of the 100 MHz source. An example of the generation of the microwave reference can be found in [START_REF] Nyman | [END_REF]. The contribution to the interferometer noise of this system was found to be 0.6 mrad per shot for τ R = 10µ s, 2T = 10 s and T c = 15 s, which is negligible with respect to the contribution from the 100 MHz source.

3

. 5. Propagation in the fiber . -In most of the experiment, the Raman lasers are generated by two independent laser sources. The two beams are mixed using a polarizing beam splitter cube, so that they have orthogonal polarizations. A small fraction of the total power can be sent to one of the two exit ports of the cube, where a fast photodetector detects the beat frequency, in order to phase lock the lasers. Both beams are finally guided towards the atoms with a polarization maintaining fiber. Since the Raman beams have orthogonal polarization, fiber length fluctuations induce phase fluctuations, due to the birefringence of the medium. The phase noise induced by the propagation in the fiber can be measured by comparing the beat signal measured after the fiber with the one we use for the phase lock.

Figure 2 displays the power spectral density of the phase noise induced by the propagation, which is dominated by low frequency noise due to acoustic noise and thermal fluctuations. This source of noise can be reduced by shielding the fiber from the air flow of the air conditioning, surrounding it with some packaging foam. An alternative technique consists in using identical linear polarizations for the Raman beams. This can be achieved using a polarizer after the mixing, by generating the Raman lasers by phase modulatuion of a single laser, or by injecting with two independent lasers a power amplifier. In this case, the noise is efficiently suppressed, as shown in figure 2, down to a level where it is negligible.

3

. 6. Detection noise. -Quantum projection noise limit constitutes the intrinsic limit of sensitivity and scales as 1/ √ N , where N is the number of detected atoms per shot. Typical values are below 1 mrad for MOT sources and 3 mrad for ultra cold atomic sources. Thus for long interrogation times, as for space applications, the dominant source of phase noise is expected to be due to the stability of the 100 MHz source.

-The case of parasitic vibrations

The same formalism can be used to evaluate the degradation of the sensitivity to inertial forces caused by vibrations of the retroreflecting mirror.

The sensitivity of the interferometer is then given by where S z (ω) is the power spectral density of position noise. Introducing the power spectral density of acceleration noise S a (ω), the previous equation can be written

σ 2 Φ (τ ) = k 2 ef f τ ∞ n=1 |H(2πnf c )| 2 S z (2πnf c ) (10)
σ 2 Φ (τ ) = k 2 ef f τ ∞ n=1 |H(2πnf c )| 2 (2πnf c ) 4 S a (2πnf c ) (11)
It is important to note here that the acceleration noise is efficiently filtered by the transfer function for acceleration, which decreases as 1/f 4 .

In the case of white acceleration noise S a , and to first order in τ R /T , the limit on the sensitivity of the interferometer is given by :

σ 2 Φ (τ ) = k 2 ef f T 4 2 2T c 3T -1 S a τ (12)
To put this into numbers, we now calculate the requirements on the acceleration noise of the retroreflecting mirror in order to reach a sensitivity of ×10 -10 m.s -2 at 1 s. For negligible dead time (T c ≃ 2T ), the amplitude noise should lie below 2.5 × 10 -10 m.s -2 .Hz -1/2 .

-Systematic effects

Regarding the targeted ultra-high accuracies for atom inertial sensors, many environmental parameters have to be controlled with great precision. An important tool, that allows us to suppress many of the remaining systematic phase shifts, relies on the fact that we can distinguish two classes of systematic effects.

Phase shifts that are independent of the direction of --→ k ef f and others that change sign when reversing --→ k ef f . This allows us to separate and distinguish their influence on the atomic signal in constantly repeated k ↑ -k ↓ measurements. Assuming that the trajectories of the atoms remain constant, the half difference 1 2 △ ↑↓ of the measured interferometer phase ∆Φ int of a consecutive k ↑ -k ↓ measurement will contain only phase terms depending on --→ k ef f . Whereas --→ k ef f independent phase terms are contained in the half sum 1 2 ↑↓ :

1 2 △ ↑↓ = -cos θ| -→ k ef f | • | -→ g | • T 2 + ∆Φ LS2 + ∆Φ Coriolis + ∆Φ Aberr + ... (13) 1 2
↑↓ = ∆Φ RF + ∆Φ LS1 + ∆Φ gradB + ... (14)

5

.

1. k-independent Phase Shifts. -Three major contributions to large phase shifts in the interferometer, that can be rejected from the k-dependent acceleration signal, ---→ k ef f • -→ g T 2 , arise from the presence of magnetic field gradients, 1-photon light shifts [39] and RF-phase shifts (see equ. ( 14)).

We have experimentally demonstrated this rejection on our cold atom gravimeter. To test the rejection of k-independent phase shifts from the actual atomic signal, we have performed differential measurements with 4 configurations (two sets of parameters and k ↑ -k ↓ each). The difference in their half-sums ∆( 12 ) gives us the additionally introduced phase shift. The difference in their half-differences ∆( 12 △) shows the remaining phase error, that is not rejected by the k ↑ -k ↓ measurement.

With this method, we demonstrated rejection efficiencies typically better than 99%, limited by the resolution of the measurement (in the case of the RF and 1-photon light shift), and by the imperfect overlap between the atomic trajectories for the k ↑ -k ↓ interferometers (in the case of the magnetic field gradient). For a space interferometer, the influence of the magnetic field gradients are expected to be reduced with respect to the case of ground experiment, as i) atoms travel along shorter distances, and ii) the two separated wave-packets experience the same phase shifts, as they alternately travel along identical trajectories.

5

.

2. k-dependent Phase Shifts. -As shown above, large systematic error contributions can be removed from the atomic signal by systematic k ↑ -k ↓ measurements. The remaining phase shifts in equation ( 13) are due to effects that are inherently sensitive to the direction of k.

The first term in equ. [START_REF]Hyper-Precision Cold Atom Interferometry in Space (HYPER)[END_REF] represents the actual atomic signal due to the acceleration of the atoms. The atomic transition wavelength, which is determining k ef f , is known to better than 1 kHz [40].

Three major sources of phase errors add onto this acceleration induced phase shift. Similar to the 1-photon light shift, the Eigen-energies of the atomic states are modified in higher order terms by 2-photon transitions of the Raman beams themselves (LS2). This 2-photon light shift introduces a phase shift similar to this of the usual 1-photon light shift, but where the shift of the resonance transition is δν AC (LS2) = Ω ef f 2 /δ Doppler , and δ Doppler = k ef f • v is the Doppler-detuning of the stimulated Raman transition from resonance.

Term 3 and 4 in equ.( 13) are phase terms that both depend on the transverse velocities of the atoms and are difficult to clearly distinguish. Any residual offset velocity v T perpendicular to --→ k ef f will lead to phase errors due to Coriolis forces and imperfect plane wave fronts of the Raman beams (aberrations) [41]. As indicated in figure 3(e), a finite v T will lead to a finite area spanned by the interferometer and thus will make it sensitive to rotations.

Contrary to that, phase shifts due to aberrations in the phase fronts of the Raman beams scale with the temperature T atoms of the atomic ensemble, depending on the shape of distortion in the wave fronts. For a in first order parabolic curvature of the wave front

∆Φ = c • r 2 , we obtain ∆Φ Aberr = c (v 2 T + 2 k B T atoms /m)T 2 .
Here, the interferometer is only sensitive to perturbations introduced in the retro-reflected beam path, which is noncommon to both Raman laser beams. In case of a pure curvature of the retro-reflected phase fronts, ∆a = 10 -10 m/s 2 would require c = 17 rad/m 2 , which corresponds to a phase front radius of R = 240 km, or a flatness of ≈ λ/1500 over a beam diameter of 10 mm. This flatness is difficult to reach, especially with retardation plates.

Irregularly shaped phase distortions, as introduced by mirror and λ/4 plate, can average out and the requirement for the surface flatness is less stringent. To deduce the effect on the atomic signal, the wave front distortion introduced by the individual optical elements have to be measured with a Shack-Hartmann sensor or ZYGO interferometer, at the level of ≈ λ/1000 or better. The collected phase shifts along the classical trajectories of the atoms can then be calculated. Measurements of the phase shift of the interferometer versus initial position and velocity, and versus temperatures can be useful to study this systematic effect, and correct for it.

In space, Coriolis acceleration could also be an issue, depending on the details of the mission. One should though keep in mind that this bias is not intrinsic to the device, but is a part of the signal and is related to the actual trajectory of the satellite. Wavefront aberrations appear as a very important systematic effect, which depends on a non trivial way with the experimental parameters (size and temperature of the atomic cloud and interaction time). To overcome this problem, extremely high flatness optics are required.

-Conclusion

Thanks to a careful study and optimization of all sources of noise and systematics on our gravimeter, we can evaluate the sensitivity and the accuracy for a space accelerometer based on classical cold atom sources and stimulated Raman transition to manipulate the atomic wavepackets. As proposed in [START_REF] Wolf | Submitted to publication to Experimental Astronomy[END_REF], such a accelerometer might be developed easily taking advantage of the similarities with the atomic clock prototype PHARAO for which key components have already been realized and tested on ground for the Engineer Model [35]. Short term sensitivity should not be better than state of the art classic proof mass accelerometers [44], but an atom interferometer should reach a much better long term stability and accuracy, without the need of spinning the satellite. Moreover, the main sources of noise and systematics cancel in a differential measurement, as in the cases of gradiometers [45] or gyroscopes.

Finally, the use of ultra-cold atoms (like BEC) will take full advantage of the space environment, by allowing to increase the measurement time and reducing the system-atics (the effect of wavefront distortion for example) [START_REF]Proposal to Cosmic vision[END_REF]. Such experiments are more complicated and still need to demonstrate their possibilities on ground. This is why several projects are currently developed to improve the knowledge and technology for zero-g environment: the QUANTUS project [46] carry on in Germany, which study the realization of BEC in the free falling Bremen tower, the Frech project ICE [START_REF] Nyman | [END_REF], which tests the realization of an atom interferometer in the zero-g airplane of CNES, or the ESA project "Atom Interferometry Sensors for Space Applications". All these developments will soon give the possibilitiy of using such a device in fundamental physics experiments or in other applications in the field of Earth observation, navigation or exploration of the solar system. * * *

3. -Sensitivity of the interferometer 3 . 1 .

 31 Sensitivity function.

Fig. 1 .

 1 Fig. 1. -Specifications for the phase noise spectral density of different quartz oscillators, transposed at 100 MHz. The phase noise of the source developed for the Pharao project is also displayed as a solid black line.

Fig. 2 .

 2 Fig.2. -Contribution of the propagation in the optical fiber to the Raman phase noise. Grey curve displays the noise in the case of a polarization maintaining fiber and two orthogonal polarizations for the two Raman lasers. Black curve displays the noise in the case of a polarizing fiber and two parallel polarizations for the two Raman lasers.

Fig. 3 .

 3 Fig. 3. -Some major contributions to systematic phase shifts in our cold atom gravimeter. (a)-(c) depend on the direction of k ef f and can be rejected to a residual r, (d)-(e) are indistinguishable from the acceleration signal and need to be controlled within a high precision.

  .

				Source 1	Source 2	Source 3
	Tc	2T	τR	σΦ(Tc)	σΦ(Tc)	σΦ(Tc)
	(s)	(s)	(µs)	(mrad)	(mrad)	(mrad)
	0.25	0.1	10	1.2	3.5	2.2
	3	2	10	21	6.5	4.4
	15	10	10	110	37	19

We would like to thank all our collaborators of the "Inertial Sensors Team" of SYRTE laboratory for their contributions to the experimental results. We also supports the Institut Francilien pour la Recherche sur les Atomes Froids (IFRAF), the European Union (FINAQS), the Délégation Générale pour l'Arrmement and the Centre Nationale d'Etudes Spatiales for their financial supports.