domain: math sdv spi

text:
Maureen Clerc

Juliette Leblond

Jean-Paul Marmorat

Christos Papageorgakis Dedicated

Uniqueness result for an inverse conductivity recovery problem with application to EEG

Keywords: elliptic and Laplace-Poisson PDE, inverse conductivity recovery problem, spherical harmonics, EEG. MS Classification 2010: 31B20, 33C55, 35J05, 35J25, 35Q61, 65R32, 92C55

Considering a geometry made of three concentric spherical nested layers, (brain, skull, scalp) each with constant homogeneous conductivity, we establish a uniqueness result in inverse conductivity estimation, from partial boundary data in presence of a known source term. We make use of spherical harmonics and linear algebra computations, that also provide us with stability results and a robust reconstruction algorithm. As an application to electroencephalography (EEG), in a spherical 3-layer head model (brain, skull, scalp), we numerically estimate the skull conductivity from available data (electrical potential at electrodes locations on the scalp, vanishing current flux) and given pointwise dipolar sources in the brain.

Introduction

We study an inverse conductivity recovery problem in the particular case of a spherical 3D domain Ω (a ball in R 3 ) and for piecewise constant conductivity functions, of which one value is unknown. More precisely, we assume Ω to be made of 3 nested spherical layers, whose conductivity values are known in the innermost and outermost layers. We assume that the elliptic partial differential conductivity equation (conductivity PDE) holds with a given source term in divergence form supported in the innermost layer.

Provided a single measurement as a pair of Cauchy data on the boundary (open subset of the sphere ∂Ω), we will establish uniqueness and stability properties together with a reconstruction algorithm for the intermediate conductivity. We will also perform some analysis in order to investigate robustness of the reconstruction with respect to available measurements and sources information.

We face a very specific version of the many inverse conductivity issues for second order elliptic PDE under study nowadays. This one is related to piecewise constant conductivities in a spherical geometry in R 3 , and set from a single (Cauchy pair of partial) boundary measurement. Similar inverse conductivity recovery problems may be formulated in more general (Lipschitz smooth) domains of arbitrary dimension, with more general conductivities. They are often considered from (several or) infinitely many boundary measurements (pairs of Cauchy data, related Dirichlet-to-Neumann operator), and are called after Calderón, or after medical imaging processes (Electrical Impedance Tomography). Uniqueness and stability conductivity recovery issues are deeply discussed in [START_REF] Alessandrini | Singular solutions of elliptic equations and the determination of conductivity by boundary measurements[END_REF][START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF][START_REF] Alessandrini | Generic uniqueness and size estimates in the inverse conductivity problem with one measurement[END_REF][START_REF] Alessandrini | Lipschitz stability for the inverse conductivity problem[END_REF][START_REF] Isakov | On uniqueness of recovery of a discontinuous conductivity coefficient[END_REF][START_REF] Uhlmann | Electrical impedance tomography and Calderón's problem[END_REF][START_REF] Kenig | The Calderón problem with partial data[END_REF].

More general inverse problems for elliptic PDEs, in particular transmission issues, are discussed in [START_REF] Isakov | Inverse problems for partial differential equations[END_REF][START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF]. Stability properties of Cauchy boundary value problems are described in [START_REF] Alessandrini | The stability for the Cauchy problem for elliptic equations[END_REF] (see also references therein).

A fundamental problem in experimental neuroscience is the inverse problem of source localization, which aims at locating the sources of the electric activity of the functioning human brain using non-invasive measurements, such as electroencephalography (EEG), see [START_REF] Baratchart | Recovery of pointwise sources or small inclusions in 2D domains and rational approximation[END_REF][START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF][START_REF] Badia | An inverse source problem in potential analysis[END_REF][START_REF] Hämäläinen | Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain[END_REF][START_REF] Hämäläinen | Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data[END_REF][START_REF] Kandaswamy | Analytic sensing for multilayer spherical models with application to EEG source imaging[END_REF].

EEG measures the effect of the electric activity of active brain regions through values of the electric potential obtained by a set of electrodes placed at the surface of the scalp [START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF] and serves for clinical (location of epilepsy foci) and cognitive studies of the living human brain.

The inverse source localization problem in EEG is influenced by the electric conductivities of the several head tissues and especially by the conductivity of the skull [START_REF] Vallaghé | A global sensitivity analysis of three-and fourlayer EEG conductivity models[END_REF]. The human skull is a bony tissue consisting of compact and spongy bone compartments, whose distribution and density varies across individuals, and according to age, since humidity of tissues, and therefore their conductivity tends to decrease [START_REF] Nunez | Electric fields of the brain: the neurophysics of EEG[END_REF]. Therefore conductivity estimation techniques are required to minimize the uncertainty in source reconstruction due to the skull conductivity.

Typically, an inverse conductivity estimation problem aims at determining an unknown conductivity value inside a domain Ω from measurements acquired on the boundary ∂Ω. In the EEG case, the measurements can be modeled as pointwise values obtained on a portion of the boundary ∂Ω (the upper part of the scalp) but they are also affected by noise and measurement errors. The questions arising are: the uniqueness of the skull conductivity for known sources inside the brain; the stability of this estimation; and a constructive estimation method.

Quite frequently, for piecewise constant conductivities, the sub-domain (supporting the unknown conductivity value) is also to be determined, in some cases more importantly than the constant conductivity value itself (for example for tumor detection, see [START_REF] Ammari | Polarization and moment tensors, with applications to inverse problems and effective medium theory[END_REF]Ch. 3] and references therein, [START_REF] Kang | Layer potential technique for the inverse conductivity problem[END_REF][START_REF] Kang | A note on uniqueness and stability for the inverse conductivity problem with one measurement[END_REF]). But in the case of EEG, the sub-domains containing the various tissues can be considered known, because they can be extracted from magnetic resonance images. And for simplicity, we only consider the inverse skull conductivity estimation problem in a three-layer spherical head geometry, using partial boundary EEG data. The dipolar sources positions and moments will be considered to be known. This may appear to be an unrealistic assumption because sources reconstruction is itself a difficult inverse problem. But in fact, in some situations there are prior assumptions as to the positions of the sources (in primary evoked electrical potentials), and the position of a source also constrains its orientation, because to the laminar organization of pyramidal neurons in the grey matter.

The overview of this work is as follows. In Section 2, we precise the model and the considered inverse conductivity recovery issue. Our main uniqueness and stability results are stated and proved in Section 3, while an application to EEG and a numerical study are given in Section 4. We then provide a short conclusion in Section 5.

Model, problems

Domain geometry, conductivity

We consider the inverse conductivity estimation problem in a spherical domain Ω ⊂ R 3 made of 3 concentric spherical layers (centered at 0), a ball Ω 0 , and 2 consecutive surrounding spherical shells Ω 1 , Ω 2 . Their respective boundaries are the spheres denoted as S 0 , S 1 , and S 2 , with S i of radius r i such that 0 < r 0 < r 1 < r 2 . We also put Ω

3 = R 3 \ Ω = R 3 \ (Ω ∪ S 2 ).
For i = 0, 1, 2, we assume that σ is a real valued piecewise constant conductivity coefficient with values σ i > 0 in Ω i . Let also σ 3 = 0.

Note that in the present work, the values σ i of the conductivity in Ω i for i = 1 outermost layers Ω 0 , Ω 2 are assumed to be known.

In the EEG framework and for spherical three-layer head models, the domains Ω i respectively represent the brain, the skull and the scalp tissues for i = 0, 1, 2, as shown in Figure 1, see [START_REF] Hämäläinen | Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain[END_REF][START_REF] Hämäläinen | Realistic conductivity geometry model of the human head for interpretation of neuromagnetic data[END_REF]. There, under isotropic assumption, it holds that 0 < σ 1 < σ 0 σ 2 .

Throughout the present work, the geometry Ω and the conductivity σ will be assumed to satisfy the above assumptions.

More general situations are briefly discussed in Remark 3.2 and Section 5.

Figure 1: Spherical head model, with one source C q , p q .

PDE, source terms, statement of the problem

We consider conductivity Poisson equations

∇ • σ∇u = S or div (σ grad u) = S in R 3 , (1)
(in the distributional sense), with a source term S taken to be a distribution on R 3 compactly supported in Ω 0 . We investigate situations where source terms S are of divergence form:

S = ∇ • J P = div J P ,
for distributions J P made of Q pointwise dipolar sources located at C q ∈ Ω 0 with (non zero) moments p q ∈ R 3 :

J P = Q q=1 p q δ Cq , whence S = Q q=1 p q • ∇δ Cq , (2)
where δ Cq is the Dirac distribution supported at C q ∈ Ω 0 . Therefore, in R 3 ,

∇ • σ∇u = Q q=1 p q • ∇δ Cq . (3)
For the EEG case, under the quasi-static approximation and modeling the primary cerebral current J P as in (2), Maxwell's equations imply that the conductivity PDE (3) drives the behaviour of the electric potential u [START_REF] Hämäläinen | Magnetoencephalography theory, instrumentation, and applications to noninvasive studies of the working human brain[END_REF].

In this work, we consider the following inverse conductivity estimation problem in the 3-layered spherical framework of Section 2.1.

From (a single pair of) Cauchy boundary data u = g in a (non-empty) open subset Γ of ∂Ω = S 2 and ∂ n u = 0 on S 2 of a solution to (3), and from a (known) source term S given by (2), we want to recover the constant value σ 1 of the conductivity σ in the intermediate layer Ω 1 .

In Section 3, we establish uniqueness properties of σ 1 from Cauchy data u on Γ ⊂ S 2 , ∂ n u on S 2 and from the source term S. A stability result is also given for Γ = S 2 and equation (1) for more general source terms S will be discussed as well.

Before, we still need to describe the PDE and associated boundary value problems in each of the consecutive layers Ω i .

Laplace-Poisson PDE and transmission issues

For i = 0, 1, 2, 3, write u |Ωi = u i for the restriction to Ω i of the solution u to (3). We put ∂ n u i for the normal derivative of u i on spheres in Ω i , the unit normal vector being taken towards the exterior direction (pointing to Ω i+1 ). In the present spherical setting, we actually have ∂ n = ∂ r .

For i = 1, 2, 3, the following transmission conditions hold on S i-1 , in particular in L 2 (S i-1 ), see [START_REF] Baratchart | Recovery of pointwise sources or small inclusions in 2D domains and rational approximation[END_REF][START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF][START_REF] Badia | An inverse source problem in potential analysis[END_REF] (and Section 2.4):

u i-1 = u i , σ i-1 ∂ n u i-1 = σ i ∂ n u i .
Linked by those boundary conditions, the solutions u i to [START_REF] Alessandrini | Generic uniqueness and size estimates in the inverse conductivity problem with one measurement[END_REF] in Ω i satisfy the following Laplace and Laplace-Poisson equations:

     ∆ u i = 0 in Ω i , i > 0 , ∆ u 0 = 1 σ 0 Q q=1 p q • ∇δ Cq in Ω 0 . (4)
We will see (in Section 3.2.2) that the transmission from

u i ∂ n u i on S i to u i-1 ∂ n u i-1
on S i-1 , for i = 1, 2, may be written

u i-1 ∂ n u i-1 |Si-1 = 1 0 0 σi σi-1 T (S i-1 , S i ) u i ∂ n u i |Si .
for some operator T (S i-1 , S i ) that accounts for the harmonicity of u i in Ω i and that we will express using spherical harmonics.

Provided (a single pair of non identically vanishing smooth enough) Cauchy boundary data u, ∂ n u on Γ ⊂ S 2 , whenever Γ = ∅ is open, and for given conductivity values σ 0 , σ 1 , σ 2 , uniqueness of u holds on S 2 , then on S 1 and S 0 , as ensured by the above formulation and Holmgren's theorem. It is enough to assume that u ∈ W 1,2 (Γ) (the Sobolev-Hilbert space of L 2 (Γ) functions with first derivative in L 2 (Γ)) and ∂ n u ∈ L 2 (Γ), see [START_REF] Brezis | Analyse fonctionnelle[END_REF][START_REF] Chen | Boundary elements methods[END_REF][START_REF] Dautray | Analyse mathématique et calcul numérique[END_REF].

We then face a preliminary data transmission issue from Γ to S 0 , a Cauchy boundary value problem for Laplace equation, which needs to be regularized in order to be well-posed [START_REF] Isakov | Inverse problems for partial differential equations[END_REF]. This is usually done by Thykonov regularization or the addition of an appropriate constraint and may be solved using boundary elements methods, see [START_REF] Atfeh | Bounded extremal and Cauchy-Laplace problems on the sphere and shell[END_REF][START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF] and references therein.

In EEG, data are provided as pointwise values of g at points in Γ (electrodes measurements), and yet another extension step is needed in order to compute an estimate of g on S 2 , using best constrained approximation, see Section 4.

Concerning the source term S, note that it's knowledge only determines u 0 on S 0 up to the addition of a harmonic function in Ω 0 . Indeed, by convolution with a fundamental solution of Laplace equation in R 3 , we see that

u s (x) = 1 4π Q q=1 < p q , x -C q > |x -C q | 3 , x ∈ {C q } , (5)
satisfies u s (x) → 0 at |x| → ∞, ∆ u s = Q q=1 p q • ∇δ Cq , in R 3
, whence in Ω and Ω 0 , and ∆ u s = 0 outside Ω 0 . Solutions u 0 to (4) in Ω 0 are then provided by u s /σ 0 up to the addition of a harmonic function in Ω 0 . The later is in fact (uniquely) determined by the (transmitted) boundary conditions, see [10, Sec. 1.2], [START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF][START_REF] Badia | An inverse source problem in potential analysis[END_REF][START_REF] Kandaswamy | Analytic sensing for multilayer spherical models with application to EEG source imaging[END_REF] where inverse source problems in the EEG setting are discussed, together with reconstruction algorithms.

Associated forward Neumann problem

Let φ ∈ L 2 (S 2 ) (actually it is enough to take φ ∈ W -1/2,2 (S 2 )) of vanishing mean value on S 2 . Then, there exists a solution u to (3) in Ω, Hölder continuous in Ω \ {C q }, which satisfies ∂ n u = φ on S 2 ; it is unique up to an additive constant. In particular, the associated Dirichlet boundary trace u |S2 is Hölder continuous on S 2 . Indeed, looking to u -u s as a (weak) solution to a strictly elliptic PDE in a bounded smooth domain Ω or to a sequence of Laplace equations in the domains Ω i , variational formulation and Lax-Milgram theorem imply that u -u s ∈ W 1,2 (Ω) and the uniqueness property, see [START_REF] Brezis | Analyse fonctionnelle[END_REF][START_REF] Chen | Boundary elements methods[END_REF][START_REF] Dautray | Analyse mathématique et calcul numérique[END_REF].

Hence u-u s belongs to W 1/2,2 (S 2 ) and actually to W 1,2 (S 2 ). That u possesses yet more regularity properties is established in [10, Prop. 1], see also [START_REF] Alessandrini | The stability for the Cauchy problem for elliptic equations[END_REF] for stability results of Cauchy boundary transmission problems.

Conductivity recovery

Uniqueness result

Recall that the geometry Ω and the conductivity coefficients satisfy the hypotheses of Section 2.1. Let Γ ⊂ S 2 a (non empty) open set. Assume the source term S given by ( 2) to be known, and not to be reduced to a single dipolar pointwise source located at the origin (S = p • ∇δ 0 ). Theorem 3.1. Let σ, σ be piecewise constant conductivities in Ω associated to two values σ 1 , σ 1 in Ω 1 and equal values σ 0 , σ 2 in Ω 0 , Ω 2 . If two solutions u, u to (3) associated with σ, σ and such that

∂ n u = ∂ n u = 0 on S 2 coincide on Γ: u |Γ = u |Γ , then σ 1 = σ 1 .
This implies that a single pair of partial boundary Dirichlet data u |Γ on Γ and Neumann data ∂ n u = 0 (vanishing) on S 2 of a solution u to (3) uniquely determines σ 1 > 0.

As the proof in Section 3.3 will show, source terms S that guarantee uniqueness are such that associated Dirichlet data u |Γ on Γ do not identically vanish. Notice also that if no source is present, uniqueness fails (boundary data identically vanish on S 2 ). However, Theorem 3.1 would also hold true for non identically vanishing Neumann on S 2 . We will discuss more general statements the Theorem in Remark 3.2 after the proof, see also Section 5.

In order to establish the result, we use spherical harmonics expansions that we now precise.

Spherical harmonics expansions

In order to express harmonic functions in the spherical shells and balls Ω i and their boundary values on S i , we use the spherical harmonics basis r k Y km (θ, ϕ), r -(k+1) Y km (θ, ϕ), k ≥ 0, |m| ≤ k, in the spherical coordinates (r, θ, ϕ). These are homogeneous harmonic and anti-harmonic polynomials for which we refer to [9, Ch. 9, 10],[15, Ch. II, Sec. 7.3] as for their properties. (the basis functions Y km (θ, ϕ) are products beteween associated Legendre functions of indices k ≥ 0, |m| ≤ k, applied to cos θ and elements of the Fourier basis of index m on circles in ϕ (real or complex valued, cos mϕ, sin mϕ or e ±imϕ ).

Source term, boundary data

The decomposition theorem [9, Thm 9.6], [15, Ch. II, Sec. 7.3, Prop. 6], is to the effect that the restriction u i of u to Ω i for i = 1, 2 may be expanded on the spherical harmonics basis as follows, at (r, θ, ϕ) ∈ Ω i :

u i (r, θ, ϕ) = ∞ k=0 k m=-k α ikm r k + β ikm r -(k+1) Y km (θ, ϕ) ∈ Ω i , (6)
where α ikm and β ikm are the spherical harmonic coefficients of the harmonic and anti-harmonic parts of u i , respectively (harmonic inside or outside

∪ j≤i Ω i ).
Similarly, because it is harmonic in a spherical layer surrounding S 0 , the restriction u 0 of u to Ω 0 is given at points (r, θ, ϕ) with r > max q |C q | > 0 by

u 0 (r, θ, ϕ) = ∞ k=0 k m=-k α 0km r k Y km (θ, ϕ) + u s (r, θ, ϕ) ,
where u s given by ( 5) is expanded there as: r -(k+1) Y km (θ, ϕ):

u s (r, θ, ϕ) = k,m β 0km r -(k+1) Y km (θ, ϕ) . (7)
Here, β 0km are the spherical harmonic coefficients of the anti-harmonic (harmonic outside Ω 0 ) function u S . The normal derivative of u i , i = 0, 1, 2, is then given in Ω i (with r > max q |C q | for i = 0) by:

∂ n u i (r, θ, ϕ) = k,m α ikm kr k-1 -β ikm (k + 1)r -(k+2) Y km (θ, ϕ) (8)
On S i , we put (because u i ∈ L 2 (S i ) where the spherical harmonics form an orthogonal basis [9, Thm 5.12]):

u i (r i , θ, ϕ) = ∞ k=0 k m=-k γ ikm Y km (θ, ϕ) , ∂ n u i (r i , θ, ϕ) = ∞ k=0 k m=-k δ ikm Y km (θ, ϕ) ,
with l 2 summable coefficients γ ikm , δ ikm (that may be real or complex valued depending on the choice for Y km ).

In particular, once the boundary data u 2 = g is extended from Γ to S 2 (see [START_REF] Atfeh | Bounded extremal and Cauchy-Laplace problems on the sphere and shell[END_REF][START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF] and the discussion in Section 2.3), we have:

u 2 (r 2 , θ, ϕ) = k,m γ 2km Y km (θ, ϕ) = k,m g km Y km (θ, ϕ) ,
with g km = γ 2km , whereas the corresponding δ 2km = 0 since ∂ n u 2 = 0 on S 2 (because σ 3 = 0).

Preliminary computations

Below, we write for sake of simplicity, for i = 0, 1, 2:

α ik = α ikm , β ik = β ikm , γ ik = γ ikm , δ ik = δ ikm , g k = g km ,
for all k ≥ 0, and every |m| ≤ k (we could also take the sums over |m| ≤ k).

Recall from Section 2.3 that the following transmission conditions hold on S i-1 for i = 1, 2, 3:

Σ i-1 u i-1 ∂ n u i-1 |Si-1 = Σ i u i ∂ n u i |Si-1 , (9)
with

Σ i = 1 0 0 σ i hence Σ -1 i = 1 0 0 1 σi and σ i Σ -1 i = σ i 0 0 1 .
By projection of ( 6), [START_REF] Atfeh | Bounded extremal and Cauchy-Laplace problems on the sphere and shell[END_REF], onto (the orthogonal L 2 (S i ) basis of) spherical harmonics, and with

T k (r i ) = r k i r -(k+1) i kr k-1 i -(k + 1)r -(k+2) i ,
we obtain for all k ≥ 0 the following relations on S i :

γ ik δ ik = T k (r i ) α ik β ik .
In particular:

β ik = r k+1 i 2k + 1 (k γ ik -δ ik ) . (10)
The transmission conditions (9) through S i-1 express as:

Σ i-1 γ i-1k δ i-1k = Σ i T k (r i-1 ) α ik β ik .
Because T k (r i ) is invertible (r i > 0), this implies that:

γ i-1k δ i-1k = Σ -1 i-1 Σ i T k (r i-1 ) T k (r i ) -1 γ ik δ ik .
Therefore, in the spherical geometry,

T (S i-1 , S i ) = T k (r i-1 ) T k (r i ) -1
for the operator T (S i-1 , S i ) introduced at the end of Section 2.3.

CLERC, LEBLOND, MARMORAT AND PAPAGEORGAKIS

Hence, because γ 2k = g k and δ 2k = 0:

δ 0k γ 0k = Σ -1 0 Σ 1 T k (r 0 ) T k (r 1 ) -1 Σ -1 1 Σ 2 T k (r 1 ) T k (r 2 ) -1 g k 0 , (11)
while

β 0k = [0 1] T k (r 0 ) -1 δ 0k γ 0k . ( 12
)
These formula express a linear relation between the source term coefficients β 0k and the boundary Dirichlet data with coefficients g k , which is studied in Appendix and gives rise to (13) below. We already see the particular role of σ 1 that appears through Σ -1 1 and Σ 1 . This explains why, after multiplication by σ 1 and algebraic manipulations, we obtain in (13) a polynomial of degree 2 in σ 1 .

Algebraic equations

As computed in Appendix, equations ( 11), ( 12) can be rewritten, for all k ≥ 0, as:

B 1 (k) σ 1 β 0k = A 2 (k) σ 2 1 + A 1 (k) σ 1 + A 0 (k) g k , (13)
with non negative quantities A i (k), i = 0, 1, 2, B 1 (k) that depend only on the geometry, on the given conductivity values σ 0 , σ 2 , and on k. Actually, A 1 (k), B 1 (k) > 0 for all k ≥ 0 while A 0 (k), A 2 (k) > 0 for k > 0 but A 0 (0) = A 2 (0) = 0. In particular, for all k ≥ 0 and for σ 1 > 0, we have

A 2 (k) σ 2 1 + A 1 (k) σ 1 + A 0 (k) > 0.
This implies that β 0k = 0 ⇔ g k = 0 and that for all k such that g k = 0, β 0k /g k is real valued positive: the spherical harmonics basis diagonalizes the transmission relations.

Uniqueness proof

Proof. (Theorem 3.1) Assume that there exists another value σ 1 > 0 of the conductivity in Ω 1 that gives rise to the same potential (and vanishing current flux) on Γ ⊂ S 2 , from the same source term u s (same boundary measurements and coefficients g k = g km , same sources term coefficients β 0k = β 0km , given). Equation ( 13) then holds for both σ 1 , σ 1 > 0. We thus get that either β 0k = g k = 0 or

β 0k g k = A 2 (k) σ 2 1 + A 1 (k) σ 1 + A 0 (k) B 1 (k) σ 1 = A 2 (k) σ ,2 1 + A 1 (k) σ 1 + A 0 (k) B 1 (k) σ 1 , whence A 2 (k) σ 2 1 + A 1 (k) σ 1 + A 0 (k) B 1 (k) σ 1 - A 2 (k) σ ,2 1 + A 1 (k) σ 1 + A 0 (k) B 1 (k) σ 1 = 0 , hence multiplying by B 1 (k) σ 1 σ 1 > 0: (σ 1 -σ 1 ) (A 2 (k) σ 1 σ 1 -A 0 (k)) = 0 .
Thus either σ 1 = σ 1 and uniqueness holds or, for all values of k ≥ 0 such that

β 0k = 0, A 0 (k) = σ 1 σ 1 A 2 (k) .
This holds for k = 0 but for k > 0 it implies that

A 0 (k) A 2 (k) = σ 1 σ 1 ,
which could not be true for more than a single value of k > 0. Indeed, the product σ 1 σ 1 is constant while A 0 (k)/A 2 (k) stricly increases with k, as we now show. We have:

A 0 (k) A 2 (k) = σ 0 σ 2 k 1 -r1 r2 2k+1 (k + 1) r1 r2 2k+1 + k = σ 0 σ 2 k 1 -2k+1 (k + 1) 2k+1 + k , (14)
with = r 1 /r 2 < 1, and we put:

E(k) = 1 σ 0 σ 2 A 0 (k) A 2 (k) = 1 -2k+1 1 + k+1 k 2k+1 , k > 0 , E(0) = 0 .
Because for k > 0, k + 2/(k + 1) < (k + 1)/k and 2k+3 < 2k+1 , the numerator of E(k) strictly increases with k while its denominator strictly decreases. Thus, E is a strictly increasing function of k, which converges to 1 as k → ∞.

Hence, among the k > 0, the equation E(k) = σ1σ 1 σ0 σ2 admits at most one solution, and pairs σ 1 , σ 1 > 0 cannot solve A 2 (k)σ 1 σ 1 -A 0 (k) = 0 for more than 1 value of k > 0 (actually, a necessary condition for σ 1 , σ 1 to solve A 2 (k)σ 1 σ 1 -A 0 (k) = 0 for 1 value of k > 0 is that σ 1 σ 1 ∈ (0 , σ 0 σ 2 )). So we must have σ 1 = σ 1 , as soon as β 0k (or g k ) does not vanish for at least 2 distinct values of k.

Finally, we show that potentials u s associated to pointwise dipolar source terms S = p • ∇δ 0 have at least 2 non-null coefficients β 0k in their spherical harmonic expansion. Indeed, assume that all the coefficients β 0km are 0, except for a single value of k > 0, say k 0 . From [START_REF] Ammari | Polarization and moment tensors, with applications to inverse problems and effective medium theory[END_REF], the function u s is then a antiharmonic homogeneous polynomial of degree k 0 and for r > max q |C q | > 0,

u s (r, θ, ϕ) = |m|≤k0 β 0k0m r -(k0+1) Y k0m (θ, ϕ) = 1 r 2k0+1 |m|≤k0 β 0k0m r k0 Y k0m (θ, ϕ) .
From [9, Ch. 5], the distribution S = ∆u s also coincides far from {0} with a polynomial divided by an odd power of r. This contradicts the assumptions on S, which has a pointwise support in Ω 0 not reduced to {0}. Remark 3.2. Theorem 3.1 is in fact valid for solutions to Equation (1) with more general source terms S. To ensure uniqueness, it is indeed enough to assume that u s does not coincide with some homogeneous anti-harmonic polynomial of positive degree, so that it admits on S 0 at least two coefficients β 0k = 0.

In the present spherical geometry, note that u s is equal on S 0 to a homogeneous harmonic polynomial if and only if so is u on S 2 .

The last part of the proof actually implies that potentials u s issued from pointwise dipolar source terms S with support in Ω 0 not reduced to {0} have infiniteley many coefficients β 0k = 0.

Stability properties

We now establish a stability result for the inverse conductivity estimation problem with respect to the source term whenever Γ = S 2 . Proposition 3.3. Assume the source terms S, S and the conductivities σ, σ to satisfy the assumptions of Theorem 3.1. Let u s , u s be the associated potentials through [START_REF] Alessandrini | The stability for the Cauchy problem for elliptic equations[END_REF]. Let u, u be the associated solutions to (3) such that ∂ n u = ∂ n u = 0 on S 2 . Put g, g for their boundary values on S 2 . Then, there exist c, c s > 0 such that

|σ 1 -σ 1 | ≤ c g -g L 2 (S2) + c s u s -u s L 2 (S0) . Whenever 0 < s m ≤ σ 1 , σ 1 ≤ s M for constants s m , s M , then c, c s do not depend on σ 1 , σ 1 but on s m , s M .
Remark 3.4. For ordered lists of sources (p q , C q ) , (p q , C q ) with length Q, we can define the geometric distance

d(S, S ) = Q q=1 |p q -p q | + |C q -C q | .
If the sources are located far enough from S 0 in the sense that max(|C q |, |C q |) ≤ ρ < r 0 , and because u s is on S 0 a continuous function of p q , C q , we can rewrite the inequality in Proposition 3.3 as:

|σ 1 -σ 1 | ≤ c g -g L 2 (S2) + c s d(S, S ) ,
with c s = K(ρ) c s for some constant K(ρ) which depends on ρ. Hence, the conductivity σ 1 depends continuously on the (complete) Dirichlet boundary data g (in L 2 (S 2 )) and on the source term S, with appropriate topology.

Notice also the relation:

β 0km = 1 2k + 1 Q q=1 p q , ∇ r k Y km (θ, ϕ) (C q ) L 2 (S0) .
Finally, observe that the constants c, c s , c s in the above inequalities also depend on the data g whence on S . The dependence between Dirichlet data g on S 2 and the source term S can be precised by using, for instance, the last equality together with relation (13) between their coefficients (g k ) and (β 0k ), and then recalling the assumption

s m ≤ σ 1 ≤ s M . Proof. Let ε k (σ 1 , β 0k , g k ) = B 1 (k) σ 1 β 0k -A 2 (k) σ 2 1 + A 1 (k) σ 1 + A 0 (k) g k , (15)
It follows from (13) that

σ 1 ε k (σ 1 , β 0k , g k ) -σ 1 ε k (σ 1 , β 0k , g k ) = 0 , whence we get 0 = σ 1 ε k (σ 1 , β 0k -β 0k , g k -g k ) + σ 1 ε k (σ 1 , β 0k , g k ) -σ 1 ε k (σ 1 , β 0k , g k ) .
But

σ 1 ε k (σ 1 , β 0k , g k ) -σ 1 ε k (σ 1 , β 0k , g k ) = g k (σ 1 -σ 1 ) [A 2 (k)σ 1 σ 1 -A 0 (k)] , so g k (σ 1 -σ 1 ) [A 2 (k)σ 1 σ 1 -A 0 (k)] = -σ 1 B 1 (k)σ 1 (β 0k -β 0k ) -(A 2 (k)σ 2 1 + A 1 (k)σ 1 + A 0 (k))(g k -g k ) .
Recall that A 2 (k) > 0 for k > 0 and other arguments of Section 3.3 together with the computations in Appendix show that

A 0 (k)/A 2 (k), A 1 (k)/A 2 (k), r k+1 0 B 1 (k)/A 2 (
k) are uniformely bounded in k from above and from below (by strictly positive constants). We can then divide by A 2 (k), in order to obtain

|σ 1 -σ 1 |   k,m |g k | 2 σ 1 σ 1 - A 0 (k) A 2 (k) 2   1 2 ≤ √ 2 σ 1   k,m σ 2 1 + A 1 (k) A 2 (k) σ 1 + A 0 (k) A 2 (k) 2 |g k -g k | 2   1 2 + √ 2 σ 1 σ 1   k,m B 2 1 (k) A 2 2 (k) |β 0k -β 0k | 2   1 2
.

In order to establish an upper bound, note that

k,m B 2 1 (k) A 2 2 (k) |β 0k -β 0k | 2 ≤ sup k r 2(k+1) 0 B 2 1 (k) A 2 2 (k) k,m r -2(k+1) 0 |β 0k -β 0k | 2 , with k,m r -2(k+1) 0 |β 0k -β 0k | 2 = u s -u s 2 L 2 (S0) .
Moreover, since g k (as β 0k ) possess non vanishing values for infinitely many (at least two) values of k, and σ 1 σ 1 -A 0 (k)/A 2 (k) can vanish for at most one value of k, it holds that:

k,m |g k | 2 |σ 1 σ 1 -A 0 (k)/A 2 (k)| 2 > 0 ,
can be bounded from below under the assumptions on σ 1 σ 1 .

Application to EEG

In order to illustrate Proposition 3.3, we now perform a short numerical analysis of the inverse conductivity estimation problem in the spherical domain and the EEG setting described in Section 2. Measurements of the Dirichlet data g on the scalp S 2 (pointwise values at electrodes locations) and known sources activity are expanded on the spherical harmonics basis, using the FindSources3D software1 (FS3D), see also [START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF]. We therefore have at our disposal the spherical harmonics coefficients (g km , β 0km ) for 0 ≤ k ≤ K for some K > 0 and |m| ≤ k.

Reconstruction algorithm

As the reconstruction of the conductivity σ 1 does not depend on the spherical harmonics indices m, in order to increase the robustness of our reconstruction algorithm, the following normalization is applied over the different spherical harmonics indices k:

       gk = |m|≤k g km β0km , β0k = |m|≤k β 0km β0km = |m|≤k |β 0km | 2 .
There, β0km is the complex conjugate number to β 0km (indeed, β 0km could be complex valued if the basis elements Y km are taken in their complex valued form).

The procedure is a least square minimization of the error equation obtained from [START_REF] Dautray | Analyse mathématique et calcul numérique[END_REF] as a truncated finite sum for K > 0:

σ est 1 = arg min s K k=0 ε k (s, β0k , gk ) 2 .
(16)

Numerical illustrations

We consider the EEG framework in the spherical three-layer head model, as described in Section 2.1, where the layers represent the brain, the skull and the scalp tissues, respectively. The radii of the spheres used in the numerical analysis are normalized to the values r 0 = 0.87, r 1 = 0.92 and r 2 = 1. In the present analysis, the brain and scalp tissue conductivities are set to σ 0 = σ 2 = 0.33 S/m, while the skull conductivity σ 1 is to be recovered. When generating simulated EEG data through the associated forward simulation, we will set σ 1 = 0.0042 S/m. Our study uses simulated data associated to a single dipole and the minimization of [START_REF] Badia | An inverse source problem in potential analysis[END_REF] for the conductivity estimation. The algorithm is written as a MATLAB code and the forward simulations are run with the FS3D software.

We validate our reconstruction algorithm using simulated EEG data by FS3D (for solving the direct EEG problem). Of course, the EEG data are subject to some ambient noise and measurements errors, and the a priori knowledge on the sources is not perfect. The inverse conductivity estimation problem is sensitive to such perturbations though it possesses the stability property described in Proposition 3.3.

To investigate the stability of our algorithm with respect to the source term, we select a source term S made of a single dipole located at C 1 = (0.019, 0.667, 0.1), mimicking an EEG source at the frontal lobe of the brain, with moment p 1 = (0.027, 0.959, 0.28). The associated spherical harmonics CLERC, LEBLOND, MARMORAT AND PAPAGEORGAKIS coefficients gk and β0k are computed for 0 ≤ k ≤ K = 30. The original source location C 1 is replaced by inexact locations C n 1 for n = 1, • • • , 20 located at a constant distance from C 1 (a percentage of the inner sphere radius r 0 ), as illustrated in Figure 2, while the source moment p 1 is retained. For each new dipole location C n 1 , the associated spherical harmonics coefficients βn 0k are simulated. We perform conductivity estimation from the pairs gk , βn 0k (recall that gk correspond to the actual β0k ). The effect of the source mislocation on the conductivity estimation is summarized in Figure 3 and Table 4.2 which respectively shows and lists the values and other characteristics of the estimated conductivities with respect to the distance between actual and inexact sources.

These preliminary results illustrate the influence of source mislocation on conductivity estimation, and the robustness character of our algorithm, in accordance with the stability result of Proposition 3.3. In order to penalize high frequencies and to get more accurate estimations, we will in particular introduce in the above criterion ( 16) appropriate multiplicative weights (decreasing with the index k). Dipole mislocation (% of radius r 0 ) σest

Conclusion

Observe that our uniqueness result, Theorem 3.1, may be expressed as an identifiability property of the conductivity value (model parameter) σ 1 in the relation (transfer function) from boundary data to sources (control to observation), [START_REF] Chavent | The output least squares identifiability of the diffusion coefficient from an H 1 observation in a 2D elliptic equation[END_REF][START_REF] Leblond | Identifiability properties for inverse problems in EEG data processing and medical engineering, with observability and optimization issues[END_REF]. This could be useful in order to couple EEG with additional modalities, like EIT (where ∂ n u = 0 is known on Γ) or even MEG (magnetoencephalography, which measures the magnetic field outside the head), and to simultaneously estimate both σ 1 and the source term S in situations where the latter is (partially) unknown. Following Remark 3.2, we may also wish to recover possibly unknown information about the (spherical) geometry of Ω 1 (like r 1 or/and r 0 ).

Situations with more than 3 spherical layers could be described similarly, which may help to consider more general conductivities (smooth but non constant) by piecewise constant discretization.

As Theorem 3.1, Proposition 3.3 would still hold true under a weaker sufficient condition for the source terms, according to which the associated potential on S 0 through (5) should admit at least 2 non-null coefficients (this is equivalent to the same property for g on S 2 , see Remark 3.2). Moreover, it could be extended to a stability property with respect to boundary data with close and non vanishing Neumann data on S 2 . However, stability properties for situations with partial Dirichlet boundary data only (on Γ ⊂ S 2 ) would be weaker, see e.g. [START_REF] Atfeh | Bounded extremal and Cauchy-Laplace problems on the sphere and shell[END_REF][START_REF] Alessandrini | The stability for the Cauchy problem for elliptic equations[END_REF][START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF].

We have also begun to study the same uniqueness and stability issues in more general (non-spherical) nested geometries, see [START_REF] Ammari | Polarization and moment tensors, with applications to inverse problems and effective medium theory[END_REF][START_REF] Kang | Layer potential technique for the inverse conductivity problem[END_REF][START_REF] Kang | A note on uniqueness and stability for the inverse conductivity problem with one measurement[END_REF][START_REF] Nédélec | Acoustic and electromagnetic equations integral representations for harmonic problems[END_REF], and also [START_REF] Alessandrini | Open issues of stability for the inverse conductivity problem[END_REF] for a number of open problems.

Acknowledgement

This work was partly supported by the Région Provence-Alpes-Côte d'Azur (France) and the company BESA GmbH (Germany). The authors thank the reviewer for providing constructive suggestions.

Appendix: More computations related to Section 3.2.2 From (11), [START_REF] Chavent | The output least squares identifiability of the diffusion coefficient from an H 1 observation in a 2D elliptic equation[END_REF], we get for all k ≥ 0,

β 0k = g k × [0 1] T k (r 0 ) -1 Σ -1 0 Σ 1 T k (r 0 ) T k (r 1 ) -1 Σ -1 1 Σ 2 T k (r 1 ) T k (r 2 ) -1 1 0 .
The matrices T k (r i ) and T k (r i ) -1 can be written:

T k (r i ) = 1 0 0 1 ri 1 1 k -(k + 1)
r k i 0 0 r -(k+1) i , T k (r j ) -1 = 1 2k + 1 r -k j 0 0 r (k+1) j k + 1 1 k -1 1 0 0 r j .
Their products that give an expression of T (S i-1 , S i ) in the spherical geometry are then such that:

T k (r i-1 ) T k (r i ) -1 = 1 2k + 1 × 1 0 0 r -1 i-1 1 1 k -(k + 1)    ri-1 ri k 0 0 ri ri-1 k+1    k + 1 1 k -1 1 0 0 r i = 1 2k + 1 r i r i-1 k+1 × 1 0 0 r -1 i-1 1 1 k -(k + 1)
ri-1 ri

2k+1 0 0 1 k + 1 1 k -1 1 0 0 r i .
. We can write

T k (r i-1 ) T k (r i ) -1 = ρ (i) k   a (i) k b (i) k c (i) k d (i) k   , with ρ (i) k = 1 2k + 1 r i r i-1 k+1 , i = 1, 2 , ρ (0) k = r k+1 0 2k + 1 ,
and the real valued quantities, with their equivalent asymptotic behaviours as k → ∞:

                               a (i) k = (k + 1) ri-1 ri 2k+1 + k ∼ k , b (i) k = r i ri-1 ri 2k+1 -1 ∼ -r i , c (i)
k = k (k+1) ri-1 ri b (i) k ∼ -k 2 ri-1 , d (i) k = ri ri-1 k ri-1 ri 2k+1 +k + 1 ∼ kri ri-1 .
Define also the real valued quantities e (0)

k , f (0)
k :

e (0) k = k , f (0)
k = f (0) = -r 0 .

We have

[0 1] T k (r 0 ) -1 = ρ (0) k e (0) k f (0) k .

Figure 2 :

 2 Figure 2: Locations (in Ω 0 ) of C 1 (red bullet) and of the 20 points C n 1 (blue cross) surrounding it, for |C 1 -C n 1 | equal to 10% of r 0 .

1 Figure 3 : 1 among n = 1 ,

 1311 Figure 3: Conductivity estimation results for various mislocations of the actual dipole used to simulate the EEG data: 20 dipole locations C n 1 are selected by displacing C 1 by a constant distance, computed as a percent of the brain radius r 0 (on the abcissa axis). Displayed are: σ 1 , the actual conductivity value used in the EEG data simulation, σ est 1 , the estimated conductivity value for each dipole position C n 1 , and σest 1 , the mean value of σ est 1 among n = 1, • • • , 20.

Table 1 :

 1 Conductivity estimation results, continued; Listed in columns are: (i) the distance between C 1 and C n 1 , (ii) the mean estimated conductivity value σest 1 , (iii) the standard deviation of σ est 1 , (iv) the mean value of the relative errors between σ 1 and σ est 1 .

			Standard	Mean of
		1	deviation	relative errors
	0	4.200e-03	0	1.858e-15
	0.1	4.195e-03 1.450e-05	3.123e-03
	1	4.187e-03 1.629e-04	3.318e-02
	5	4.160e-03 7.703e-04	1.511e-01
	10	4.741e-03 1.512e-03	3.350e-01

See http://www-sop.inria.fr/apics/FindSources3D/.

CLERC, LEBLOND, MARMORAT AND PAPAGEORGAKIS

Then, equation [START_REF] Chen | Boundary elements methods[END_REF] holds true with:

k a

(2)

Observe that r k+1 0 B 1 (k) acts on r -(k+1) 0 β 0k that are members of an l 2 sequence (see Section 3.2.1 and equation ( 10) with i = 0).

As in [START_REF] Clerc | Source localization using rational approximation on plane sections[END_REF], one can show with the above expressions that the behaviours of the ratios A i (k)/A 2 (k), B 1 (k)/A 2 (k) ensure that they all are uniformly bounded from below on from above by positive constants, for k > 0 .

Note also that

where Ãi , Ãij , B1 only depend on the spherical geometry.


Ramin Khayatzadeh
email: ramin.khayatzadeh@univ-amu.fr

Mathieu Collombon

Didier Guyomarc'h

Didier Ferrand

Gaëtan Hagel

Marie Houssin
email: marie.houssin@univ-amu.fr.

Olivier Morizot

Caroline Champenois

Martina Knoop

Noise characterization of an Optical Frequency Comb using Offline Cross-Correlation

Keywords: Mode-locked lasers, Phase noise, Digital signal processing, Cross correlation, Difference frequency generation

Using an offline cross-correlation technique, we have analyzed the noise behavior of a new type of optical frequency comb (OFC), which is carrier envelope offset (CEO) free by configuration, due to difference frequency generation. In order to evaluate the instrument's ultimate noise floor, the phase and amplitude noise of a stabilized OFC are measured simultaneously using two analog-to-digital converters. Carrier recovery and phase detection are done by post-processing, eliminating the need for external phase-locked loops and complex calibration techniques. In order to adapt the measurement noise floor and the number of averages used in cross correlation, an adaptive frequency resolution for noise measurement is applied. Phase noise results are in excellent agreement with measurements of the fluctuations of the repetition frequency of the OFC obtained from optical signal.

I. INTRODUCTION

O PTICAL frequency combs (OFCs) find a multitude of applications in different scientific and engineering domains. One of the prime applications of OFCs is optical frequency measurement where an OFC, thanks to its large frequency span, behaves as a ruler to measure frequency differences [START_REF] Udem | Optical frequency metrology[END_REF], [START_REF] Pastor | Absolute frequency measurements of the 2 3 S 1 → 2 3 P 0,1,2 atomic helium transitions around 1083 nm[END_REF]. Frequency combs that emit at telecommunication wavelengths are especially attractive for microwavephotonics and radio-over-fibre communication systems in order to generate millimetre wave and low-phase-noise carriers which can hugely increase communication data rates [START_REF] Stöhr | 60 GHz radio-over-fiber technologies for broadband wireless services[END_REF], [START_REF] Khayatzadeh | Impact of phase noise in 60-GHz radio-over-fiber communication system based on passively modelocked laser[END_REF]. Recently, offset-free OFCs have been introduced in which the carrier-envelope offset (CEO) frequency, f 0 , cancels out in a non-linear process of difference frequency generation (DFG) [START_REF] Nakamura | Offset-free broadband Yb:fiber optical frequency comb for optical clocks[END_REF], [START_REF] Puppe | Characterization of a DFG comb showing quadratic scaling of the phase noise with frequency[END_REF], [START_REF] Krauss | All-passive phase locking of a compact Er: fiber laser system[END_REF]. Ultimate performances on all OFC modes can be reached by locking the OFC to an ultra-stable frequency reference [START_REF] Mcferran | Low-noise synthesis of microwave signals from an optical source[END_REF] by means of a single parameter, the repetition frequency f rep . Because of the practical importance of OFCs for the mentioned applications, the characterization of noise in OFC lasers has received considerable attention. There are many techniques in both RF and optical domains for measuring the phase and amplitude noise of OFCs. Traditionally, phase noise measurements are made by analogue techniques [START_REF]Agilents phase noise measurement solutions[END_REF], [START_REF] Walls | Extending the range and accuracy of phase noise measurements[END_REF], [START_REF] Howe | Properties of signal source and measurement methods[END_REF]. Recently new approaches based on digital measurements have attracted many attentions. For example, in [START_REF] Khayatzadeh | Accurate millimeter-wave laser phase noise measurement technique[END_REF], the authors present an accurate phase-noise measurement technique based on sampling and post-processing which is able to separate amplitude and phase noise measurements. However, when the phase noise of the device under test (DUT) is much lower than the one of the measurement device, the measure is limited by the measurement device noise. This impact can be removed by implementing a cross-correlation technique [START_REF] Rubiola | Correlation-based phase noise measurements[END_REF], [START_REF] Feldhaus | A 1MHz to 50 GHz direct down-conversion phase noise analyzer with cross-correlation[END_REF]. In [START_REF] Grove | Direct-digital phase-noise measurement[END_REF], a real-time direct digital phase noise measurement device based on cross correlation is developed. This set-up uses two analog-to-digital converters (ADC) in order to sample signal from the DUT, and then measures phase noise by using a digital phase detector and applying cross correlation via a FPGA processor.

In this paper, a noise detection and cross-correlation technique is presented to measure the amplitude and phase noise of a stabilized OFC using only two ADCs. The described method has been implemented completely offline. Comparing this technique to the real-time analog techniques, there is no need of phase-locked loops and phase detectors at the inputs since the phase and amplitude detection are performed by post-processing. Other advantages of this technique can be mentioned as follows: the oscillator noise can be compared at different frequencies, amplitude and phase noise can be measured simultaneously and the cost of this cross correlator is largely inferior to commercial devices. The measured results are compared with the direct optical evaluation presented in [START_REF] Puppe | Characterization of a DFG comb showing quadratic scaling of the phase noise with frequency[END_REF] and a very good agreement is found.

II. SET-UP SCHEME AND PRELIMINARY TESTS A. Configuration of the cross-correlation

Figure 1 presents the schematic diagram of the measurement set-up. In this diagram, the output signal of the DUT is divided into two parts and after passing through anti-harmonics and anti-aliasing filters (LPF1 and LPF2), each part is sampled using an ADC. Both ADCs used in this set-up are independent and they work using different clocks. After sampling, a complete post-processing approach is used to detect the phase noise of the DUT using an I/Q phase detector and cross correlation technique. These processes are presented in the dashed box in Fig. 1. In the first step, the quadrature component of the sampled signal is determined using the Hilbert transfer function (HT). Then, knowing the in-phase and quadrature components, the phase of the sampled signal (φ 1 and φ 2 ) can be determined by using unwrapped arctangent function. In the next step, the carrier frequency is removed from the detected phase whose general expression is

φ(t) = 2πf c t + ϕ(t) (1)
where f c is the carrier frequency and ϕ(t) is the detected phase noise.

We need to mention that since the detected phase using the I/Q detector in the previous step is unwrapped, its value contains the angular frequency multiplied by time values (see Eq. 1). Thus, considering that the carrier frequency is much higher than the frequency fluctuations, one can find the exact value of the carrier frequency utilizing a linear polynomial curve fitting technique.

After determining f c and removing the term 2πf c t, one can extract the phase noise value, ϕ(t), from Eq.( 1). The cross correlation technique is applied in the next step to reduce the noise floor of the measurement system.

Additionally, the amplitude noise (t) can be measured using an envelope detector based on in-phase and quadrature components as follow:

(t) = I 2 + Q 2 -A A (2)
where I and Q are the in-phase and quadrature components and A is the mean amplitude of the signal.

B. Simulation and calibration

Prior to measurements on the OFC, a simulation is performed in order to verify the accuracy of phase and amplitude detectors used in the cross correlation technique. To check the phase noise detector, a sine function signal is corrupted with phase and amplitude noise and by using the phase and amplitude detectors the power spectral density (PSD) of the detected noise and the injected noise are compared. The phase noise is generated by integrating a white noise signal giving a slope of the PSD of phase noise at -20 dB per decade which presents the impact of 1/f noise. The amplitude noise is considered to be white Gaussian noise. Figure 2 illustrates the PSD of the injected phase noise (black bold curve) and the detected phase noise (gray dashed curve) and the PSD of the injected (black dash curve) and detected (dashed dot curve) amplitude noise. As can be seen, the noise detected by the described measurement set-up and the injected noise are identical, which confirms the accuracy of the phase and amplitude detectors. This simulation permits us to calibrate the phase and amplitude noise detectors.

III. EXPERIMENTAL SET-UP

The block diagram of the experimental set-up of the measurement system is presented in Fig. 3. Passive mode-locking [START_REF] Mollenauer | The soliton laser[END_REF] of an Erbium fiber laser at 1550 nm generates femtosecond pulses which are the basis of the OFC. After spectral broadening in a photonic crystal fiber [START_REF]Difference frequency comb[END_REF] and DFG [START_REF] Nakamura | Offset-free broadband Yb:fiber optical frequency comb for optical clocks[END_REF], [START_REF] Puppe | Characterization of a DFG comb showing quadratic scaling of the phase noise with frequency[END_REF], [START_REF] Krauss | All-passive phase locking of a compact Er: fiber laser system[END_REF], which both preserve the phase relationship between all modes, a super-continuum light spanning from 1450 nm to 1625 nm is obtained. It is composed by equally spaced modes separated by the repetition frequency f rep without offset frequency, mode N has a frequency f N = N ×f rep . Other phasecoherent processes, allow to access the wavelengths used in our set-up designed for precision spectroscopy of trapped Ca + -ions [START_REF] Champenois | Terahertz frequency standard based on three-photon coherent population trapping[END_REF]. For long-term stability, the tenth harmonics of the repetition frequency 10f rep at 800 MHz, filtered by a band-pass filter (BPF1), is phase-locked to a low-noise crystal oscillator referenced to GPS (Global positioning system). Corrections of the servo system are applied to the piezoelectric element of the OFC laser after applying proportional control gain (K p ).

To investigate f rep noise in a cross-correlation method, one needs to detect f rep with two fully uncorrelated experimental benches. We use two photodiodes at 794 nm and 866 nm, to detect the OFC signal via 3-meter long singlemode, polarization maintaining fibers which are shorter than the coherence length of the signal. Each mode beats with all its neighbors, and therefore each photodiode signal is composed of a fundamental beat at 80 MHz and harmonics. After bandpass filtering (BPF2 and BPF3 in Fig. 3) a beat note is obtained at f rep =80 MHz exhibiting a 40 dB signal to noise ratio. In the next step, the two signals are frequency down-converted to lower frequencies (250 Hz) in order to be sampled by two analog to digital converters (ADC). The local oscillators (LO1 and LO2) are two independent signal generators (HP 8657A and Marconi 2022C) and the ADCs are two National Instrument acquisition cards (16 bit, 2 Msps) driven by two different clock signals. After signal acquisition, a complete off-line post-processing cross-correlation analysis is applied to eliminate all non-common noises (PD, LO, down-conversion, ADC noises) which reduces the noise floor of the measurement and allows to determine phase noise of the repetition frequency f rep .

In order to verify the accuracy of the measurement results, the PSD of the phase noise of the repetition frequency is deduced from a phase noise measurement directly in the optical domain. The block diagram of the set-up is depicted in Fig. 4 and the results are compared with those of the crosscorrelation technique (See Fig. 5). The optical measurement is performed utilizing an ultra-stable frequency Titaniumsapphire (Ti:Sa) laser as a reference. The frequency of the laser is stabilized using an ultra low expansion (ULE) cavity and applying Pound-Drever-Hall technique. This laser has a linewidth of 2 Hz (measured at resolution frequency of 1 Hz) and a frequency stability of better than 2.10 -14 at 1 second. In the optical measurement technique, we record the beat note signal between the reference laser signal at 729 nm and the nearest optical mode of the CEO free OFC at f 729nm (f 729nm = N 729nm × f rep , N 729nm = 5137500). This beat-note is filtered using BPF4 (around 48 MHz) and then sampled using a digital storage oscilloscope (DSO) in order to measure the PSD of its phase noise. Since its phase noise is N 729nm times higher than that of the repetition frequency and so its PSD is N 2 729nm times the PSD of the repetition frequency, as it is shown and proved in [START_REF] Puppe | Characterization of a DFG comb showing quadratic scaling of the phase noise with frequency[END_REF] for the same type of OFC, the PSD of the beat note can be measured directly.

For comparison to the cross-correlation data, this measured PSD is then mathematically divided by N 2 729nm .

IV. RESULTS AND DISCUSSION

The PSD of the phase and amplitude noise of the repetition frequency f rep of the OFC under test is presented in Fig. 5. Curve (A) shows the phase noise PSD measured by the direct optical measurement technique. The curves (C) and (B) present the PSD of the measured phase noise for two different frequency resolutions and different number of averages used in cross correlation analysis. For the curve (C), since the phase-noise values are higher for the offset frequencies close to the carrier (below 10 kHz), a lower number of averages (2 × 10 5 ) is sufficient to measure the real phase-noise value, the frequency resolution is then 150 Hz. In contrary, since the phase-noise values for the offset frequencies higher than 10 kHz are much lower than those of frequencies closer to zero Hz, the number of averages needs to be increased to 5 × 10 6 in order to reach the phase noise values of repetition frequency. The resolution is then 5 kHz. The problem that arises when using a higher number of averages is that the frequency resolution will be less due to the finite number of samples. It appears that an adaptation of the number of averages as a function of the noise frequency range is essential. As can be seen, the PSDs obtained from the cross correlation technique and the indirect measurement result are identical for offset frequencies higher than 1 kHz. Phase noise value of -157 dB rad 2 /Hz is measured at 10 kHz. The differences below 1 kHz are a consequence of harmonics of the carrier signal (250 Hz) at 500 Hz, 750 Hz and 1 kHz in the cross correlation measurement chain which originate from the frequency down conversion and ADC sampling process. Curve (D) presents the PSD of the amplitude noise measured using the cross correlation technique. The impact of the harmonics can be seen on this curve, too. The amplitude noise values for the In order to illustrate the correlation between the amplitude and phase noise a coherence function is calculated using the following equation [START_REF] Tsuchida | Correlation between amplitude and phase noise in a modelocked Cr:LiSAF laser[END_REF]:

γ 2 (f ) = S * φ (f )S φ (f ) S (f )S φ (f ) (3)
where γ 2 is a normalized parameter such that the values 1 and 0 mean perfect and no correlation, respectively. In this equation, S (f ), S φ (f ), and S φ (f ), are the PSDs of the amplitude noise, phase noise and the cross spectral density between the amplitude and phase noise, respectively.

The coherence function result is presented in Fig. 6. As can be seen, there is no significant coherence between the phase and amplitude noise for frequencies higher than 1 kHz. This lack of coherence confirms that the amplitude noise has no impact on the phase noise due to its very low level. On the contrary, relatively large peaks can be found at harmonics of 250 Hz. The strong coherence at these frequencies confirms that these harmonics on amplitude and phase noise have the same origin, as has been observed on Fig. 5 for both curves (C) and (D). Since these harmonics are missing in the optical measurement results, there is evidence that they are artifacts due to the frequency down-conversion and sampling process in the cross-correlation chain.

V. CONCLUSION

In this paper, a complete offline cross-correlation phase and amplitude noise measurement technique is presented to characterize a CEO free optical frequency comb. The output signal of the device under test is sampled by standard data acquisition cards and all the measurement processes such as phase detection, amplitude detection, and cross correlation are performed by post processing. The cost and complexity of the proposed technique are lighter compared to commercial digital and analog techniques. The accuracy of the measurement results is verified using an optical measurement technique and a very good agreement between the results is observed. Cross correlation analysis is also applied to observe the coherence between the amplitude and phase noise measurements, showing a high degree of independence for frequencies higher than 1 kHz.

Fig. 1 :

 1 Fig. 1: Block diagram of offline noise detector base on cross correlation, LPF: Low Pass Filter, ADC: Analog to Digital Converter, HT: Hilbert Transfer.

Fig. 2 :

 2 Fig. 2: PSD of the injected and detected phase and amplitude noise.

Fig. 3 :

 3 Fig. 3: Experimental set-up of the noise measurement system, OFC: Optical Frequency Comb, BPF: Band Pass Filter, PD: Phase Detector, ADC: Analog to Digital Converter, OCXO: Oven-Controlled Crystal Oscillator, PZT: Piezoelectric Transducer, LO: Local Oscillator.

Fig. 4 :

 4 Fig. 4: Experimental set-up of the optical noise measurement system, OFC: Optical Frequency Comb, BPF: Band Pass Filter, PD: Phase Detector, DSO: Digital Storage Oscilloscope, OCXO: Oven-Controlled Crystal Oscillator, EOM: Electro Optic Modulator, LO: Local Oscillator, ULE: Ultra Low Expansion.

Fig. 5 :

 5 Fig. 5: Power spectral density of the phase noise for (A) Optical measurement (B) Cross correlation 5.10 6 average (C) Cross correlation 2.10 5 average and (D) the amplitude noise results 2.10 5 average.

Fig. 6 :

 6 Fig. 6: Coherence function between amplitude and phase noise measurements.

ACKNOWLEDGMENTS

This work has been carried out thanks to the support of the A*MIDEX project (n • ANR-11-IDEX-0001-02), EquipEx Refimeve+ (n • ANR-11-EQPX-0039), and Labex First-TF (n • ANR-10-LABX-48-01), all funded by the "Investissements d'Avenir" French Government program, managed by the French National Research Agency (ANR).