Jean-Claude Fort

Thierry Klein

Agnès Lagnoux

Béatrice Laurent

Estimation of the Sobol indices in a linear functional multidimensional model

Keywords: Mathematics Subject Classification: Keywords Karhunen-Loève expansion, fractional Gaussian process, semiparametric efficient estimation, sensitivity analysis, quadratic functionals

We consider a functional linear model where the explicative variables are stochastic processes taking values in a Hilbert space, the main example is given by Gaussian processes in L 2 ([0, 1]). We propose estimators of the Sobol indices in this functional linear model. Our estimators are based on U -statistics. We prove the asymptotic normality and the efficiency of our estimators and we compare them from a theoretical and practical point of view with classical estimators of Sobol indices.

Introduction

Many mathematical models encountered in applied sciences involve a large number of poorly-known parameters as inputs. It is important for the practitioner to assess the impact of this uncertainty on the model output. An aspect of this assessment is sensitivity analysis, which aims to identify the most sensitive parameters, that is, parameters having the largest influence on the output. In global stochastic sensitivity analysis (see for example [START_REF] Oakley | Probabilistic sensitivity analysis of complex models: a Bayesian approach[END_REF] and [START_REF] Saltelli | Sensitivity analysis[END_REF] and references therein) the input variables are assumed to be independent random variables. Their probability distributions account for the practitioner's belief about the input uncertainty. This turns the model output into a random variable, whose total variance can be split down into different partial variances (this is the socalled Hoeffding decomposition see [START_REF] Van Der | Asymptotic statistics[END_REF]). Each of these partial variances measures the incertitude on the output induced by each input variable uncertainty. By considering the ratio of each partial variance to the total variance, we obtain a measure of importance for each input variable that is called the Sobol index or sensitivity index of the variable [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]; the most sensitive parameters can then be identified and ranked as the parameters with the largest Sobol indices.

Once the Sobol indices have been defined, the question of their effective computation or estimation remains open. In practice, one has to estimate (in a statistical sense) those indices using a finite sample of evaluations of model outputs [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. Many Monte Carlo or quasi Monte Carlo approaches have been developed by the experimental sciences and engineering communities. This includes the FAST methods (see for example [START_REF] Ri Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF], [START_REF] Tissot | A bias correction method for the estimation of sensitivity indices based on random balance designs[END_REF] and references therein) and the Sobol pick-freeze (SPF) scheme (see [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF][START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]). Nevertheless, those methods require many evaluations of model outputs which can be a strong limitation when thoses evaluations are expensive. Many approaches have been developped to overcome this issue. The most popular are Bayesian approach (see for example [START_REF] Oakley | Probabilistic sensitivity analysis of complex models: a Bayesian approach[END_REF]) or the construction of metamodels. As mentioned in Kleijnen [START_REF] Jack | Sensitivity analysis and related analyses: A review of some statistical techniques[END_REF] (see equation ( 1) page 121) one can use functional linear regression as metamodel.

In this paper, we study the particular context of the functional linear regression and propose a different way of estimation. We consider nonparametric estimators of quadratic functionals by projection methods, which are related to the procedures developed by Laurent (see [START_REF] Laurent | Efficient estimation of integral functionals of a density[END_REF][START_REF] Laurent | Adaptive estimation of a quadratic functional of a density by model selection[END_REF]) in a density model and by Da Veiga and Gamboa in [START_REF] Veiga | Efficient estimation of nonlinear conditional functionals of a density[END_REF] in a regression model. This method allows us to estimate simultaneously all the Sobol indices with a single sample of reasonable size.

More precisely we consider a separable Hilbert space H endowed with the scalar product <, > and X 1 , . . . , X p , p independent centered, H-valued, stochastic processes. The model that we consider is a linear regression model :

Y = µ + p k=1 < β k , X k > +ε.
where β i , 1 ≤ i ≤ p are elements of H, µ is in R and ε is a centered noise independent of the processes X 1 , . . . , X p .

Our approach is based on the so-called Karhunen-Loève decomposition of the processes X k ( [START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF][START_REF] Bronski | Small ball constants and tight eigenvalue asymptotics for fractional Brownian motions[END_REF]). Thanks to this decomposition we construct natural estimators of the Sobol indices for whom we prove asymptotic normality and efficiency. (Asymptotic efficiency is a natural property which generalizes the notion of minimum variance unbiased estimator, see [START_REF] Van Der | Asymptotic statistics[END_REF] chapters 8 and 25 or [START_REF] Ibragimov | Statistical estimation-asymptotic theory[END_REF] for more details.)

This paper is organized as follows: in the first section, we set up the notations for the model, review the definition of Sobol indices and present our estimators. In the second section, we prove asymptotic normality and efficiency when we consider a simple functional linear regression model. These two properties are generalized in the third section in the general setting of the multiple functional linear regression. In Section 3, we also compare this method with the classical SPF. The fourth section gives numerical illustrations on a benchmark model.

Setting and notations

Let H be a separable Hilbert space endowed with the scalar product <, > and let X 1 , . . . , X p be p independent centered, H-valued, stochastic processes.

The generic model that we consider is the following

Y = µ + p k=1 < β k , X k > +ε, (1)
where β k , 1 ≤ k ≤ p are elements of H, µ is in R and ε is a centered noise with variance η 2 independent of the processes X 1 , . . . , X p , and with finite fourth order moment.

We define, as usual in a finite dimensional setting, the Sobol index with respect to the entry (explicative variable) number k :

S k = Var(E(Y |X k )) Var(Y ) .
From the observation of X 1 i , . . . , X p i , Y i 1≤i≤n i.i.d. obeying to Model (1), our aim is to estimate the vector S 1 , . . . , S p . Since Var(Y ) can be easily estimated by the empirical variance based on (Y 1 , . . . , Y n ), the main purpose of the paper is to estimate for all k the quantity

V k = Var(E(Y |X k )).
We will assume that we know the distributions of the input processes X k 1≤k≤p . In the next section we consider the simple case where p = 1. In this setting the Sobol index is of less interest, but the computations then easily extend to the generic model.

Simple functional linear regression model

Using the same notations as in Section 1, we consider, in this section, the case p = 1, which leads to the model

Y = µ+ < β , X > +ε. (2)
We observe a n-sample of (X, Y ), that we denote

(X i , Y i ), 1 ≤ i ≤ n. We have E(Y |X) = µ+ < β, X > and V X = Var (E(Y |X)) = Var(< β , X >).
The Sobol index S X is defined by

S X = V X Var(Y )
.

We assume that E X 2 < ∞, so that the covariance operator defined for all f ∈ H by Γ(f ) = E [< X, f > X] is Hilbert-Shmidt and is diagonalizable via the Karhunen-Loève expansion in an orthonormal basis of eigenfunctions (ϕ l , 1 ≤ l), with decreasing eigenvalues (λ l , 1 ≤ l) such that ∞ l=1 λ l < ∞. See e.g. [START_REF] Loève | Fonctions aléatoires de second ordre[END_REF][START_REF] Loève | Probability theory. I[END_REF][START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF]. We set < X, ϕ l >= √ λ l ξ l . The variables (ξ l ) l≥1 are centered, uncorrelated, with variance 1. When X is a Gaussian process, the variables (ξ l ) l≥1 are i.i.d. standard Gaussian variables.

λ l = 1 (π(l -1 2 )) 2 and ϕ l (t) = √ 2 sin( t √ λ l ), l ≥ 1. 2. Let W = (W t ) t∈[0,1] a Brownian bridge and H = L 2 ([0, 1], dt). λ l = 1 (πl) 2 and ϕ l (t) = √ 2 sin(πlt), l ≥ 1.
3. Let X = (X t ) t∈[0,1] a fractional Brownian motion with Hurst exponent H.

Then the asymptotic of the eigenvalues is given by

λ l = sin(πH)Γ(2H + 1) (lπ) 2H+1 + o l (2H+2)(4H+3) 4H+5
+δ .

In the case H = 1/2, the Brownian case, the result agrees with the exact result

λ l = 1 (π(l-1/2)) 2 .
In the basis (ϕ l ) l≥1 , β has the expansion β = ∞ l=1 γ l ϕ l , and we have :

E(Y X) = E(< X, β > X) = Γ(β) = ∞ l=1 λ l γ l ϕ l .
Setting g = ∞ l=1 λ l γ l ϕ l , we have : γ l =< g, ϕ l > /λ l . The index V X is then given by

V X = Var(E(Y |X)) = E(E(Y |X) 2 ) = E(< β, X >< β, X >) =< β, Γ(β) > . which expends V X = ∞ l=1 λ l γ 2 l .
We consider the empirical unbiased estimator of γ l :

γ l = 1 λ l 1 n n i=1 < X i , ϕ l > Y i .
For all m ∈ N * , we introduce the U -statistics of order 2 :

V X m = m l=1 1 λ l 1 n(n -1) 1≤i =j≤n < X i , ϕ l > Y i < X j , ϕ l > Y j . (3)
We have

E V X m = m l=1 λ l γ 2 l , hence V X m is a biased estimator of V X , with bias B m = ∞ l=m+1 λ l γ 2 l .
In the following section, we study the asymptotic behaviour of V X m , for a suitable choice of the parameter m.

Central limit theorem

Using Hoeffding's decompsition of the U -statistics V X m (see [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF]), we straightforwardly get the following result.

Proposition 1. V X m can be rewritten as the sum of a totally degenerated Ustatistics of order 2, a centered linear term and a deterministic term in the following way :

V X m = U n K + P n L + m l=1 γ 2 l λ l ( 4
)
where

U n K := m l=1 1 λ l 1 n(n -1) 1≤i =j≤n (Y i < X i , ϕ l > -γ l λ l ) (Y j < X j , ϕ l > -γ l λ l ) P n L := 2 n m l=1 n i=1 γ l (Y i < X i , ϕ l > -γ l λ l ) .
As a consequence, we have

V X m -V X = U n K + P n L - l>m γ 2 l λ l = U n K + P n L -B m . Theorem 2. Let (X i , Y i ) 1≤i≤n be i.i.d.
observations with the same distribution as (X, Y ) from Model (2). We assume that E( X 4 ) < +∞ and that E(ε 4 ) < +∞. We consider the Karhunen-Loève expansion of X :

X = l≥1 λ l ξ l ϕ l .
We assume

sup l≥1 E(ξ 4 l ) < +∞. (5)
We consider the estimator V

X m of V X defined by (3) with m = m(n) = √ nh(n), where h(n) satisfies : h(n) → 0 and ∀α > 0, n α h(n) → +∞ as n → +∞.
We assume that there exist C > 0 and δ > 1 such that

∀l ≥ 1, λ l ≤ Cl -δ Then we have √ n( V X m -V X ) L → n→∞ N (0, 4Var(Y < X, β >)). (6)
Comments :

1. We may assume that h(n) = 1/ log(n), and hence

m(n) = √ n/log n, to fill the condition ∀α > 0, lim n→∞ n α h(n) = +∞.
The estimator V X m converges at the parametric rate 1/ √ n, for any β.

2. The nature of the problem of estimating the quadratic functional V X m is completely different from the estimation of the signal β, where nonparametric rates are obtained (see for example [START_REF] Comte | Adaptive estimation in circular functional linear models[END_REF] for the estimation of β in a circular functional linear model).

3. We will prove in the next section the asymptotic efficiency of the estimator V X m . 4. Note that the condition ( 5) is verified when X is a Gaussian process, since in the case, the variables (ξ l ) l≥1 are i.i.d. standard Gaussian variables.

5. In Theorem 2 we have assumed that there exist C > 0 and δ >

1 such that ∀l ≥ 1, λ l ≤ Cl -δ .
Let us recall that since E X 2 < +∞, l≥1 λ l < +∞ hence this assumption is not very strong.

Proof. In order to prove (6), we will show that

     B 2 m = o (1/n) E (U n K) 2 = o (1/n) √ nP n L L → n→∞ N (0, 4Var(Y < X, β >))
.

a) Bias term

Since we have assumed that ∀l ≥ 1, λ l ≤ Cl -δ for some C > 0 and δ > 1, and since l≥1 γ 2 l < +∞, we get

B m = l>m γ 2 l λ l ≤ C l>m γ 2 l l -δ ≤ m -δ l>m γ 2 l .
Recalling the definition of m = m(n), we obtain

B 2 m = o 1 n . b) Term U n K One has E((U n K) 2 ) = E 1 + E 2 + E 3 where E 1 = 2 (n(n -1)) 2 m l1,l2=1 1 λ l1 λ l2 1≤i1 =j1≤n E (Y i1 < X i1 , ϕ l1 > -λ l1 γ l1 ) (Y i1 < X i1 , ϕ l2 > -λ l2 γ l2 )(Y j1 < X j1 , ϕ l1 > -λ l1 γ l1 )(Y j1 < X j1 , ϕ l2 > -λ l2 γ l2 ) , E 2 = 4 (n(n -1)) 2 m l1,l2=1 1 λ l1 λ l2 i1,j1,j2all = E (Y i1 < X i1 , ϕ l1 > -λ l1 γ l1 ) (Y i1 < X i1 , ϕ l2 > -λ l2 γ l2 )(Y j1 < X j1 , ϕ l1 > -λ l1 γ l1 )(Y j2 < X j2 , ϕ l2 > -λ l2 γ l2 ) , E 3 = 1 (n(n -1)) 2 m l1,l2=1 1 λ l1 λ l2 i1,j1,i2,j2all = E (Y i1 < X i1 , ϕ l1 > -λ l1 γ l1 ) (Y i2 < X i2 , ϕ l2 > -λ l2 γ l2 )(Y j1 < X j1 , ϕ l1 > -λ l1 γ l1 )(Y j2 < X j2 , ϕ l2 > -λ l2 γ l2 ) .
One easily see that

E 2 = 0 and E 3 = 0,
since the variables (Y < X, ϕ l > -λ l γ l ) are centered and independent.

Let us now compute

E 1 . Denote Zj,k := Y j < X j , ϕ k > -λ k γ k . E 1 = 2 (n(n -1)) 2 m l,k=1 1 λ l λ k 1≤i =j≤n E Z1,l Z2,l Z1,k Z2,k = 2 n(n -1) m l,k=1 1 λ l λ k E Z1,l Z2,l Z1,k Z2,k = 2 n(n -1) m l,k=1 1 λ l λ k E 2 Z1,l Z1,k ≤ 2 n(n -1) m l,k=1 1 λ l λ k E (Y < X, ϕ l >) 2 E (Y < X, ϕ k >) 2 = 2 n(n -1) m l,k=1 E Y 2 ξ 2 l E Y 2 ξ 2 k ≤ 2 n(n -1) E Y 4 m l,k=1 E ξ 4 l 1/2 E ξ 4 k 1/2 ≤ 2 n(n -1) E Y 4 m l=1 E ξ 4 l 1/2 2
By assumption [START_REF] Gamboa | Statistical inference for sobol pick freeze monte carlo method[END_REF] we know that sup k E(ξ 4 k ) ≤ K, hence we have

E 1 ≤ 2Km 2 n(n -1) E Y 4
and we still obtain that

E (U n K) 2 = o 1 n . c) Term P n L First, let P n L := 2 n n i=1 Y i < X i , β > -E(Y < X, β >) (7) = 2 n n i=1   Y i < X i , β > - l≥1 γ 2 l λ l   (8) since E(Y < X, β >) = E(< X, β > 2 ) = E   ∞ l=1 γ l λ l ξ l 2   = ∞ l=1 γ 2 l λ l ,
we write

P n L = 2 n m l=1 n i=1 γ l (Y i < X i , ϕ l > -γ l λ l ) = 2 n n i=1 Y i < X i , β > m - m l=1 γ 2 l λ l = P n L + (P n L -P n L ),
where

< X i , β > m denotes m l=1 γ l < X i , ϕ l >.
In the one hand,

P n L -P n L = 2 n n i=1 Y i [< X i , β > -< X i , β > m ] -2 l>m γ 2 l λ l = 2 n n i=1 l>m γ l [Y i < X i , ϕ l > -γ l λ l ] Let Z := l>m γ l [Y < X, ϕ l > -γ l λ l ]
and let us give an upper bound for Var(Z).

Var(Z) = Var l>m γ l [Y < X, ϕ l > -γ l λ l ] = Var l>m γ l Y < X, ϕ l > ≤ E   Y 2 l>m γ l < X, ϕ l > 2   = E Y 2 < X, β > 2 m ⊥ ≤ E Y 2 ||X|| 2 m ⊥ ||β|| 2 m ⊥ by Cauchy-Schwarz inequality ≤ E Y 2 ||X|| 2 ||β|| 2 m ⊥ ,
where < X, β > m ⊥ and |β|| m ⊥ respectively denote l>m γ l < X i , ϕ l > and l>m γ 2 l . The first expectation does not depend on n and the second term tends to 0 when n (and therefore m) tends to ∞. As a consequence,

• Var(Z) = o(1), • Var(P n L -P n L ) = o 1 n , • P n L -P n L = o P 1 √ n .
In the other hand, we establish a central limit theorem for P n L , since it is an empirical sum of iid centered variables. As a conclusion,

√ nP n L L → n→∞ N (0, 4Var(Y < X, β >)) (9)
and by Slutsky theorem since

P n L -P n L = o P 1 √ n , √ nP n L L → n→∞ N (0, 4Var(Y < X, β >)). (10)
This concludes the proof.

Our aim is to estimate the Sobol index S X which is defined by S X = V X /Var(Y ).

It is therefore natural to introduce the estimator of S X defined by V X m / V , where V is the empirical variance estimating Var(Y ):

V := 1 n -1 n i=1 (Y i -Ȳn ) 2 .
In the following theorem, we obtain the asymptotic behaviour of this estimator.

Theorem 3. Under the same assumptions as in Theorem 2, we have

√ n V X m V -S X L → n→∞ N 0, Var(U ) (Var(Y )) 2 (11)
where

U := 2Y < X, β > -S X (Y -E(Y )) 2 . Proof. (i) First, let W i = (2 [Y i < X i , β > -E(Y < X, β >)] , Y 2 i , Y i ) t (i = 1, . . . , n) and W = (2 [Y < X, β > -E(Y < X, β >)] , Y 2 , Y ) t . Then W n = 1 n n i=1 W i = (P n L , Y 2 n , Y n ) and E W := E(W ) = (0, E(Y 2 ), E(Y )).
Now from the Central Limit Theorem,

√ n W n -E W L → n→∞ N 3 (0, Σ) with Σ =   4Var(Y < X, β >) 2Cov(Y i < X, β >, Y 2 ) 2Cov(Y < X, β >, Y ) 2Cov(Y < X, β >, Y 2 ) Var(Y 2 ) Cov(Y 2 , Y ) 2Cov(Y < X, β >, Y ) Cov(Y 2 , Y ) Var(Y )   (ii) Then let W i = W i + (U n K + (P n L -P n L ) + m l=1 γ 2 l λ l , 0, 0). We have W n = 1 n n i=1 W i = ( V X m , Y 2 n , Y n ). Since, U n K + (P n L -P n L ) + l>m γ 2 l λ l = o P 1 n , we still have a Central Limit Theorem √ n W n -E W L → n→∞ N 3 (0, Σ)
where

E W = (E(Y < X, β >), E(Y 2 ), E(Y )).
(iii) Let Φ the mapping from R 3 to R defined by Φ(x, y, z) = x y-z 2 . We want to apply the Delta method (cf. [START_REF] Van Der | Asymptotic statistics[END_REF]) to W and Φ. Easily,

       ∂Φ ∂x (x, y, z) = 1 y-z 2 ∂Φ ∂x (E W ) = 1 Var(Y ) ∂Φ ∂y (x, y, z) = -x (y-z 2 ) 2 ∂Φ ∂y (E W ) = -E(Y <X,β>) (Var(Y )) 2 ∂Φ ∂z (x, y, z) = 2xz (y-z 2 ) 2 ∂Φ ∂x (E W ) = 2E(Y <X,β>)E(Y ) (Var(Y )) 2 so that √ n Φ(W n ) -Φ(E W ) L → n→∞ N 0, Φ (E W )ΣΦ (E W ) t . But Φ(W n ) = V X m V , Φ(E W ) = S X and Φ (E W )ΣΦ (E W ) t = 4 Var(Y < X, β >) (Var(Y )) 2 + 2 E(Y < X, β >) (Var(Y )) 3 4E(Y )Cov(Y < X, β >, Y ) -2Cov(Y < X, β >, Y 2 ) + 2(E(Y )) 2 E(Y < X, β >) + (E(Y < X, β >)) 2 (Var(Y )) 4 Var(Y 2 ) -4E(Y )Cov(Y 2 , Y )
One can check that

Φ (E W )ΣΦ (E W ) t = Var(U ) (Var(Y )) 2 .
Remark 2.1.

1. In the case where the variance of Y is known and plug in the estimator, the result is

√ n V X m Var(Y ) -S X L → n→∞ N (0, 4 Var(Y < X, β >) (Var(Y )) 2 ). (12)
2. In the general case where the variance of Y is unknown and Y is centered, we thus have an extra term (due to the estimation of the variance), namely the difference between the two asymptotic variances is

Var(U ) -4Var(Y < X, β >) (Var(Y )) 2 = E(Y < X, β >) (Var(Y )) 4 E(Y < X, β >)Var(Y 2 ) -4Var(Y )Cov(Y < X, β >, Y 2 )
It is worth to notice that this term is not always positive. Namely, if < X, β > is a centered variable with second moment s 2 and fourth moment k 1 s 4 and if ε is a centered variable with second moment η 2 and fourth moment k 2 η 4 , then the latter extra term is

- 6s 4 (s 2 + η 2 ) 2 + s 4 (s 2 + η 2 ) 4 ((3 -k 1 )(3s 4 + 4η 2 s 2 ) -(3 -k 2 )η 4 ).
In particular, when < X, β > L ∼ N (0, s 2 ) and

L ∼ N (0, η 2 ), k 1 = k 2 = 3
and one gets an extra term equal to -6s 4 (s 2 + η 2 ) 2 (Var(Y )) 4 .

In this particular case, the asymptotic variance of the estimator of the Sobol index is smaller when the variance of Y is estimated that when it is known and plug in the estimator.

Asymptotic efficiency

In this section we prove that V X m / V is asymptotically efficient for estimating the Sobol index S X (see [START_REF] Van Der | Asymptotic statistics[END_REF], Section 25 for the definition of asymptotic efficiency). This notion somewhat extends the notion of Cramér-Rao bound and enables to define a criterion of optimality for estimators, called asymptotic efficiency.

Proposition 2.1 (Asymptotic efficiency). Under the same assumptions as in Theorem 2, V X m / V is asymptotically efficient for estimating S. Proof. Note that

S = Φ(V X , Var(Y )) and V X m V = Φ( V X m , V )
where Φ is defined by Φ(x, y) = x/y. First let us prove that V X m is asymptotically efficient for estimating V X . The result directly comes from decomposition (4) and the Hoeffding's decomposition V X m -V X = U n K + P n L -B m , P n L is asymptotically efficient (example 25.24 [efficiency of the empirical distribution]), the other terms are o P 1 √ n . From which we conclude the asymptotic efficiency of V X m .

Second using again example 25.24, we get the asymptotic efficiency of V for estimating Var(Y ).

Then theorem 25.50 (efficiency in product space) in [START_REF] Van Der | Asymptotic statistics[END_REF] gives that ( V X m , V ) is asymptotically efficient for estimating (V X , Var(Y )). Now since Φ is differentiable in R 2 \ {y = 0} then by theorem 25.47 (efficiency and Delta method), (Φ( V X m , V )) n is asymptotically efficient for estimating Φ(V X , Var(Y )).

Multiple functional linear regression model

We recall the model (1):

Y = µ + p k=1 < β k , X k > +ε.
In this setting we observe a n-sample (Y i , X 1 i , . . . , X p i ) of (Y, X 1 , . . . , X p ).

Estimation by U-statistics of order 2

We assume that the processes X k 1≤k≤p are independent. Under this assumption, the previous section whole applies and the Sobol indices can be estimated from the U statistics of order 2. For all m ≥ 1, let

V k m = m l=1 1 λ k l 1 n(n -1) 1≤i =j≤n < X k i , ϕ l > Y i < X k j , ϕ l > Y j . For a suitable choice of m = m(n), V k m is a consistent estimators of V k . Actually by independence E(Y |X k ) =< β k , X k >,
and denoting ε = p j=1,j =k < β j , X j > +ε, we retrieve the simple regression model. In fact we can obtain a central limit theorem for V m where

V m = ( V 1 m , V 2 m , . . . , V p m ) T . ( 13
) Let V j = Var(E(Y |X j )), j = 1 . . . p and V = (V 1 , V 2 , . . . , V p ) T .
Theorem 4. Let (Y i , X 1 i , . . . , X p i ) be i.i.d. observations from Model (1). We assume that for all 1 ≤ j ≤ p, E( X j 4 ) < +∞ and that E(ε 4 ) < +∞. We consider the Karhunen-Loève expansion of X j :

X j = l≥1 λ j l ξ j l ϕ j l .
We assume

∀1 ≤ j ≤ p, sup l≥1 E((ξ j l ) 4 ) < +∞. ( 14
)
We consider the estimator V m of V defined by ( 13)

with m = m(n) = √ nh(n), where h(n) satisfies : h(n) → 0 as n → +∞ and ∀α > 0, n α h(n) → +∞ as n → +∞. We assume that ∀ 1 ≤ j ≤ p, ∀l ≥ 1, λ j l ≤ C j l -δj for some C j > 0 and δ j > 1.
The following result holds :

√ n( V m -V ) L → n→∞ N p 0, 4 Cov(Y < X k , β k >, Y < X l , β l >) k,l=1...p . ( 15
)
Proof. Using the same decomposition

V m -V = U n K + P n L -B m
and the fact that U n K , P n L-P n L and B m 2 are o P 1 n and o 1 n , we just need to check that P n L = (P n L 1 , P n L 2 , . . . , P n L p ) T converges in distribution towards a N p 0, 4 Cov(Y < X k , β k >, Y < X l , β l >) k,l=1...p . We have

P n L = 2 n n i=1 W i
where

W i = (W 1 i , . . . , W p i ) T and W k i = Y i < X k i , β k > -E(Y i < X k i , β k >).
Since the W i 's are i.i.d., the result follows from the central limit theorem in

R p . Now define S m as V m / V = ( V 1 m / V , . . . , V p m / V ) T
, where V ) denotes the empirical estimator of the Var(Y ). Let us remind that S = V /Var(Y ) = (S 1 , . . . , S p ) T . Using once again the Delta method, we easily get the following result. Proposition 3.1. Under the same assumptions as in Theorem 4, we have

√ n S m -S L → n→∞ N p 0, Cov(U k , U l ) (Var(Y )) 2 k,l=1...p . ( 16
)
where

U k := 2(Y -E(Y )) < X k , β k > -S k (Y -E(Y )) 2 .

Comparison with Sobol estimators

Sobol [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] proposed an empirical method, based on a particular design of experiments, to estimate Sobol indices. This method called Sobol Pick and Freeze (SPF) is also studied in [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF]. In [START_REF] Janon | Asymptotic normality and efficiency of a two sobol index estimators[END_REF][START_REF] Gamboa | Statistical inference for sobol pick freeze monte carlo method[END_REF], the authors prove asymptotic normality and efficiency of Sobol estimators. In this section, our aim is to compare the method proposed in this paper, which is based on U statistics of order 2, and the method proposed by Sobol in order to see which experimental set up should be used by the practitionner. It is important to notice that the two procedures are not based on the same design of experiments. Both methods lead to asymptotically normal and efficient estimators, but the asymptotic variances are different due to the designs of experiments. We therefore want to compare the asymptotic variances obtained by the two methods, for a similar number of experiments in both cases. Let us first recall the design of experiments and the estimators proposed by Sobol.

Let (X, Y ) obey to the model ( 1), where X = (X 1 , . . . , X p ). Let X = (X 1 , . . . , X p ) an i.i.d. copy of X. For all k ∈ {1, . . . , p}, let Y k be defined by

Y k =< X k , β k > + p j=1,j =k < X j , β j > +ε , where ε is an i.i.d. copy of ε. Let (Y 1 , . . . , Y n ) be i.i.d. copies of Y and for all k ∈ {1, . . . , p} let (Y k 1 , . . . , Y k n ) be i.i.d. copies of Y k . The estimator proposed by Sobol to estimate V k = Var(E(Y |X k )) is defined by V k SP F = 1 n n i=1 Y i Y k i - 1 n(p + 1) n i=1 Y i + Y 1 i + . . . + Y p i 2 . Now let S SP F = V SP F / V = ( V 1 SP F / V , . . . , V p SP F / V ) T
, where V is the estimator of Var(Y ) in the S.P.F. method defined by

V = 1 n(p + 1) n i=1 (Y i ) 2 + Y 1 i 2 + . . . + (Y p i ) 2 - 1 n(p + 1) n i=1 Y i + Y 1 i + . . . + Y p i 2 .
Theorem 5. We have (see [START_REF] Janon | Asymptotic normality and efficiency of a two sobol index estimators[END_REF][START_REF] Gamboa | Statistical inference for sobol pick freeze monte carlo method[END_REF])

√ n S SP F -S L → n→∞ N p 0, Cov(V k , V l ) (Var(Y )) 2 k,l=1...p . ( 17
)
where

V k := Y Y k -S k Y 2 + p k=1 (Y k ) 2 /(p + 1).
The asymptotic variance appearing in the central limit theorem for the estimator 15)). Our aim is to compare Γ and Γ SP F .

S SP F = ( S 1 SP F , . . . , S p SP F ) (resp. S m = ( S 1 m , . . . , S p m )) is denoted by (Var(Y )) -2 Γ SP F (resp. (Var(Y )) -2 Γ (cf. (
For sake of simplicity we assume ε = 0. Then we define

W i =< X i , β i >, σ 2 i = Var(W i ), σ 2 = p i=1 σ 2 i and S i = σ 2 i /σ 2 . After computations, one obtains that Γ is defined by ∀k, Γ(k, k) = σ 4 k σ 4 p i=1 E(W 4 i ) -3σ 4 i + 4 σ 2 -σ 2 k σ 2 E(W 4 k ) + 4σ 2 σ 2 k + 12 σ 6 k σ 2 -14σ 4 k ∀k = l, Γ(k, l) = - 2 σ 2 E(W 4 k )σ 2 l + E(W 4 l )σ 2 i + 2σ 2 k σ 2 l 3 σ 2 k + σ 2 l σ 2 -1 + σ 2 k σ 2 l σ 4 p i=1 E(W 4 i ) -3σ 4 i and Γ SP F by ∀k, Γ SP F (k, k) = 1 - 4 p + 1 σ 2 k σ 2 E(W 4 k ) + σ 4 k σ 4 p 2 -2p + 5 (p + 1) 2 p i=1 E(W 4 i ) + σ 2 k (p + 1)σ 2 4 + σ 2 k σ 2 -2p 2 + 5p -13 p + 1 p i=1 σ 4 i +σ 4 + -p 2 -7p + 6 (p + 1) 2 σ 4 k + 4 p + 1 σ 6 k σ 2 - 4 p + 1 σ 2 k σ 2 ∀k = l, Γ SP F (k, l) = - 2 p + 1 σ 2 k E(W 4 l ) + σ 2 l E(W 4 k ) σ 2 + σ 2 k σ 2 l σ 4 p 2 -2p + 5 (p + 1) 2 p i=1 E(W 4 i ) + 2(σ 2 k + σ 2 l ) + -2p 2 + 5p -13 (p + 1) 2 σ 4 p i=1 σ 4 i +σ 4 + σ 2 k σ 2 l 2p 2 -p + 9 (p + 1) 2 - p + 3 p + 1 σ 2 (σ 2 k + σ 2 l ) + 2 p + 1 σ 2 k σ 2 l σ 2 (σ 2 k + σ 2 l ).
D := (p + 1)Γ SP F -Γ.
In order to compare the two methods, we evaluate the eigenvalues of the matrix D in order to determine whether it is positive-definite or not. If not we study only the sign of the diagonal terms of D that correspond to the difference of the asymptotic variances obtained with both methods.

First example: We asssume that for i = 1, . . . , p with p ≥ 2 < X i , β i >∼ N (0, 1).

We have

D i,i = p 5 + 2p 4 -5p 3 -2p 2 + 6p -6 p(p + 1) > 0 ∀ 1 ≤ i ≤ p, D i,j = p 5 + 4p 4 + 3p 3 + 2p -6 p(p + 1) ∀ 1 ≤ i = j ≤ p.
Let α = p 5 + 2p 4 -5p 3 -2p 2 + 6p -6, β = p 5 + 4p 4 + 3p

+ p 3 -8p 2 -2) > 0 ∀p ≥ 2.
For p ≥ 2, the matrix D is not positive nor negative. Hence, as mentioned before, we simply study the sign of the diagonal terms of the matrix D. Since D i,i > 0 for any p ≥ 2 and i = 1 . . . p, the asymptotic variance of each estimator of Sobol indices with SPF method is larger than the one obtained by the method proposed in this paper using the same number of observations.

Second example : In the case p = 2 and

σ 2 = σ 2 1 + σ 2 2 = 1, the model is Y =< X 1 , β 1 > + < X 2 , β 2 > . Then if W 1 =< X 1 , β 1 > and W 2 =< X 2 , β 2 >, D 1,1 = E(W 4 1 ) -1 + 2 3 σ 4 1 + 2 3 σ 4 1 E(W 4 2 ) + 3 -8σ 6 1 + 10σ 4 1 -8σ 2 1 + σ 2 1 (σ 4 1 + σ 4 2 ) 4 - 2 3 σ 2 1 D 2,2 = E(W 4 2 ) -1 + 2 3 σ 4 2 + 2 3 σ 4 2 E(W 4 1 ) + 3 -8σ 6 2 + 10σ 4 2 -8σ 2 2 + σ 2 2 (σ 4 1 + σ 4 2 ) 4 - 2 3 σ 2 2 D 2,1 = σ 2 1 σ 2 2 2 3 (E(W 4 1 ) + E(W 4 2 )) + 3 + σ 4 1 + σ 4 2 7 3 + 3σ 2 1 σ 2 2 -2
In the particular case where W 1 ∼ N (0, x) and W 2 ∼ N (0, 1 -x), for some

x ∈]0, 1[, we obtain

D 1,1 (x) = 3 -4x + 1 3 x 2 -8x 3 + 8x 4 , D 2,2 (x) = 3 -4(1 -x) + 1 3 (1 -x) 2 -8(1 -x) 3 + 8(1 -x) 4 , D 1,2 (x) = 1 3 + 10 3 x - 40 3 x 2 + 20x 3 -10x 4 .
Let λ 1 (x) and λ 2 (x) be the eigenvalues of D(x) and x 1 and 1 -x 1 the real zeros of its determinant (x 1 ≈ 0.6701). Then

x 0 1 -x 1 x 1 1 λ 1 (x) + λ 2 (x) + | + | + λ 1 (x)λ 2 (x) - 0 + 0 - For x ∈]x 1 , x 2 [, the matrix D(x) is positive-definite.
We can also study the signs of the diagonal terms of the matrix D(x) for x ∈]0, 1[. Now let x 0 the zero of D 1,1 (x) (x 0 ≈ 0.6738). The folowing tabular gives the signs of the diagonal terms of the matrix D(x) :

x 0 1 -x 0 x 0 1 D 1,1 (x) + | + 0 - D 2,2 (x) - 0 + | +
We deduce from this tabular that for x ∈]1-x 0 , x 0 [, our method leads to smaller variances for both estimators.

Numerical experiments

We perform a simulation study to evaluate the performances of the procedure proposed in this paper for estimating Sobol indices in a functional linear model, and to compare this procedure with Sobol's procedure. We consider the model :

Y i = p k=1 < β k , X k i > +ε i , (18)
where for all k = 1, . . . , p, X k and β k are defined from the coefficients of their expansions onto an orthonormal basis of L 2 ([0, 1]). For all k, we consider the basis corresponding to the eigenfunctions of the Karhunen-Loève expansion of the process X k and we define the function β k by the coefficients of its expansion onto this basis. For the simulations, the basis that we consider is either the one associated to the Karhunen-Loève expansion of the Brownian motion (that has been recalled in Section 2) or the one associated to the fractional Brownian motion. It is important to recall that the processes X k and the related functions β k are expanded onto the same basis and that the estimators proposed in this paper will perform well in the case where the coefficients of the functions β k in those basis are non-increasing. The variables ε i are i.i.d. Gaussian centered variables with variance σ 2 .

First example

Using Karhunen-Loève expansion of the Brownian motion

We consider the model ( 18) with p = 2 and ε i = 0 for all i. The processes X 1 and X 2 are given by a truncated expansion of the Brownian motion onto its Karhunen-Loève basis, with coefficients:

λ l = 1 π(l -1/2) 2 , 1 ≤ l ≤ L; λ l = 0, l > L.
Concerning the functions β 1 , β 2 , they are respectively caracterized by the coefficients (γ l ) l≥1 of their expansion onto the same basis:

Conclusion

We have proposed an estimator of the sobol indices in a simple functional regression model. When the entries of the model are independent our estimator is asymptotically efficient. This is also the case of the so-called SPF estimator.

As the designs underlying the two estimators are not the same, we numerically compared the two estimators in the same conditions, that is with the same number of calls to the model. The results are strongly dependent on the behavior of the regressors. When the regressors are well represented onto the Karhune-Loève basis associated to the entries then our estimator behave well and takes advantage of the information given by the knowledge of the entries, it performs better than the SPF. A contrario, if the regressors have small coefficients on the first elements of the Karhunen-Loève basis, the SPF is more precise than our estimator, since the knowledge of the entries is of less use. These conclusions are brought in the case of a model with independent entries. One possible advantage of our estimator is that we can generalize it to the case of dependent entries, which is not the case of the SPF which deeply relies on the independence of the entries. This will be the topic of a forthcoming research paper.

Example 2 . 1 (

 21 Examples of Karhunen-Loève expansions ([START_REF] Luschgy | Functional quantization of Gaussian processes[END_REF][START_REF] Bronski | Small ball constants and tight eigenvalue asymptotics for fractional Brownian motions[END_REF])).

1 .

 1 Let B = (B t ) t∈[0,1] a Brownian motion and H = L 2 ([0, 1], dt). Then

  Sobol experiment requires (p + 1)n observations (or computations of the blackbox code) to estimate the p indices. In order to have a fair comparison of both methods, we consider that we have n(p+1) i.i.d. observations from Model (1) to estimate the Sobol indices by our methods. With n(p + 1) observations instead of n, the asymptotic variance matrix Γ m is devided by p + 1, hence, we have to study the matrix

  First Model: S = [0.5438, 0.2639, 0.0635, 0.1288] n Bias( Ŝm) Bias( ŜSP F ) 10 2 10 -3 [3.39, 1.30, 0.56, 0.08] 10 -3 [5.32, 5.41, 8.87, 4.70] 10 3 10 -4 [7.66, 3.60, 0.33, 2.07] 10 -4 [3.84, 0.59, 0.35, 5.65] [5.54, 4.28, 2.35, 3.23] 10 -2 [9.91, 9.79, 9.71, 9.62] 10 3 10 -2 [1.82, 1.36, 0.72, 0.99] 10 -2 [3.13, 3.12, 3.11, 3.06] Second Model: S = [0.7561, 0.1871, 0.0112, 0.0456] n Std( Ŝm) Std( ŜSP F ) 10 2 10 -2 [6.11, 3.72, 1.22, 2.01] 10 -1 [1.07, 1.00, 1.01, 0.99] 10 3 10 -2 [1.97, 1.17, 0.33, 0.60] 10 -2 [3.36, 3.16, 3.14, 3.13] First Model: S = [0.5438, 0.2639, 0.0635, 0.1288] n RMSE( Ŝm) RMSE( ŜSP F ) 10 2 10 -2 [5.55, 4.29, 2.35, 3.22] 10 -2 [9.92, 9.80, 9.75, 9.63] 10 3 10 -2 [1.82, 1.36, .72, 0.99] 10 -2 [3.13, 3.12, 3.11, 3.06] Second Model: S = [0.7561, 0.1871, 0.0112, 0.0456] n RMSE( Ŝm) RMSE( ŜSP F ) 10 2 10 -2 [6.14, 3.72, 1.22, 2.01] 10 -1 [1.07, 1.00, 1.01, 0.99] 10 3 10 -2 [1.97, 1.17, 0.33, 0.60] 10 -2 [3.36, 3.16, 3.14, 3.13]

		Second Model: S = [0.7561, 0.1871, 0.0112, 0.0456]
	n	Bias( Ŝm)	Bias( ŜSP F )
	10 2	10 -3 [5.72, 0.81, 0.17, 0.18]	10 -3 [1.01, 5.64, 5.68, 05.44]
	10 3	10 -4 [6.27, 1.05, 0.08, 0.61]	10 -4 [5.38, 2.92, 7.14, 7.56]
		First Model: S = [0.5438, 0.2639, 0.0635, 0.1288]
	n	Std( Ŝm)	Std( ŜSP F )
	10 2	10 -2

Acknowledgements This work has been partially supported by the French National Research Agency (ANR) through COSINUS program (project COSTA-BRAVA no ANR-09-COSI-015).

1. First Model: let (γ 1 l ) l and (γ 2 l ) l be defined as γ 1 l = l -1/2-1/100 , 1 ≤ l ≤ L, γ 1 l = 0, l > L; γ 2 l = l -1 , 1 ≤ l ≤ L, γ 2 l = 0, l > L. 2. Second Model: we suppress the two first coefficients of γ defined as in the first model. Let (γ 1 l ) l and (γ 2 l ) l be defined by γ 1 l = γ 1 l+2 and γ 2 l = γ 2 l+2 . 3. Third Model: we replace the fourth first coefficients of γ defined as in the first model by 0. Let (γ 1 l ) l and (γ 2 l ) l be defined by γ 1 i = γ 2 i = 0 for i = 1, . . . 4 and for l ≥ 5 by γ 1 l = γ 1 l and γ 2 l = γ 2 l .

We denote by S = (S 1 , S 2 ) the vector of Sobol indices. We perform N sim = 5000 simulations, we set L = 100 and we study the influence of the parameter n with n = 10 2 or 10 3 . We compare the estimator Ŝm of the vector S defined in Section 3.1 with the SPF estimator defined in Section 3.2. Both estimators are based on 3n observations. The observations are i.i.d. for the estimator Ŝm and obey to the design described in Section 3.2 for ŜSP F . We set m = 3n/ log 3n (i.e. m = 7 or 19) in the definition of the estimator Ŝm .

In the following tabulars we report an empirical estimator of the bias (Bias( Ŝm ) and Bias( ŜSP F )), of the standard deviation (Std( Ŝm ) and Std( ŜSP F )) and of the square root of the Mean Squared Error (RM SE( Ŝm ) and RM SE( ŜSP F )) for each estimator. Those simulations show that as expected the method of the paper provides better MSE as soon as the signal is concentrated on the first terms of the basis which is considered. If it is not the case as in the third model, the SPF can be more efficient. Models 2 and 3 are unfavorable to our estimator, but it is worth to notice that in many cases, considering the linear functional regression model, when a signal has weak components on the first elements of the karhunen-Loève basis it tends to make the influence on the model (Sobol index) smaller.

Using Karhunen-Loève expansion of the fractional Brownian motion

We consider the model ( 18) with p = 2 and ε i = 0 for all i. The processes X 1 and X 2 are defined by the coefficients of their expansion onto the Karhunen-Loève basis of the fractional Brownian motion (we use an approximation of the eigenvalues) :

Note that no close form of the Karhunen-Loève basis of the fractional Brownian motion is known but we don't need it for the simulations, we simply use the fact that X i , ϕ l = √ λ l ξ i l , where the variables (ξ i l ) l≥1 , i = 1, 2 are i.i.d. standard normal variables. Concerning the functions β 1 , β 2 , they are respectively caracterized by the coefficients of their expansion onto the same basis:

1. First Model: let (γ 1 l ) l and (γ 2 l ) l be defined as

2. Second Model: we suppress the two first coefficients of γ defined as in the first model. Let (γ 1 l ) l and (γ 2 l ) l be defined by γ 1 l = γ 1 l+2 and γ 2 l = γ 2 l+2 .

We perform N sim = 5000 simulations, L = 100 and study the influence of the parameter n. We set H = 1/8, m = 3n/ log(3n) in the definition of Ŝm and denote by S = (S 1 , S 2 ) the vector of Sobol indices.

In the following tabulars we report the simulation results. The same conclusion holds: if the signal is not concentrated on the first terms of the basis which is considered, then the SPF can be more efficient. When this is not the case, our estimator has better performances than the SPF for a similar number of observations.

Second example

Using Karhunen-Loève expansion of the Brownian motion

We consider the model ( 18) with p = 4 and ε i = 0 for all i. The processes X k i for 1 ≤ k ≤ p and 1 ≤ i ≤ n are i.i.d. and defined as X 1 and X 2 in Section 4.1.1.

1. First Model: concerning the functions β k for k = 1, . . . , 4, they are respectively caracterized by their coefficients (γ k l ) l 's k = 1, . . . , 4 of their expansion onto the same basis:

2. Second Model: we replace the first coefficient of γ defined as in the first model by 0.

Let for k = 1, . . . , 4, (γ k l ) l be defined by γ k 1 = 0 and for l ≥ 2 by γ k l = γ k l .

3. Third Model: we suppress the two first coefficients of γ defined as in the first model. Let for k = 1, . . . , 4, (γ k l ) l be defined by γ k l = γ k l+2 .

We perform N sim = 5000 simulations, L = 100 and we study the influence of the parameter n, where 5n observations are used for both methods. We set m = √ 5n/ log(5n) in the definition of Ŝm and denote by S = (S 1 , S 2 , S 3 , S 4 ) the vector of Sobol indices.

In the following tabulars we report the simulation results. We see here that the method based on U -statistics always leads to better results than the SPF method. This is due to the fact that the SPF method needs 5n simulations to estimate the four Sobol indices. With the same number of simulations, the method proposed in the paper performs better, even in the second and third model. Let us recall that a drawback of the method proposed in this paper is that the distribution of the input processes X 1 , . . . , X p is assumed to be known which is not necessary for the method SPF.


Philippe Delaye
email: philippe.delaye@iota.u-psud.fr

Andrei Fotiadi

Gérald Roosen

RECIPROCITY AND PHASE CONJUGATION : APPLICATION TO THE STUDY OF THE DOUBLE PHASE CONJUGATE MIRROR

We develop a theoretical model that shows that the overlap integral of two counterpropagating waves is an invariant along propagation in a reciprocal optical system. We use this property to implement a method that allows a quantitative measurement of the conjugation fidelity of a phase conjugate mirror. Applied to the double phase conjugate mirror (DPCM) this invariant can be rewritten as an equality between the ratio of the conjugation fidelity of the two ports of the DPCM and the ratio of the phase conjugate transmission of the two beams. If restricted to a plane wave case this relation becomes the well known equality of the diffraction efficiencies in both directions.

We use the model to implement an experimental characterization of the DPCM, measuring its conjugation fidelity, and demonstrating that it is a reciprocal optical system.

I. INTRODUCTION

A phase conjugate mirror is a device that creates a wavefront that is the phase conjugate of an incident wavefront, i.e. it is a wave that has exactly the same structure than the incident wave but that is counterpropagating to it. This property is very important as it allows to compensate distortions encountered by an optical wave through propagation in an optical system. The phase conjugate wave follows exactly the same path as the incident wave and sees exactly the same distortions but with opposite sign and what compensates them. This is, for example, efficiently used in laser systems to generate high peak power high quality laser beams [START_REF]Optical phase conjugation[END_REF]. Due to this important property, a lot of studies were devoted to the problem of the reconstruction of a wavefront with phase conjugation from a theoretical point of view as well as experimentally [START_REF] Wolf | Phase conjugacy and symmetries in spatially bandlimited wavefields containing no evanescent components[END_REF][START_REF] Nazarathy | Phase conjugacy and symmetries in general optical systems[END_REF][START_REF] Beckwith | Image distortion in multimode fibers and restoration by polarization preserving phase conjugation[END_REF][START_REF] Tomita | Theory of polarization and spatial information recovery by modal dispersal and phase conjugation[END_REF]. Nevertheless most of these works supposed that the propagation media was lossless and that the phase conjugate mirrors were perfect, what is a particular case far from reality. It is only recently [START_REF] Gu | Partial phase conjugation, fidelity, and reciprocity[END_REF][START_REF] Campbell | Fidelity of image restoration by partial phase conjugation through multimode fibers[END_REF][START_REF] Mcmichael | Efficiency of phase conjugation for highly scattered light[END_REF], that the problem of partial phase conjugation (i.e. the case of a phase conjugate mirror that conjugates only a part of the incident wave) had been tackled.

These works are the starting points in our characterization of the Double Phase Conjugate Mirror (DPCM) [START_REF] Weiss | Double phase-conjugate mirror: analysis, demonstration, and applications[END_REF][START_REF] Zozulya | Double phase conjugate mirror: convection and diffraction[END_REF][START_REF] Orlov | Conjugation fidelity and reflectivity in photorefractive double phase-conjugate mirror[END_REF]. Our aim is to characterize the quality of this DPCM, and especially to measure quantitatively a very important parameter, the conjugation fidelity, that measures the part of the returned wave that is exactly the phase conjugate of the incident wave. In a first section of the paper we develop a theoretical model that gives the background of the study. It allows to define a set-up to quantitatively measure the conjugation fidelity, and to show that if the DPCM is a reciprocal device a relation links the ratio of conjugation fidelities on both beams to the ratio of the phase conjugate transmissions. In a second part, we experimentally measure the conjugation fidelity of the DPCM and demonstrate its reciprocity.

II. THEORETICAL STUDY

We use the Butterwerk's formalism [START_REF] Butterweck | General theory of linear, coherent optical data-processing systems[END_REF] that defines an optical system by a transfer function with specific characteristics depending on the properties of the optical system (reciprocity, losslessness, symmetry,...). However, Butterweck's formalism [START_REF] Butterweck | General theory of linear, coherent optical data-processing systems[END_REF][START_REF] Butterweck | Principles of optical data-processing[END_REF] is not exactly adapted to our problem as it considers only scalar field. That is why we use a derivation of this formalism developed to take care of vectorial fields [START_REF] Van Weert | A vectorial treatment of optical systems[END_REF][START_REF] Van Weert | Losslessness of optical systems in a vectorial treatment[END_REF]. This formalism will allow to define invariant quantities along propagation, that in the following will permit to describe the properties of our setup. This reciprocity invariant will then be used to describe a method to measure the conjugation fidelity of a phase conjugate mirror, and a method to show that the DPCM is a reciprocal element

II.1. Position of the problem

The basic scheme that we will use (Fig. 1) is an optical system that transfers optical fields between two plane 1 and 2 (with position in the planes given by variables r 1 and r 2 ). All the fields are considered monochromatic with pulsation ω. Field r E r,z

( )e i !t" kz ( ) propagates along r Z from plane 1 towards plane 2, and field r F r, z

( )e i !t + kz ( ) propagates from plane 2 towards plane 1. The fields are vectorial quantities, i.e. they are polarized. We only consider the case of paraxial waves (i.e. the highest spatial frequency is small compared to the wavenumber k=2π/λ) [START_REF] Butterweck | Principles of optical data-processing[END_REF][START_REF] Van Weert | A vectorial treatment of optical systems[END_REF], so the field vectors completely lie in the plane which is perpendicular to r Z . For each plane we define two perpendicular polarization axes r X and r Y . The complex amplitude r E and r F of the fields are written :

r E = E X r X + E Y r Y Or in matricial notation : r E = E X E Y ! " # $ %
where E X and E Y are complex quantities.

For every vector r E or r F we can also determine the only associated scalar value, which characterises the average intensity of the field I z

( ) = t r E r, z ( ) r E * r, z ( )dr

!

(integration is done on the whole plane).

According to the Butterweck's notations [START_REF] Butterweck | General theory of linear, coherent optical data-processing systems[END_REF], adapted to the vectorial case by Van Weert [START_REF] Van Weert | A vectorial treatment of optical systems[END_REF], the transfer from plane 1 towards 2 occurs according to the following law :

r E 2 r 2 ( ) = T 2 1 r 2 ,r 1 ( ) r E 1 r 1 ( )dr 1 ! (1)
With the transfer matrix :

T 2 1 r 2 ,r 1 ( ) = T 2 1XX r 2 ,r 1
( ) T 2 1XY r 2 ,r 1 ( ) T 2 1YX r 2 ,r 1
( ) T 2 1YY r 2 ,r 1 ( ) ! " # $ % &
similarly the transfer from 2 to 1 occurs according to the law:

r F 1 r 1 ( ) = T 1 2 r 1 , r 2 ( ) r F 2 r 2 ( )dr 2 ! (2)
We also define the inverse transfer as (the propagation is still from plane 1 towards plane 2) :

r E 1 r 1 ( ) = T 2 1 !1 r 1 ,r 2 ( ) r E 2 r 2 ( )dr 2 " (3)
The transfer matrix T 2 1 r 2 ,r 1 ( ) and its inverse T 2 1

!1 r 1 ,r 2
( ) are linked by :

T 2 1 r 2 ,r 1
( )T 2 1 !1 " r 1 , r 2 ( )dr 2 # = $ r 1 ! " r 1 ( )I (4)
I being the identity matrix

I = 1 0 0 1 ! " # $ % and ! r 1 " # r 1 (
) the Dirac distribution.

The general system we now consider has some specific properties, among which the most important and the most used are reciprocity and/or losslessness. In the following we will consider that the only property, that the studied systems possess, is reciprocity. The case of a lossless system (very used theoretically but generally not obtainable in practice) will be considered in the appendix. The reciprocity (using the fact that we consider only paraxial fields) can be translated into relations between transfer matrices [START_REF] Van Weert | A vectorial treatment of optical systems[END_REF][START_REF] Van Weert | Losslessness of optical systems in a vectorial treatment[END_REF] :

T 1 2 r 1 ,r 2 ( )= t T 2 1 r 2 , r 1 ( ) (5
) t T being the transpose of matrix T .

II.2. Reciprocity invariant

Equation ( 5) is our basic relation. It allows to define a quantity that is invariant along the propagation of the fields. We first calculate the overlap integral (i.e. the scalar product) of the two counterpropagating fields at the output of the system :

t r E 2 r 2 ( ) r F 2 r 2 ( ) dr 2 ! = t T 2 1 r 2 ,r 1
( ) t r E 1 r 1 ( )dr 1 ! ( ) r F 2 r 2 ( )dr 2 = t r E 1 r 1 ( ) t T 2 1 r 2 , r 1 ( ) r F 2 r 2 ( )dr 2 ! ( ) dr 1 ! ! (6)
We now use expression (5) :

t r E 2 r 2 ( ) r F 2 r 2 ( ) dr 2 ! = t r E 1 r 1 ( ) T 1 2 r 1 , r 2 ( ) r F 2 r 2 ( ) dr 2 ! ( ) dr 1 ! (7)
what gives (with definition (2)) the expected results, i.e. that the overlap integral of two waves is conserved through propagation in a reciprocal system :

t r E 1 r 1 ( ) r F 1 r 1 ( )dr 1 ! = t r E 2 r 2 ( ) r F 2 r 2 ( )dr 2 ! (8)
This invariant relation is important as it says that the scalar product of two counterpropagating waves is identical in every point of the propagation path in a reciprocal optical system. In particular, two orthogonal waves (i.e. such that t r E 1 r 1 ( )

r F 1 r 1 ( )dr 1 ! = 0
) in a plane of the propagation, remain orthogonal everywhere.

If now, we consider any plane parametrized by an abscissa z, the reciprocity invariant can be written as :

t r E r, z ( )

r F r, z ( )dr ! = c z (9)
c z being a constant independent of z, r E and r F being two fields that propagate in the optical system in the + r Z and -r Z directions respectively. This equation is the basic relation that we will use in the following, firstly to describe a method that can be used to quantitatively measure the conjugation fidelity of a phase conjugate mirror, secondly to show that the double phase conjugate mirror is reciprocal.

II.3. Application to the measurement of the phase conjugation fidelity

II.3.a. Reciprocity and phase conjugation

Relation ( 9) is verified along propagation as soon as the optical system is reciprocal, which is the case for most optical systems. In particular a multimode fiber is reciprocal, in the limit it is stable and its characteristics do not fluctuate temporally [START_REF] Lefèvre | Gyromètre interférométrique à fibre optique[END_REF]. This stationary condition can be generalized to a quasi steady state condition as soon as the time of variation of the characteristics is slow compared to the transit time of the light through the fiber, what is usually the case.

Let us consider the returning field r F (Fig. 1). It is written, at each z position, as a linear combination of two orthogonal fields (not necessarily normalized) :

r F r, z ( ) = a z ( ) r E * r, z ( ) + r F ! r, z ( ) (10)
The first term is the phase conjugate of the incident field r E , at the same position z, multiplied by a constant a(z) (a constant regarding r but eventually depending on position z). The second one is a field r F ! orthogonal to r E , i.e. that verifies :

t r E r, z ( ) r F ! r, z ( )dr " = 0 (11)
Using these definitions we can show that :

t r E r, z ( ) r F r, z ( )dr !
= a z ( )I z ( ) = c z (12)
Writing the field according to relation [START_REF] Zozulya | Double phase conjugate mirror: convection and diffraction[END_REF] shows that if we send the field r F ! r, z ( ) in the optical system, this field will never create a phase conjugate component of the incident field, as there is conservation of the quantity t r E r, z ( )

r F r, z ( )dr
! (if at a given position the reciprocity invariant between the two fields is null, it is null everywhere). If a field component strictly non phase conjugated can never reconstruct a phase conjugate beam, the inverse proposition is generally wrong in presence of losses. To illustrate this point we take the example of a wave passing through an aperture in a screen (Fig. 2). Even if in the upper part of the screen we create a perfectly phase conjugated wave and if this wave is transmitted without any loss through the aperture, it is no more a perfectly conjugated wave in the lower part of the screen since the part of the field blocked by the screen is not reconstructed.

II.3.b. Transmission of a phase conjugate wave through an optical system

The last exposed point will be confirmed by looking at the losses that are seen by the phase conjugate component. Going back to expression (10) that gives r F , we now calculate the intensity of the phase conjugate component. It writes as :

I PC z ( ) = t a z ( ) r E * r,z ( ) ( ) a z ( ) r E * r, z ( ) ( ) * dr ! = a z ( ) 2 r E r, z ( ) 2 dr ! = a z ( ) 2 I z ( ) (13)
According to this, the transmission of the phase conjugate component a z ( )

r E * r, z
( ) between a plane parametrized by the abscissa z 2 and a plane of abscissa z 1 , is defined as :

T PC z 1 ,z 2 ( ) = I PC z 1 ( ) I PC z 2 ( ) = a z 1 ( ) 2 I z 1 ( ) a z 2 ( ) 2 I z 2 ( ) (14)
Using relation [START_REF] Butterweck | General theory of linear, coherent optical data-processing systems[END_REF], we then show that the transmission of the phase conjugate component of the returned wave r F is strictly equal to T z 2 ,z 1 ( ) the transmission of the incident wave r E between the plane with abscissa z 1 and the one with abscissa z 2 :

T PC z 1 ,z 2 ( ) = I z 2 ( ) I z 1 ( ) = T z 2 ,z 1 ( ) (15)
This relation is true whatever the nature of the loss that causes the decrease of the signal. This can be either an absorption loss, a masking by a screen or the loss induced by the necessity to repolarize a depolarized beam.

II.3.c. Measure of the conjugation fidelity

We now calculate the intensity of the returning field r F , using expression [START_REF] Zozulya | Double phase conjugate mirror: convection and diffraction[END_REF] and the fact that r F ! is orthogonal to r E * . We thus have :

I R z ( ) = t r F r, z ( ) r F * r, z ( )dr ! = a z ( ) 2 I z ( ) + I " z ( ) (16)
with

I ! z ( ) = r F ! r, z ( ) 2 dr "
, the intensity of the non phase conjugate component of the wave.

We now come to the main point of our study, the quality of the phase conjugation that is characterized by the conjugation fidelity [START_REF] Zozulya | Double phase conjugate mirror: convection and diffraction[END_REF]:

CF z ( ) = t r E r,z ( ) r F r, z ( )dr ! 2 r E r,z ( ) 2 dr ! r F r,z ( ) 2 dr ! (17)
what can be rewritten taking into account what precedes (( 16) and ( 12)) as :

CF z ( ) = a z ( ) 2 I z ( ) a z ( ) 2 I z ( ) + I ! z ( ) = I PC z
( )

I PC z ( ) + I ! z ( ) (18)
The transmission of the returned field r F between a plane of abscissa z 2 and a plane of abscissa z 1 will then be, using relations ( 18) and ( 14) :

T R z 1 , z 2 ( ) = I R z 1 ( ) I R z 2 ( ) = CF z 2 ( ) CF z 1 ( ) T PC z 1 ,z 2 ( ) (19)
As the transmission of the phase conjugate component equals the transmission of the initial beam (15), we then have :

T R z 1 , z 2 ( ) = CF z 2 ( ) CF z 1 ( ) T z 2 ,z 1 ( ) (20)
This relation is the main expression that we will now use in our experimental study as it is more practical than the reciprocity invariant. It links the conjugation fidelity in two planes to the transmission of the waves in both directions between these two planes. The problem will now be to find a way to separate the phase conjugate component from the non phase conjugate one. A method using a measurement of the quality of images reconstructed by phase conjugation was previously proposed [START_REF] Gu | Partial phase conjugation, fidelity, and reciprocity[END_REF][START_REF] Campbell | Fidelity of image restoration by partial phase conjugation through multimode fibers[END_REF][START_REF] Orlov | Conjugation fidelity and reflectivity in photorefractive double phase-conjugate mirror[END_REF]. We prefer to use a single mode filter (made of a polarization maintaining single mode fiber and a polarizer) for this selection. In such a filter one mode and only one can propagate. Thus both counterpropagating waves will be strictly phase conjugated to each other. Moreover the measurement of the image quality does not strictly give the conjugation fidelity as defined here. In one case [START_REF] Gu | Partial phase conjugation, fidelity, and reciprocity[END_REF][START_REF] Campbell | Fidelity of image restoration by partial phase conjugation through multimode fibers[END_REF], the measured image to background ratio (IBR), is a parameter that should goes to infinity for a perfect phase conjugation (no background noise present) and not to a maximum of one as in our case. In the other case [START_REF] Orlov | Conjugation fidelity and reflectivity in photorefractive double phase-conjugate mirror[END_REF], it is said that the fidelity measurement does not take into account the phase of the returned beam, what is done in our case. Now we consider that we filter the phase conjugate component in the plane of abscissa z 1 , i.e. that we assure that in this plane the intensity of the non phase conjugate component is zero (I ! z ( ) = 0 ) and thus the conjugation fidelity is CF(z 1 )=1 (according to (18)). We thus obtain a relation that gives a measurement of the conjugation fidelity in the plane with abscissa z 2 :

CF z 2 ( ) = T R z 1 ,z 2 ( ) T z 2 ,z 1 ( ) (21)
The conjugation fidelity is in general a parameter very hard to evaluate quantitatively. We here propose a very simple way to obtain it through a simple measurement of the transmission of two beams.

The experimental set-up to evaluate the conjugation fidelity of a phase conjugate mirror is thus the following (Fig. 3). A laser beam is filtered by a single mode filter and then sent through a reciprocal optical system (i.e. in practice any classical optical system) on the phase conjugate mirror. We then measure the transmission of the incident beam through the single mode filter and the optical system, as well as the transmission of the beam returned by the phase conjugate mirror. The ratio of both transmissions directly gives the fidelity of the phase conjugate mirror, according to (20).

II.4. Reciprocity of the double phase conjugate mirror

Until now we have used the model to describe how a phase conjugate wave is transmitted through a reciprocal optical system and how we can use this to measure the conjugation fidelity of the phase conjugate mirror. We will now show how the same expression can be used to prove that an optical system is reciprocal and apply this to an important kind of phase conjugate mirror : the double phase conjugate mirror (DPCM) [START_REF] Weiss | Double phase-conjugate mirror: analysis, demonstration, and applications[END_REF][START_REF] Zozulya | Double phase conjugate mirror: convection and diffraction[END_REF][START_REF] Orlov | Conjugation fidelity and reflectivity in photorefractive double phase-conjugate mirror[END_REF]. Until now the question of its reciprocity has only be considered once and only theoretically [START_REF] Gu | Reciprocity in photorefractive wave mixing[END_REF]. It was shown that the equality of the diffraction efficiency in both directions (in a plane wave description of the DPCM) is a direct consequence of the reciprocity of the DPCM. Our model is an extension to the more general case of non plane incident waves and will give us an experimental method to show that the DPCM is reciprocal. II.4.a. Conservation of the fidelity-transmission product We consider the DPCM as our optical system (Fig. 4). The DPCM is a phase conjugate mirror in which two beams 1 and 2 that can be mutually incoherent are incident on a photorefractive crystal. The DPCM creates two beams 3 and 4 that have the phase conjugate structure of beams 1 and 2 respectively. Moreover the DPCM is a phase conjugate mirror that operates in transmission as beam 4 is issued from the diffraction of beam 1 on the grating written in the photorefractive crystal (with the same properties for beams 3 and 2). We can thus make the analogy between the different notations of Fig. 4 and Fig. 1, r E 1 is the field that characterizes wave 1, r E 2 characterizes wave 4, whereas r F 1 and r F 2 characterize waves 3 and 2 respectively. Now if the DPCM is reciprocal, relation (20) is verified between the two planes of abscissa z 1 and z 2 . We thus have :

T R z 1 , z 2 ( ) = CF z 2 ( ) CF z 1 ( ) T z 2 ,z 1 ( )
This can be rewritten as a fidelity-transmission product conservation :

CF cw T DPCM cw = CF ccw T DPCM ccw ( 22
)
where CF is the conjugation fidelity of the two waves issued from the DPCM and T DPCM is the phase conjugate transmission of the waves that go through the DPCM (i.e. ratio of the intensity of the beam that goes in the direction of the phase conjugate beam I 3 (resp. I 4 ) on the intensity of the incident beam I 2 (resp. I 1 )). Indices cw and ccw indicate that the waves propagate in the clockwise and counterclockwise directions respectively. Thus CF cw is the conjugation fidelity of wave 4 (compared to the incident wave 2) whereas T DPCMcw is the ratio of intensity of wave 4 to the intensity of wave 1. Thus if the DPCM is reciprocal the product of the conjugation fidelity by the transmission of the phase conjugate wave is the same in both propagation directions. This expression generalizes the relation derived in Ref. [START_REF] Gu | Reciprocity in photorefractive wave mixing[END_REF] to the case of waves that are not plane waves. Indeed, in this plane wave case, as two counterpropagating plane waves are necessarily phase conjugate to each others, the conjugation fidelity is necessarily 1 and relation ( 22) becomes the equality between the diffraction efficiency of the grating in both directions as stated in Ref. [START_REF] Gu | Reciprocity in photorefractive wave mixing[END_REF]. This is no more true for non plane waves as illustrated in the case of a DPCM operated with two beams having different diameters (Fig. 5). The interaction area is limited to the area of the crystal covered by both beams, which in first approximation is the area covered by the smallest beam. The small beam is totally covered by the large beam, it thus interacts with the written grating on its entire aperture, and thus diffracts with a large efficiency. On the contrary a large portion of the large beam does not participate to the grating and is thus not diffracted giving a small diffraction efficiency. We now look at the quality of the phase conjugation. On the side of the large beam, the structure of the diffracted wave corresponds to a small portion of the beam only and thus corresponds to a very bad quality of the phase conjugation, with the reverse phenomenon at the other side of the crystal. We see that on one side we have a good energetic efficiency, together with a bad quality of the phase conjugation, and, on the other side, a bad efficiency compensated by a very good quality of the phase conjugation. This example illustrates in a simple manner how relation ( 22) allows to characterize the DPCM, with the help of a new invariant that works whatever the structure of the incident waves and that thus works in all the operating cases of the DPCM.

II.4.b. Measure of the reciprocity of the DPCM The previous invariant we derived, is a simple consequence of the reciprocity of the DPCM. We now have to prove that the DPCM is a reciprocal optical system. Once again we will use expression (20), but applied to a more complicated optical system (Fig. 6). Plane 1 (at abscissa z 1 ) is a single mode filter (constituted of a polarization maintaining single mode fiber and a polarizer). Then the beam is transferred from the single mode filter to the DPCM through a multimode fiber. This multimode fiber send on the DPCM a speckled beam that allows to prevent the conical diffraction that can occur with plane waves and that is detrimental to the conjugation quality of the DPCM. The optical set-up is then duplicated for the other branch of the set-up until plane 2 (at abscissa z 2 ). The optical system is thus the whole set-up (Fig. 6) :

" single mode filter -> multimode fiber -> DPCM -> multimode fiber-> single mode filter ". With this optical system we have perfect single mode filters in the planes of abscissa z 1 and z 2 , thus we have for both conjugation fidelities CF(z 1 )=CF(z 2 )=1. Expression (20) then becomes :

T R z 1 , z 2 ( ) = T z 2 ,z 1 ( ) (23)
This relation tells us that if the system is reciprocal (i.e. if each of its elements are reciprocal [START_REF] Butterweck | Principles of optical data-processing[END_REF]) the transmission coefficients in both propagation directions measured in place where CF=1 are strictly identical whatever the exact value of these transmissions. We have thus a very simple way to verify that our DPCM is a reciprocal device (since all the other elements are considered reciprocal), and this whatever the operating regime of the DPCM.

III. EXPERIMENTAL VERIFICATION

We have developed a theoretical model that gives a way to characterize phase conjugation fidelity quantitatively and to control experimentally whether the DPCM is a reciprocal element. This is expressed by a relation that links phase conjugate fidelity on both sides to diffraction efficiency in both directions. We now present the experimental study, stressing the precautions that have absolutely to be taken in order to have a meaningful and reliable proof of the reciprocity of the DPCM and of the measurement of the conjugation fidelity.

III.1. Experimental set-up

Our aim here is to characterize a double phase conjugate mirror in a certain configuration and not to optimize its performances. We choose a configuration of the photorefractive BaTiO 3 crystal that works quite well and that easily gives high efficiencies [START_REF] Ph | Heterodyne detection of ultrasound from rough surfaces using a double phase conjugate mirror[END_REF]. It is shown in Fig. 7. Typical phase conjugate transmission was around 40%, which taking into account the ≈50% loss in the crystal (absence of antireflection coatings, and absorption of the crystal) gives a corrected diffraction efficiency around 80% at the level of the crystal.

The experimental set-up (Fig. 8) is constructed around the DPCM with the architecture proposed above (Fig. 6). The laser beam is issued from an Argon laser emitting at λ=514nm in a single line highly multimode regime, in order to have a short coherence length. It is sent on a non polarizing beam splitter BS, in order to form the two pump beams of the DPCM. Each beam is sent on a polarization maintaining single mode fiber (PMSF), after passing a polarizer (P), in order to form our single mode filter. The half wave plate (λ/2) aligns the polarization along the proper polarization axis of the fiber. Just in front of the polarizers we place 50/50 beam splitters (BS 1 and BS 2 ) to extract the returning phase conjugate signal. At the output of the single mode fibers the beam is sent into a multimode fiber (MF) (core diameter 50µm). The multimode fibers have different lengths, in order to assure that both pump beams of the DPCM are mutually incoherent. This prevents the writing of reflection gratings between the pump beams. At the output of the multimode fibers the beams are polarized with a polarization beam splitter (PBS), in order to have extraordinary polarized pump beams incident on the photorefractive crystal (PRC). A thick glass plate (GP) in one of the pump beam extracts a part of the phase conjugate beam to follow the growth of the diffracted signal in the vicinity of the photorefractive crystal. The incident power of each pump beam on the DPCM is typically between 200µW and 2mW depending on the experiment. Three detectors (D 1 , D 2 and D 3 ) allow to follow the phase conjugate signal after the single mode filter for the first two and at the level of the DPCM for the third one (the diameter of D 3 is adjusted to the diameter of the pump beam at this position, i.e. to the diameter of the expected phase conjugate beam).

III.2. Conjugation fidelity measurement

In order to measure the conjugation fidelity, we use our set-up (Fig. 8) that was made according to the principle proposed in Fig. 3. According to this and using expression (21) a simple measurement of the different intensities defined in Fig. 3, allows to evaluate the conjugation fidelity of the phase conjugate mirror. We perform this measurement and find values of the conjugation fidelity that were much greater than 1. This result indicates that there is a problem either in the model or in the measurement procedure, as CF is evidently smaller or equal to 1.

In fact there is a problem in both. Indeed, using expression (21) we suppose that at the position of the measurement, i.e. at the entrance of the single mode filter, we have CF(z 1 )=1. This is certainly true inside the single mode fiber, but no more true if we consider that the measurement is performed outside the fiber. The coupling of the incident beam inside the single mode fiber is certainly not perfect, there are losses. Thus, as we saw previously, the phase conjugate component of the wave will see the same losses, that will give a decrease of the phase conjugate fidelity, just as in the case of masking by a screen we discussed previously (Fig. 2). So we can not use expression (21) in our evaluation. We have to go back to relation (20), with the problem that CF(z 2 ) is no more evaluated directly. We have first to measure CF(z 1 ) before having a value of CF(z 2 ). Another way would be to assure that the coupling inside the single mode fiber is without any loss, but we think that it is hardly possible in practice.

In order to calculate CF(z 1 ), we use the same expression (20) but exchanging the role of both planes. We replace the phase conjugate mirror to be evaluated by a system that is a perfect phase conjugator. In such a case we have CF(z 2 )=1 and from the measurement of the intensities defined in Fig. 3, we deduce the different transmissions that give CF(z 1 ). It is equivalent to consider that we make a relative measurement in which we compare our phase conjugate mirror to a perfect phase conjugate mirror. This method is less accurate than the absolute measurement we expected to perform, but it nevertheless gives a good idea of the quality of the phase conjugation fidelity. The precision of the measurement will depend on the quality of the "perfect" phase conjugate reference we will choose. To choose the "perfect" phase conjugate mirror we remind us that at the output of the single mode fiber, we have a good quality beam, which is collimated, and can thus be considered as a plane wave. When sent on a plane mirror along its normal, this plane wave generates another counterpropagating plane wave that is its exact phase conjugate (the non unity reflection coefficient of the mirror just reduces the efficiency of the "perfect" phase conjugate mirror but not the conjugation fidelity that stays equal to 1). So we experimentally replace the multimode fiber and the phase conjugate mirror by the plane mirror perpendicular to the collimated beam. A very accurate alignment both for the collimation of the beam at the output of the single mode fiber and the autocollimation of the plane mirror has been achieved. As we couple back the reflected signal in the single mode fiber, we are sure that only the phase conjugate component arrives at the entrance of the single mode fiber, and that the optimization of I out (Fig. 3) is a characteristics of the quality of our "perfect" phase conjugate mirror. Once the returned signal optimized, we perform the measurement of CF(z 1 ), using relation (20) with CF(z 2 )=1. We deduce CF(z 1 )=0.35, what is far from the supposed value of 1 we were assuming at the beginning.

We can now evaluate the conjugation fidelity of the DPCM, using the proposed set-up, with relation (20) and the now known value of CF(z 1 ). The measurement of I out is given by detector D 2 , and I PCMout is given by detector D 3 , both independently calibrated regarding the exact value of the measured signal, whereas I inc =34.2mW and I PCMinc =1.36mW. We are now able to follow the time response of the DPCM, both for the intensity reflectivity R=(I PCMout /I PCMinc ) of the DPCM through D 3 and the conjugation fidelity CF(z 2 )=CF(z 1 )*(I out /I PCMout )*(I inc /I PCMinc ) through D 2 . The curves are presented in Fig. 9. They are taken at the build-up of the DPCM and after some time in the "steady-state" regime. The curve shows that the reflectivity is stable whereas the conjugation fidelity strongly oscillates, reaching a value close to 1, at some instants. According to these results we conclude that the DPCM diffracts a beam that has the global structure of the incident beam (mainly same diameter), and that looks like a phase conjugate speckle beam if regarded near the DPCM (Fig. 10). The conjugate beam thus "seems" a real conjugate but the measurement of the conjugation fidelity reveals that in fact it is not. As we conjugate a speckle beam we can conclude that the large structure of the beam is conjugated but that its fine speckled structure fluctuates slowly in time, being at some moment perfectly non phase conjugated (CF=0 at t=0 in Fig. 9). This shows that the usual measurement of the conjugation reflectivity is not at all a measurement of the quality of the phase conjugation. Indeed this reflectivity can be important even if the returned beam is not at all phase conjugated. This oscillating operation of the DPCM has been chosen because it showed the decorrelation between the reflectivity value and the quality of the phase conjugation. It is obtained in a specific regime of alignment and intensity ratio of the pump beams incident on the DPCM. We have obtained alignments giving stable operation (to some percents) of the DPCM both for reflectivity and conjugation fidelity for more than one hour.

III.3. The double phase conjugate mirror is reciprocal

We have just seen that the problem of the quality of the coupling of the incident beam in the single mode fiber prevents a direct measurement of the conjugation fidelity. This problem also occurs in the experimental demonstration of the reciprocity of the DPCM, in preventing the direct use of expression (23). Here again we have to go back to relation (20). Nevertheless, the problem occurs to a lesser extent as we do not need to measure exactly both CF(z 1 ) and CF(z 2 ) but only their ratio. As the value of CF(z 1 ) and CF(z 2 ) are fixed and determined by the coupling at the entrance and output of the single mode filter only, they do not depend on the nature of the system between the single mode filters. Then relation (20) will be identical whatever the optical system we place between both single mode filters, as soon as this system is reciprocal. So in order to show that the DPCM is reciprocal, we just have to compare the ratio of the transmissions of the system in both directions for the total system (multimode fiber -DPCM -multimode fiber) and for a reference reciprocal system, in which the output of the single mode fiber is directly coupled in the other fiber by mirrors (insert in Fig. 8). Thus if we can show that we have the same ratio of transmissions for the DPCM and the reference system, we will prove that the DPCM is a reciprocal system.

Experimentally, in our set-up the ratio of the transmissions is directly proportional to the ratio of the signals received by detectors D 1 and D 2 , as the incident signal ratio is constant during the measurement. Moreover, in order to have a better accuracy, we make the transission measurements for several operating points by varying the transmission of the DPCM, i.e. its efficiency (changing the alignment and superposition of the pump beams) and by making the measurement during the build-up of the grating (or by choosing an unstable operation regime as the one presented in Fig. 9). This last point is justified by the fact that the durations involved in the measurements are very long compared to the transit time of the light in the optical system, assuring that the optical system is in a quasi steady state regime. This point is also a posteriori justified by the fact that in all the experiments performed with the DPCM, both signals given by D 1 and D 2 were always exactly proportional to each other, demonstrating that the DPCM is reciprocal whatever its operation regime. For the reference system, we vary the transmission by using a variable diaphragm that cuts the beam. Moreover in order to allow a more direct comparison of both signals a neutral density is used to adjust the transmission of the reference system to the lower transmission value of the DPCM system.

Finally, to increase the accuracy of the measurement at very low transmission of the optical systems (when the returned phase conjugate signal becomes smaller than the returned signal due to the reflection on the entrance angle-cut facet of the single mode fiber), we use shutters (Sh 1 and Sh 2 in Fig. 8) that are alternatively shut. The signal received by each detector during the measurement is shown in Fig. 11. For example detector D 2 , when both shutters are opened receives light that comes from the DPCM and that is the signal we are interested in, but also receives retroreflected light that comes from the different optical elements placed between D 2 and the single mode fiber. So we first close Sh 1 , D 2 does not received light from the DPCM anymore, but it continues to receive the retroreflected light (RR 2 ), which can thus be evaluated. Second we close Sh 2 (Sh 1 being still closed), no more light is sent in the system, we thus measure the zero signal of D 2 (background light and electrical offset). Third, we open Sh 1 , D 2 receives only the light transmitted through the DPCM (without the retroreflected light) what constitutes the real signal we are interested in (PC 2 ). The same kind of signal (but inverted is seen by the other detector D 1 ). This procedure allows to measure exactly the signal transmitted through the DPCM even in the presence of retroreflected light. Nevertheless this procedure has to be applied very carefully, as we have to assure that it does not disturb the measurement. Thus some precautions have to be taken. The main one is to prevent the erasure of the grating when only one beam is present (as well as the recovery of the induced absorption we have observed in our crystal). First we close the shutters for a very short time, typically about 150ms. Second we use a very low pump beam power (typically around 200µW) to decrease at maximum the response time of the DPCM (typically several minutes). All these precautions assure that our measurement procedure is really non destructive and non perturbative.

The result of a measurement is presented in Fig. 12 where we plot the signal received by detector D 2 as a function of the signal received by detector D 1 , for the DPCM and for the reference system (each point is the phase conjugate signal only and is extracted from temporal curves identical to the one shown in Fig. 11). We see that in both cases we have a straight line and that both straight lines have exactly the same slope. This measurement was performed several times with different alignments both of the DPCM and of the reference system, with always the same result. This assures that the DPCM is a reciprocal element and that thus expression ( 22) is true for the DPCM.

IV. CONCLUSION

We have developed a model that allows to describe the transmission of a phase conjugate wave in a reciprocal optical system. Using this model we have shown that the scalar product (i.e. the overlap integral) of two counterpropagating waves is conserved through the propagation in a reciprocal system. Applied to the problem of the propagation of a phase conjugate wave in an optical system, this property allows to show that the phase conjugate component of the returning wave sees exactly the same losses than the incident wave. We also deduce from this model that the ratio of conjugation fidelity in two planes of the propagation equals the ratio of the transmission of the waves between these two planes in both directions. This gives an experimental method to measure quantitatively the phase conjugation fidelity of a phase conjugate mirror, using the filtering by a single mode filter. This method measures the conjugation fidelity of our phase conjugate mirror, compared to the fidelity of a "perfect" phase conjugate mirror used as a reference.

We also use this model to experimentally prove that the DPCM is a reciprocal element. This property allows to find a relation that links the conjugation fidelity on both ports of a DPCM to the phase conjugate transmission in both directions. This relation can be seen as a conservation of the phase conjugate fidelity-transmission product for the two incident beams of the DPCM. It generalizes the equality of the diffraction efficiencies in both directions that was previously proposed using plane wave description of the DPCM.

APPENDIX

For a lossless system the transfer matrix is characterized by :

T 2 1 !1 r 1 ,r 2 ( )= t T 2 1 * r 2 ,r 1 ( ) (A.1)
T * being the conjugate of the matrix T . The invariant that is a consequence of losslessness, expresses classically the conservation of energy through the propagation. It is calculated from expressions (1)-( 4), using relations (A.1) (integration is done on the whole space) :

t r E 1 r 1 ( ) r E 1 * r 1 ( )dr 1 ! = t r E 2 r 2 ( ) r E 2 * r 2 ( )dr 2 ! (A.2)
This relation says that the energy of the field is the same at the entrance and at the exit of a lossless optical system. This relation is well known and commonly used, even if in fact only few systems can be considered as really lossless.

FIGURE CAPTIONS

Figure 1 :

 1 Figure 1 : description of the propagation of a wave through a general optical system.

Figure 2 :

 2 Figure 2 : Illustration of the conjugation quality loss during the transmission of a wave through a mask like an aperture in a screen.

Figure 3 :

 3 Figure 3 : Schematic set-up for the measurement of the conjugation fidelity of a phase conjugate mirror (PCM). The locations where the different intensities I are measured are indicated.

Figure 4 :

 4 Figure 4 : Schematic principle of a double phase conjugate mirror. Beams 1 and 2 are the incident pump beams, whereas 3 and 4 are the conjugated beams.

Figure 5 :

 5 Figure 5 : Illustration of the conservation of the phase conjugate fidelity-transmission product in the DPCM.

Figure 6 :

 6 Figure 6 : Schematic set-up for the demonstration of the reciprocity of the double phase conjugate mirror.

Figure 7 :

 7 Figure 7 : Experimental configuration of the photorefractive crystal used for the operation of the DPCM.

Figure 8 :

 8 Figure 8 : Experimental set-up used for the conjugation fidelity measurement and demonstration of the reciprocity of the DPCM (see text for definition of the different elements). In the insert, the reference optical system used for the reciprocity measurement is shown.

Figure 9 :

 9 Figure 9 : Temporal variation of the reflectivity and of the conjugation fidelity of the DPCM.

Figure 10 :

 10 Figure10: Photographs of a screen intercepting the "phase conjugate" beam transmitted by the DPCM. (A) image of the beam fanning structure when only one pump beam is incident on the DPCM. (B) Image of the phase conjugate beam, after build-up of the grating, the beam fanning collapse in a well defined beam, whose diameter corresponds to the one of the pump beam at the same place. The fine speckled structure is not seen due to saturation of the camera. In both cases the bright spot at right is the transmitted pump beam.

Figure 11 :

 11 Figure 11 : Temporal signal registered by both detectors D 1 and D 2 , when the shutters are closed successively.Figure12: Phase conjugate signal received by detector D 2 as a function of the signal received by D 1 , for the reference system (grey signs and curve) and the DPCM (black signs and curve). Both

Figure 12 :

 12 Phase conjugate signal received by detector D 2 as a function of the signal received by D 1 , for the reference system (grey signs and curve) and the DPCM (black signs and curve). Both slopes are equal within the accuracy of the measurement (3σ error bars).

Figure 1

 1 Figure 1

Figure 6

 6 Figure 5

Figure 12

 12 Figure 7

Acknowledgments : Support from SFIM Industries for this study is greatly acknowledged.