Salah Naili

Vu-Hieu Nguyen

Mai-Ba Vu

Christophe Desceliers

Christian Soize

Ultrasound wave propagation in a stochastic cortical bone plate

Keywords: random media, ultrasound wave propagation, multilayer heterogeneous medium, finite difference scheme, time domain, bone tissue

Ultrasonic guided-wave technologies are powerful nondestructive testing techniques to characterize bone material. This work aims to evaluate the effect due to spatial heterogeneity of bone material properties on its ultrasound response using axial transmission technique. A probabilistic model is introduced to describe the mechanical behavior of bone material. The numerical results focused on studying of FAS (First Arriving Velocity) showing that this quantity strongly depends on the dispersion induced by statistical fluctuations of stochastic elasticity field.

Introduction

It is well-known that cortical bone is a highly complex composite material formed by a hierarchical and multiscale constituents. Due to the fluctuation of pore distribution and physical properties of mineralization of bone tissues, cortical bone at the vascular scale is a heterogeneous and random medium. These factors would not be neglected when performing diagnostics of cortical bone. One of the most usual techniques used for diagnostics of long bone is known as ultrasonic pulsed through-transmission or axial transmission technique (ATT). This technique, which measures the wave velocity in the bone longitudinal direction, has been shown particularly suitable to predict mechanical as well as geometrical characteristics of the bone [START_REF] Lowet | Ultrasound velocity measurements in long bones: measurement method and simulation of ultrasound wave propagation[END_REF]. In the past, most models of ATT considered cortical long bone as a medium with homogeneous properties along its longitudinal direction. Hence, homogeneous or functionally graded material properties have been used to model the cortical bone plate [START_REF] Bossy | Bidirectional axial transmission can improve accuracy and precision of ultrasonic velocity measurement in cortical bone: a validation on test materials[END_REF][START_REF] Haiat | Influence of a gradient of material properties on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF][START_REF] Naili | Influence of viscoelastic and viscous absorption on ultrasonic wave propagation in cortical bone: Application to axial transmission[END_REF][START_REF] Nguyen | Ultrasonic wave propagation in viscoelastic cortical bone plate coupled with fluids: a spectral finite element study[END_REF]. In practice, the exploitation of measured signal data naturally needs to also take into account the uncertainty of material characteristics. However, most of parametric studies of wave propagation in bone are mainly limited to deterministic media. Some studies have recently been carried out to investigate the influence of random properties of cortical bone in the context of ultrasound characterization. A probabilistic model based on the maximum entropy principle has been constructed for considering cortical bone plate as homogeneous or multi-layered media [START_REF] Macocco | Elastoacoustic model with uncertain mechanical properties for ultrasonic wave velocity prediction; application to cortical bone evaluation[END_REF][START_REF] Desceliers | Determination of the random anisotropic elasticity layer using transient wave propagation in a fluidsolid multilayer: model and experiments[END_REF][START_REF] Desceliers | Probabilistic model of the human cortical bone with mechanical alterations in ultrasonic range[END_REF]. It has been shown that a simplified mechanical model with an additional stochastic modeling of bone elasticity properties are able to represent the in vivo measurements in the statistical sense. However, in these works, material inhomogeneities along longitudinal direction of bone have always been neglected.

This paper presents a stochastic model to consider the random fluctuation of material properties in both radial and longitudinal directions in a cortical bone plate. Section 2 presents the description of the problem and the numerical method to compute the time domain solution of the acoustic response of a heterogeneous bone plate coupled with fluid. Next, we provides a procedure to construct a parametric model for the elasticity properties of bone plate. Then, Section 3 shows a numerical test to study the sensitivity of V FAS due to dispersion of the stochastic elasticity field in the bone. Last, Section 4 gives some conclusions of this work. ) is an orthonormal basis. Figure 1 shows a common configuration for modeling an ultrasound axial transmission test. It consists of an anisotropic elastic layer sandwiched between two idealized acoustic fluids. The bone material is assumed to be heterogeneous in the plane (e 1 , e 2 ) but homogeneous along e 3 direction. An acoustic line source, which is parallel to e 3 , produces an excitation at x s = (x s 1 , x s 2 ) inside the upper fluid layer. Hence, the problem may be reduced into a two-dimensional plane strain problem in the (O; e 1 , e 2 ) plane. In Figure 1, the infinite bone layer occupying the domain Ω b with a constant thickness h (Ω b = {x 1 ∈ R, 0 ≥ x 2 ≥ -h}). This bone plate is loaded on its upper and lower surfaces by two fluid halfspaces which represents the soft tissues. The upper fluid domain is denoted by

Ω f 1 (Ω f 1 = {x 1 ∈ R, x 2 ≥ 0}
) and the lower one is denoted by

Ω f 2 (Ω f 2 = {x 1 ∈ R, x 2 ≤ -h}).
The interfaces between the bone (Ω b ) and the fluids (Ω f 1 and Ω f 2 ) are denoted by Γ bf 1 and Γ bf 2 , respectively (see Fig. 1). In that follows, the whole of domain is denoted by

Ω = Ω f 1 ∪ Ω b ∪ Ω f 2 .

Governing equations and numerical methods

For the sake of simplifying, the described trilayer system, which consists of a heterogeneous anisotropic elastic layer Ω b and two homogeneous inviscid acoustic fluid layers, will be all modeled as an elastic medium. At each point x in Ω, the constitutive law equation is expressed by where σ(x, t) and (x, t) denote the stress and strain tensors, respectively; u(x, t) denotes the displacement vector; c(x) is the fourth-order elasticity tensor of the fluid or solid depending on x:

σ = c , with = 1 2 (grad u + (grad u) T ), (1) 6
- e 2 e 1 Fluid 2 (Ω f 2 ) source receivers Fluid 1 (Ω f 1 ) Solid (Ω b ) s s s s s s s s Γ bf 1 Γ bf 2 (x s 1 , x s 2 ) O Γ f∞ 1 Γ f∞
c(x) = c f for x ∈ Ω f 1 ∪ Ω f 2 and c(x) = c(x) for x ∈ Ω b . The elasticity tensor of the fluid c f , which is homogeneous in Ω f 1 ∪ Ω f 2
, is as an isotropic elastic tensor without shear modulus. The elasticity tensor of the solid c depends on x and is determined by 6 elastic constants. Note that in the next section, the probabilistic model only concerns the tensor of elasticity of the domain Ω b .

By neglecting the body force field, the equation of motion in domain Ω reads

div σ = ρ ü, ∀x ∈ Ω, ( 2
)
where ρ(x) is the mass density and div designates the divergence operator. The domain Ω is at rest for t < 0. Note that no interface condition is required for this elastoacoustic problem because both of fluid and bone domains are modeled as elastic media.

An explicit time domain finite difference scheme, which is based on a staggered grid formulation for the velocity and stress components, has been implemented [START_REF] Levander | Fourth-order finite-difference P-SV seismograms[END_REF] to solve the two-dimensional problem (2). The scheme is second-order in time and fourth-order in space. In order to avoid the wave reflected due to the boundaries at finite distances, the PML Perfectly Matched Layer (PM) technique has been is used.

The acoustic response will be captured at a linear array of receivers as shown in Fig. 1. This emitter and receiver configuration is typical one in ultrasonic axial transmission devices for characterizing the cortical layer of bone. When using the axial transmission technique, the earliest event or wavelet (also called by First Arriving Signal, FAS) of the multicomponent signal recorded at the receivers has been the most often investigated because it is considered as a relevant index of bone status [START_REF] Stegman | Cortical ultrasound velocity as an indicator of bone status[END_REF][START_REF] Barkmann | A new method for quantitative ultrasound measurements at multiple skeletal sites[END_REF][START_REF] Hans | Does combining the results from multiple bone sites measured by a new quantitative ultrasound device improve discrimination of hip fracture?[END_REF]. In this work, the quantity of interest to be studied will be the FAS velocity.

Construction of a parametric probabilistic model for an anisotropic elastic medium

Parametric probabilistic models are built by modeling the local physical properties of the medium. In the present study, only the uncertainty of elasticity tensor C(x), of which the mean value is defined by c(x), will be considered. At each point x ∈ Ω b , tensor C depends on 21 independent random variables. In practice, it is impossible to have a sufficient large set of experimental data, especially for bone material, to estimate the probability distribution of the random elasticity tensor. Moreover, the random spatial variation of the elastic property should be considered. To overcome this difficulty, we use Soize's model whereby the probability distribution is built by full-filling the maximum entropy principle. Using this principle, the stochastic elasticity tensor C(x) may be parameterized by its mean values c(x) via its matrix representation, and by a minimal set of essential parameters which consists of only 4 parameters: one scalar dispersion level δ and one vector λ that contains three spatial correlation lengths. The detail of essential steps to estimate the random field C(x; c, δ, λ) at every point x ∈ Ω b may be found in [START_REF] Soize | Random matrix theory for modeling uncertainties in computational mechanics[END_REF][START_REF] Soize | Non-Gaussian positive-definite matrix-valued random fields for elliptic stochastic partial differential operators[END_REF] .

Numerical results

This section will present some illustrative results on studying the ultrasound wave propagation through a random 4mm-thickness cortical bone layer. For this simulation, the pressure source acts on a spatial length of 0.75 mm which is horizontally placed in the fluid at 2 mm from the upper interface of the solid layer. The excitation signal is a Gaussian with a center frequency f c = 1 MHz [START_REF] Desceliers | Three dimensional transient elastic waves in multilayer semi-infinite media solved by a timespace-spectral numerical method[END_REF]. The responses will be calculated at 14 receivers are regularly spaced with a pitch of 0.8 mm and a distance of emitter to closest receiver equal to 11 mm.

Fluid and solid material properties Both fluid layers, which represent the soft tissues and bone marrow, are considered to be an inviscid water with a mass density ρ f = 1 000 kg.m -3 and a bulk modulus K f = 2.25 GPa.

The mean model of the bone material is assumed to be a transversely isotropic medium. This behavior has been experimentally shown by different authors [START_REF] Dong | The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity[END_REF][START_REF] Reilly | The mechanical properties of cortical bone[END_REF][START_REF] Rho | An ultrasonic method for measuring the elastic properties of human tibial cortical and cancellous bone[END_REF] to be a realistic approximation of cortical bone properties. This example uses data given by [START_REF] Dong | The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity[END_REF] who measured the homogenized bone properties by performing tensile and torsional tests with a mechanical testing system on 18 different human femoral bone specimens. The components of the mean elasticity tensor are given by: a 11 = 23.05 GPa, a 22 = a 33 = 15.08 GPa, a 12 = a 13 = 8.72 GPa, a 44 = 3.3 GPa, a 55 = a 66 = 4.7 GPa.

Parameters for the uncertain elasticity model

As discussed before, this study is restricted to only consider the uncertain heterogeneous in the plane (x 1 , x 2 ). Thus only the dispersion δ and two correlation lengths (λ 1 , λ 2 ) need to be introduced to control the random variation of elasticity field. For this example, a correlation length λ 1 = λ 2 = 2×10 -4 m, which may be seen as a typical centerto-center distance between osteons in cortical bone, is used. Two different values of the dispersion δ = 0.1 and δ = 0.3 will be investigated. Figure 2 shows a map illustrating the spatial variation of A ij . Numerical parameters for the finite-difference solver The finite difference simulation is performed for each of 800 Monte-Carlo realizations. A rectangular domain with size 0.03 × 0.01 m is used to be able to contain all details of the model, i.e three layers with the source and receivers. The grid steps in both x 1 -and x 2 -direction are chosen to be identical: ∆x 1 = ∆x 2 = 2 × 10 -5 m to ensure that grid intervals are smaller than about 1/8 of the smallest wavelength in whole domains. As a result, the grid has 1501 × 501 points (about 1.4 × 10 6 degrees of freedom). The time step is chosen by using CFL stability condition ∆t < α∆/c, where α is a constant, ∆ is the smallest space interval and c is the maximum wave velocity in the domain. For this study, the time step size is fixed by ∆t = 10 -9 s.

Results and discussions Fig. 3 (left) presents the V FAS obtained for all of 800 Monte-Carlo realizations. Two levels of dispersion δ = 1 and δ = 0.3 have been considered. It can be seen that when δ = 0.1, the values of V FAS oscillate around the one corresponding to the mean model which is shown as the red continuous line (V mean FAS = 3606 m.s -1 ). When the dispersion is higher (δ = 0.3), the measured values of V FAS are globally decrease. It may be explained by the fact that when the heterogeneity of the medium is greater, the scattering phenomenon becomes more significant and may modify the global wave velocity in the domain. In addition, we obtained a greater dispersion of V FAS is more important in the case δ = 0.3 than the one in the case δ = 0.1.

In Fig. 3 (right), the probability density functions of V FAS are computed for both cases. One may state that the probability density function of V FAS strongly depends on the dispersion δ. It means that in the practice, neglecting the uncertain heterogeneity in bone may lead to a poor prediction of the mean value of mechanical properties by using V FAS as a index.

Conclusion

The heterogeneity of bone's mechanical properties is significant and need to be considered when studying the sensitivity of ultrasound response in cortical long bone. As the statistical data on real bone material is hardly found, a parametric probabilistic method, which is based on the maximum entropy principle, has been used to generate an optimal probabilistic model for taking into account the uncertainties of bone elasticity. A explicit FDTD solver has been developed for simulating the wave propagation in a transversely isotropic heterogeneous medium in the time domain. It has been shown that the FAS velocity is very sensitive to dispersion of the bone's elasticity tensor in statistic sense. Detailed study need to be carried out and will be presented in a forthcoming paper.

Although the uncertain elastic model of the plate presented in this paper was developed to study behavior of cortical long bone, this procedure is of wider interest and may be applied to the characterization of other materials.

Fig. 1 .

 1 Fig. 1. Geometrical configuration of the trilayer model for ultrasound axial transmission test

Fig. 2 .

 2 Fig. 2. Spatial variation of the elasticity tensor components A ij

Fig. 3 .

 3 Fig. 3. (left) Values of V FAS evaluated at each realization; the solid line is obtained for the homogeneous material properties; (right) probability density functions of VFAS

  , e 2 , e 3 ) be a reference Cartesian frame where (e 1 , e 2 , e 3

	2 Modeling of transient wave propagation in a
	heterogeneous bone immersed in fluid
	2.1 Geometrical configuration
	Let R(O; e 1


Karine Beauchard
email: karine.beauchard@cmla.ens-cachan.fr

Paulo Sérgio Pereira Da Silva

Pierre Rouchon
email: pierre.rouchon@mines-paristech.fr

Stabilization of an arbitrary profile for an ensemble of half-spin systems ⋆

Keywords: nonlinear systems, Lyapunov stabilization, LaSalle invariance, quantum systems, Bloch equations, ensemble controllability, infinite dimensional system

We consider the feedback stabilization of a variable profile for an ensemble of non interacting half spins described by the Bloch equations. We propose an explicit feedback law that stabilizes asymptotically the system around a given arbitrary target profile. The convergence proof is done when the target profile is entirely in the south hemisphere or in the north hemisphere of the Bloch sphere. The convergence holds for initial conditions in a H 1 neighborhood of this target profile. This convergence is shown for the weak H 1 topology. The proof relies on an adaptation of the LaSalle invariance principle to infinite dimensional systems. Numerical simulations illustrate the efficiency of these feedback laws, even for initial conditions far from the target profile.

Introduction

Ensemble controllability as introduced in [START_REF] Li | Ensemble control of Bloch equations[END_REF] is an interesting control theoretic notion well adapted to nuclear magnetic resonance (NMR) systems (see, e.g., [START_REF] Li | Control of inhomogeneous quantum ensembles[END_REF] and the reference herein). In [START_REF] Beauchard | Controllability issues for continuous spectrum systems and ensemble controllability of Bloch equations[END_REF] some controllability issues of such NMR systems are investigated using openloop controls involving Dirac-combs. In [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF] such open-loop Dirac-combs are combined with Lyapunov stabilizing feedback to ensure closed-loop convergence towards a target profile that is one of the two steady-states, the south and north poles of the Bloch sphere. In this note, we extend this Lyapunov design to arbitrary target profiles and prove its local convergence for weak H 1 topology when the target pro-file lies entirely in the south hemisphere or in the north hemisphere.

We consider an ensemble of non interacting half-spins in a static field (0, 0, B 0 ) t in R 3 , subject to a transverse radio frequency field (ũ 1 (t), ũ2 (t), 0) t in R 3 (the control input). The ensemble of half-spins is described by the magnetization vector M ∈ R 3 depending on time t but also on the Larmor frequency ω = -γB 0 (γ is the gyromagnetic ratio). It obeys to the Bloch equation:

∂M ∂t (t, ω) = (ũ 1 (t)e 1 + ũ2 (t)e 2 + ωe 3 ) ∧ M (t, ω), (1)
where -∞ < ω * < ω * < +∞, ω ∈ (ω * , ω * ), (e 1 , e 2 , e 3 ) is the canonical basis of R 3 , ∧ denotes the wedge product on R 3 . The equation ( 1) is an infinite dimensional bilinear control system. The state is the ω-profile M , where, for every ω ∈ (ω * , ω * ), M (t, ω) ∈ S 2 (the unit sphere of R 3 ). The two control inputs ũ1 and ũ2 are real valued.

We propose here a first answer to the local stabilization of an arbitrary profile: given an arbitrary target profile M f : (ω * , ω * ) → S 2 , define an explicit control law (ũ 1 (t, M ), ũ2 (t, M )), a neighborhood U of M f (in some space of functions to be determined), a diverging sequence of times (t n ) n∈N , such that, for every initial condition M 0 ∈ U , the solution of the closed loop system is uniquely defined and satisfies

lim n→+∞ M (t n , .) -M f (.) L ∞ (ω * ,ω * ) = 0.
In this note, the Lyapunov feedback proposed in [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF] is adapted to provide a constructive answer to this question. Section 2 is devoted to control design and closed-loop simulations. In section 3 we state and prove the main convergence result, theorem 1.

2 Lyapunov H 1 approach

Some preliminaries

Let us recall the concept of a solution for (1) when the control input u contains Dirac distributions. When ũ1 , ũ2 ∈ L 1 loc (R), then, for every initial condition

M 0 ∈ L 2 ((ω * , ω * ), R 3 ), the equation (1) has a unique weak so- lution M ∈ C 0 ([0, +∞), L 2 ((ω * , ω * ), R 3 )). Denote by δ(t -a) the Dirac distribution located at t = a. When ũ1 = αδ(t -a) + u ♯
1 and ũ2 = u ♯ 2 where u ♯ j ∈ L 1 loc (R), α > 0 and a ∈ (0, +∞), then the solution is the classical solution on [0, a) and (a, +∞), it is discontinuous at the time t = a, with an explicit discontinuity given by an instantaneous rotation of angle α around the axis Re 1

M (a + , ω) =     1 0 0 0 cos(α) -sin(α) 0 sin(α) cos(α)     M (a -, ω).
The symbol . (resp. ., . ) denotes the Euclidian norm (resp. scalar product) on R 3 and the associated operator norm on M 3 (R).

Transformation into a driftless system

As in [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF] we consider a control with an "impulse-train" structure

ũ1 = u 1 + +∞ k=1 π δ(t -kT ), ũ2 = (-1) ǫ(t) u 2 (2)
where ǫ(t) := E(t/T ), for some period T > 0 and E(γ) denotes the integer part of the real number γ. The new controls u 1 , u 2 belong to L 1 loc (R). Considering the change of variable M 1 (t, ω) := P (t)M (t, ω) where

P (t) :=     1 0 0 0 ǫ(t) 0 0 0 ǫ(t)    
(3) one gets the following dynamics

∂M 1 ∂t (t, ω) = [u 1 (t)e 1 + u 2 (t)e 2 + ǫ(t)ω] ∧ M 1 (t, ω). (4)
The application of impulses at t = kT , by changing the sense of rotation of the null input solution, is expected to reduce the dispersion in the closed loop system. Since M (t, ω) = M 1 (t, ω) for every t ∈ [2kT, (2k + 1)T ], any convergence result on M 1 (t) when t → +∞ provides a convergence result on M .

The first step of the control design consists in putting the system (4) in driftless form. The new function

M 2 (t, ω) := exp[σ(t)ωS]M 1 (t, ω)
where

σ(t) := t 0 ǫ(s)ds, S :=     0 1 0 -1 0 0 0 0 0     , (5)
solves

∂M 2 ∂t (t, ω) = 2 i=1 u i (t) exp(σ(t)ωS)e i ∧ M 2 (t, ω). (6)
Since σ(2kT ) = 0, ∀k ∈ N, any convergence on M 2 (t) when t → +∞ provides a convergence on M 1 (2kT ) when k → +∞.

Transformation of the target profile

The second step of the control design consists in transforming a convergence to a variable profile M f into a convergence to the constant profile -e 3 , for which we developed tools in the previous work [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF]. It relies on the following proposition.

Proposition 1 There exists C > 0 such that, for all

M f ∈ H 1 ((ω * , ω * ), S 2 ), there exists R ∈ H 1 ((ω * , ω * ), SO 3 (R)) satisfying R(ω)M f (ω) = -e 3 , ∀ω ∈ [ω * , ω * ], (7)
R H 1 C M f H 1 . (8) Proof: Let M f ∈ H 1 ((ω * , ω * ), S 2 ) and set f (ω) := M ′ f (ω) ∧ M f (ω). Denote by A(ω) the skew-symmetric operator defined by R 3 ∋ M → f (ω) ∧ M ∈ R 3 . Con- sider the Cauchy problem d dω R = RA(ω) on [ω * , ω * ] with R(ω * ) = R * where R * is any rotation sending M f (ω * ) to -e 3 : R * M f (ω * ) = -e 3 . Since ω → A(ω) is L 2 the solution R is well defined, unique and belongs to H 1 ((ω * , ω * ), S 2 ). Direct computations show that d dω (RM f ) = 0. Thus R(ω)M f (ω) ≡ -e 3 . Moreover, R(ω) = 1 and R ′ (ω) = A(ω) = M ′ f (ω) for all ω ∈ [ω * , ω * ], which proves (8).
✷.

Let us consider a target profile

M f ∈ H 1 ((ω * , ω * ), S 2 ). Take R ∈ H 1 ((ω * , ω * ), SO 3 (R))
given by the above proposition. To any solution M 2 of (6), we associate the function

N (t, ω) := R(ω)M 2 (t, ω), ∀ω ∈ [-ω * , ω * ]. ( 9
)
This function solves the equation

∂N ∂t (t, ω) = 2 i=1 u i (t) F (t, ω)e i ∧ N (t, ω) (10)
where

F (t, ω) := R(ω) exp(σ(t)ωS). ( 11
)
The convergence of N (t, ω) to -e 3 as t → +∞ is equivalent to the convergence of M 2 (t, ω) to M f (ω) as t → +∞.

Lyapunov feedback

Let us consider the following Lyapunov-like functional

L(N ) := N + e 3 2 H 1 2 = ω * ω * 1 2 ∂N ∂ω 2 + 1 + N, e 3 dω. (12)
The function L is defined for any N ∈ H 1 ((ω * , ω * ), S 2 ) and takes its minimal value on this space at the point N = -e 3 with L(-e 3 ) = 0. For any solution of (10), some computations show that

dL dt [N (t)] = 2 i=1 u i (t)H i [t, N (t)]
where, for i = 1, 2 one has

H i [t, N ] := ω * ω * dN dω (ω), ∂F ∂ω (t, ω)e i ∧ N (ω) + e 3 , F (t, ω)e 1 ∧ N (ω) dω.
Hence, with the feedback laws

u i (t, N ) := -H i [t, N ], ∀i ∈ {1, 2}, (13)
it follows that

dL dt [N (t)] = -u 1 (t, N ) 2 -u 2 (t, N ) 2 0. ( 14
)
As in [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF], we have the following result.

Proposition 2 For every initial condition N 0 ∈ H 1 ((ω * , ω * ), S 2 ), the closed loop system (10), ( 13) has a unique solution

N ∈ C 1 [0, ∞), H 1 (ω * , ω * ), R 3 such that N (0) = N 0 .

Closed-loop simulations

We assume here ω * = 0, ω * = 1 and we solve numerically the T -periodic system (1) with the feedback law (ũ 1 , ũ2 ) given by ( 2), ( 13). The closed-loop simulation is per-

formed for t ∈ [0, T f ], T f = 20T and T = 2π/(ω * -ω * ). The ω-profile [ω * , ω * ] ∋ ω → (x(t, ω), y(t, ω), z(t, ω)) is discretized {1, . . . , N + 1} ∋ k → (x k (t), y k (t), z k (t))
with a regular mesh of step ǫ N = ω * -ω *

N

with N = 100. In other words, one has a set of discrete values {ω i , i = 1, . . . , N + 1}, where

ω i = ω * + (i -1)ǫ N .
We have checked that the closed-loop simulations are almost identical for N = 100 and N = 200. In the feedback law ( 16), the integral versus ω is computed assuming that (x, y, z) and (x ′ , y ′ , z ′ ) are constant over ](k

-1 2 )ǫ N , (k + 1 2 )ǫ N [, their values being (x k , y k , z k ) and x k+1 -x k-1 2ǫN , y k+1 -y k-1 2ǫN , z k+1 -z k-1 2ǫN
. The obtained differential system is of dimension 3(N + 1). It is integrated via an explicit Euler scheme with a step size h = T /1000. We have tested that h = T /2000 yields almost the same numerical solution at t = T f = 20T . After each time-step the new values of (x k , y k , z k ) are normalized to remain in S 2 . The initial ω-profile M 0 (ω) of (x, y, z) ∈ S 2 is given by x 0 = 0, y 0 = -1z 2 0 , where

z 0 = -cos( π 8 ) + 0.05 1 -cos( π 8 ) cos(ω π 2 ) . The desired final profile M f (ω) is given by x f = -1 -z 2 f , y f = 0, where z f = -cos( π 16 ) + 0.1 1 -cos( π 16 ) sin(ω π 4 ) .
The map R(ω) is constructed for the discrete set {ω i , i = 1, . . . , N + 1}, in the following way. For Figures 1 and 2 summarize the main convergence issues for these choices of initial profile M 0 and of the desired final profile M f . The convergence speed is rapid at the beginning and tends to decrease at the end. We start with L(0) ≈ 0.1929. We get L(20T ) ≈ 0.0032. This numerically observed convergence is confirmed by Theorem 1 here below. 3 Main Result

Local stabilization

Theorem 1 For every

M f ∈ H 1 ((ω * , ω * ), S 2 ) with M f (ω), e 3 = 0, ∀ω ∈ [ω * , ω * ], (15)
there exists δ 1 > 0 such that, for every N 0 ∈ H 1 ((ω * , ω * ), S 2 ) with N 0 + e 3 H 1 ≤ δ 1 , the solution of the closed loop system (10), ( 13) with initial condition N (0, ω) = N 0 (ω) satisfies N (t) ⇀ -e 3 weakly in H 1 (ω * , ω * ) when t → +∞.

The above theorem has the following corollary.

Corollary 1 For every M f ∈ H 1 ((ω * , ω * ), S 2 ) with ( 15), there exists δ 2 > 0 such that, for every M 0 ∈ 2. Initial (t = 0) and final (t = T f ) ω-profiles for x, y and z solutions of the closed-loop system 1 with the feedback (2), ( 13).

H 1 ((ω * , ω * ), S 2 ) with M 0 -M f H 1 < δ 2 , the solution of the system (1) with the initial condition M (0, ω) = M 0 (ω) and the feedback law given by ( 2), ( 13) satisfies M ((2kT ) + ) ⇀ M f weakly in H 1 when k → +∞. In particular,

lim k→+∞ M ((2kT ) + , .) -M f L ∞ (ω * ,ω * ) = 0.
The remaining part of this section is devoted to the proof of Theorem 1

LaSalle invariant set

The first step of our proof consists in checking that, locally, the invariant set is reduced to {-e 3 }.

Proposition 3 For every M f ∈ H 1 ((ω * , ω * ), S 2 ) with ( 15), there exists δ > 0 such that, for every

N 0 ∈ H 1 ((ω * , ω * ), S 2 ) with N 0 + e 3 < δ, the map t → L[N (t)] is constant on [0, +∞) if and only if N 0 = -e 3 .
Proof: Let us assume that t → L[N (t)] is constant. Then u 1 = u 2 = 0 and N (t, ω) ≡ N 0 (ω) (see ( 14) and ( 10)). Thus, for every j ∈ {1, 2} and t ∈ [0, +∞)

0 = ω * ω * N ′ 0 (ω), ∂F ∂ω (t, ω)e j ∧ N 0 (ω) + e 3 , F (t, ω)e j ∧ N 0 (ω) dω. ( 16
) For t ∈ [0, T ], σ(t) = t so F (t, ω) = ∞ k=0 t k ω k k! R(ω)S k and ∂F ∂ω (t, ω) = ∞ k=0 t k ω k k! R ′ (ω)S k + ∞ k=1 t k ω k-1 (k-1)! R(ω)S k .
Developing ( 16) in power series expansions of t and using (5), we obtain, for every j ∈ {1, 2} and k 1,

ω * ω * N ′ 0 (ω), ω k k! R ′ (ω) + ω k-1 (k-1)! R(ω) e j ∧ N 0 (ω) + e 3 , ω k k! R(ω)e j ∧ N 0 (ω) dω = 0.
By linearity, the following equality holds, for every

Q ∈ R[X] and j ∈ {1, 2} ω * ω * N ′ 0 (ω), (Q(ω)R ′ (ω) + Q ′ (ω)R(ω)) e j ∧ N 0 (ω) + e 3 , Q(ω)R(ω)e j ∧ N 0 (ω) dω = 0.
(17) Thanks to the density of polynomial functions in H 1 ((ω * , ω * ), C), the previous equality holds for ev-

ery Q ∈ H 1 ((ω * , ω * ), R). Let us recall the relations X, Y ∧ Z = Y, Z ∧ X and M X, Y = X, M ⊤ Y , ∀X, Y, Z ∈ R 3 and M ∈ M 3 (R)
, where M ⊤ denotes the transposed matrix of M . Then, the equality (17) may also we written

ω * ω * Q(ω)e j , R ′ (ω) ⊤ [N 0 (ω) ∧ N ′ 0 (ω)] + Q ′ (ω)e j , R(ω) ⊤ [N 0 (ω) ∧ N ′ 0 (ω)] + Q(ω)e j , R(ω) ⊤ [N 0 (ω) ∧ e 3 ] dω = 0. (18)
for every Q ∈ H 1 ((ω * , ω * ), R) and j ∈ {1, 2}. By linearity, we deduce that

ω * ω * Q(ω), R ′ (ω) ⊤ [N 0 (ω) ∧ N ′ 0 (ω)] + R(ω) ⊤ [N 0 (ω) ∧ e 3 ] + Q ′ (ω), R(ω) ⊤ [N 0 (ω) ∧ N ′ 0 (ω)] dω = 0 (19) for every Q ∈ H 1 ((ω * , ω * ), V)
where V := Span(e 1 , e 2 ). Let P : R 3 → V be the orthogonal projection on V. The previous equality is equivalent to

PR(ω) ⊤ [N 0 (ω) ∧ e 3 -N 0 (ω) ∧ N ′′ 0 (ω)] = 0 in H -1 R(ω) ⊤ [N 0 (ω) ∧ N ′ 0 (ω)] = 0 at ω = ω * and ω * .
(20) Here, H -1 denotes the dual space of H 1 0 (ω * , ω * ) for the L 2 -scalar product; the first equation has to be understood in the distribution sens. Thanks to N 0 (ω) ≡ 1, we have R(ω) ⊤ [N 0 (ω) ∧ N ′ 0 (ω)] ≡ N ′ 0 (ω) . Thus, the second line of ( 20) is equivalent to N ′ 0 = 0 at ω * and ω * . Notice that PR(ω) ⊤ | V is bijective on V for every ω ∈ [ω * , ω * ]. Indeed, thanks to (7) and (15), we have

Range[PR(ω) ⊤ | V ] ⊥ = Ker[PR(ω)| V ] = {v ∈ V; R(ω) ∈ Re 3 } = V ∩ RM f (ω) = {0}. Moreover, (PR(ω) ⊤ | V ) -1 ∈ H 1 , thus (20) gives -N ′′ 0 ∧ e 3 + N 0 ∧ e 3 = g in H -1 ((ω * , ω * ), V), N ′ 0 ∧ e 3 = 0 at ω * , ω * , (21)
where

g(ω) := -(PR(ω) ⊤ | V ) -1 PR(ω) ⊤ [N ′′ 0 (ω) ∧ (N 0 (ω) + e 3 )].
Therefore, there exists

C 1 = C 1 (ω * , ω * ) > 0 such that N 0 (ω) ∧ e 3 H 1 C 1 g H -1 . (22)
Thanks to (8), there exists

C 2 = C 2 (ω * , ω * , M f H 1 ) > 0 such that g H -1 C 2 N 0 + e 3 2 H 1 .
When N 0 is close enough to -e 3 in H 1 , then N 0 ∧e 3 H 1 and N 0 + e 3 H 1 are equivalent norms and then ( 22) gives

N 0 (ω) ∧ e 3 H 1 C 3 N 0 (ω) ∧ e 3 2 H 1 for some constant C 3 = C 3 (ω * , ω * , M f H 1 ) > 0. This implies N 0 ∧ e 3 = 0, i.e. N 0 = -e 3 . ✷
Remark 1 For M f ≡ e 1 , any constant function N 0 with values in Span(e 2 , e 3 ) belongs to the invariant set (see (20)). Thus, an assumption of the type ( 15) is required for our strategy to work.

Convergence proof

For the proof of Theorem 1, we need the following result.

Proposition 4 Take M f ∈ H 1 ((ω * , ω * ), S 2 ) and R ∈ H 1 ((ω * , ω * ), SO(3)) as in Proposition 1. Let (N 0 n ) n∈N a sequence of H 1 ((ω * , ω * ), S 2 ) and N 0 ∞ ∈ H 1 ((ω * , ω * ), S 2 ) such that N 0 n ⇀ N 0 ∞ weakly in H 1 and N 0 n → N 0 ∞ strongly in L 2 . Let α ∈ [0, 2T ] and (τ n ) n∈N be a sequence of [0, 2T ) such that τ n → α. Let N n (resp. N ∞ )
be the solutions of the closed loop system (10), ( 13) associated to the initial condition +∞) and n ∈ N so there exists M 0 > 0 such that N n (t) H 1 M 0 , for every t ∈ [τ n , +∞) and n ∈ N. The function t ∈ R → F (t, .) ∈ H 1 ((ω * , ω * ), M 3 (R)) defined by ( 11) is continuous and 2T -periodic, thus, there exists M 1 > 0 such that F (t, .) H 1 M 1 , for every t ∈ R. Thanks to (13), we have |u j [t, N n (t)]| 2M 1 M 0 , for every n ∈ N and t ∈ [τ n , +∞). We deduce from (10) that ∂Nn ∂t (t) L 2 4M 1 M 0 for every t ∈ [τ n , +∞) and n ∈ N. The end of the proof is as in [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF]. ✷

N n (τ n ) = N 0 n (resp. N ∞ (α) = N 0 ∞ ). Then, we have N n (t) ⇀ N ∞ (t) weakly in H 1 , ∀t > α, and u j [t, N n (t)] → u j [t, N ∞ (t)], ∀t > α, ∀j ∈ {1, 2}. Proof: The sequence (N 0 n ) n∈N is bounded in H 1 and L[N n (t)] L[N 0 n ], for every t ∈ [τ n ,
Proof of Theorem 1: The proof is as in [START_REF] Beauchard | Stabilization for an ensemble of half-spin systems[END_REF]. One may replace Barbalat's lemma by the Lebesgue reciprocal theorem, in the following way. Thanks to (14), t → u i [t, N (t)] belongs to L 2 (0, +∞), thus, for any diverging sequence of times k n , the sequence (t ∈ (0, +∞) → u i [2k n T + t, N (2k n T + t)]) n∈N converges to zero in L 2 (0, +∞). Therefore, there exists a subset N ⊂ (0, +∞) with zero Lebesgue measure such that u j [t, N (2k n T + t) → 0 for every t ∈ (0, +∞) -N and j ∈ {1, 2}. ✷

Concluding remark

Open-loop "impulse-train" control are combined with Lyapunov feedback to steer an initial profile [ω * , ω * ] ∋ ω → M (0, ω) of the Bloch-sphere system (1) towards an arbitrary target profile [ω * , ω * ] ∋ ω → M f (ω). Convergence is proved to be local for any target profile belonging either to the south or to the north hemisphere. We guess that our convergence proof could be extended to the case where M f intersects transversely the equator and thus where M f is not confined in only one hemisphere.

Fig. 1 .

 1 Fig. 1. Lyapunov function L(t) defined by (12) and the closed-loop control (u1, u2) defined by(13)

Acknowledgements

KB and PR were partially supported by the "Agence Nationale de la Recherche" (ANR), Projet Blanc C-QUID number BLAN-3-139579 and Project Blanc EMAQS number ANR-2011-BS01-017-01. PSPS is Partially supported by CNPq -Conselho Nacional de Desenvolvimento Cientifico e Tecnologico, and FAPESP-Fundação de Amparo a Pesquisa do Estado de São Paulo.

i = 1, one takes r 3 (ω 1 ) = M 0 (ω 1 ). Now choose a vector θ among the vectors of the canonical basis in a way that θ, M f (ω 1 ) is the minimum value. Construct r 2 (ω 1 ) = 1 θ∧r3(ω1) (θ ∧ r 3 (ω 1 )). Then one may take r 1 (ω 2 ) = r 2 (ω 2 ) ∧ r 3 (ω 2 ). Now, for i = 2, 3, . . . , N + 1 one chooses r 3 (ω

, and so on. The orthogonal matrix R(ω) formed by the column vectors r 1 , r 2 , r 3 is then transposed to obtain R(ω).


Bo Hou

Philippe Grosso
email: philippe.grosso@telecom-bretagne.eu

B Hou1

P Grosso1

J L De Bougrenet De La Tocnaye1

Marc Le Menn3

Principle and implementations of a refracto-nephelo-turbidimeter for seawater measurements

come

Introduction

The study of physical properties of seawater requires to measure temperature, salinity and pressure with high resolutions and accuracies. In order to describe the medium, chemical and optical properties need also to be measured. Turbidity is one of them. Measurement sensors can be mounted on drifting floats as ARGO ones [START_REF] Stewart | Introduction to Physical Oceanography[END_REF]. In this case, the development of in situ multifunctional sensors is valuable. In general, a set of these properties at different geographical positions is most valuable for the study of oceanography. This makes the research of in-situ multi-functional sensors valuable. Among these properties, attention has been focused on salinity and turbidity. In fact, temperature, salinity, and pressure are used to calculate density which variations create ocean currents. Concerning optical properties, turbidity is one of the most frequently measured quantities as it allows to quantify living and sedimentary particles concentrations and the penetration depth of light in the medium. Oceanographic studies are generally made in coastal zones or in open ocean areas. The resolution requirements for measuring salinity and turbidity are quite different in these zones. To keep a good sensitivity, transmissiometers must have a long path-length and scatterometers a high resolution. In the same way, in order to detect small density variations, salinity measurements need a high resolution. In coastal areas, suspended particles concentrations are higher and transmission path length must be smaller to avoid multiple reactions and to keep valuable the conditions of the Beer-Lambert relation. Due to the fact that coastal areas are generally characterized by high temperature and salinity gradients, resolution specifications on salinity measurements are lower. According to these requirements and application areas, it is meaningful to study the design of a high resolution in-situ salino-turbidi-meter for different scenarios. Previous works have demonstrated an optical refractometer to measure the salinity of seawater based on laser beam deviation measurement with a 1-D Position Sensitive Device (PSD) [START_REF] Malarde | High-resolution and compact refractormeter for salinity measurements[END_REF]. This refractometer has been containerized and tested in real condition (2000 meters in open ocean) [START_REF] Le Menn | Advances in measuring ocean salinity with an optical sensor[END_REF]. The resolution for measuring salinity reaches 0.002 g/kg with a measurement range from 0 g/kg to 40 g/kg. Further research to substitute Charge-Coupled Device (CCD) to PSD improved the resolution of nearly 1.5 times compared to PSD-based refractometer [START_REF] Hou | Chargecoupled devices combined with centroid algorithm for laser beam deviation measurements compared to a position-sensitive device[END_REF]. For turbidity measurement, we have analyzed the possibility of measuring salinity and turbidity simultaneously based on our refractometer from the same water sample and studied the interaction between turbidity and salinity measurements [START_REF] Hou | A turbidi-meter based on a refractometer using a charge-coupled device[END_REF]. A part of this correlation is caused by the combination of transmitted and scattered light in turbid medium, which can't be separated by PSD. That is one of the reasons why PSD is inappropriate to measure the light beam deviation. Another artifact is related to the gaussian beam projection angle. When the gaussian beam projected on a plane at a non-perpendicular angle, it forms an asymmetric gaussian beam seen by the sensor as a center shift. According to these two interactions, it is mandatory to measure the turbidity for salinity measurement in a turbid medium. In this paper, we propose two new high resolution in-situ salino-turbidi-meters derived from these works. A brief summary of the salinity and turbidity measurements as well as their correlations are discussed in section 2. According to the requirements and challenges pointed out in section 2, two new salino-turbidi-meters are presented in section 3. A further discussion on configurations and applications of these two salino-turbidimeter is carried out in section 4.

Requirements and Technical Challenges

Principle of refractometer

UNESCO/IOC SCOR/IAPSO working group 127 (WG127) [START_REF] Ioc | The international thermodynamic equation of seawater 2010: Calculation and use of thermodynamic properties[END_REF] has defined the measurement of seawater refractive index as one of the methods to measure the absolute salinity of seawater. Figure 1 shows the principle of a basic refractometer. A light source generates a beam of light, which propagates in the prism with a known refractive index n1. At the interface between the prism and the seawater sample, the light beam is redirected into the seawater and deviated to a light position sensor at position B, which measures the light spot position. With this light spot position, according to the Snell-Descartes law, the refractive index n2 of the seawater sample can be calculated. Our previous work included the design and implementations of a refractometers with a PSD and a V-shape refraction surface [START_REF] Malarde | High-resolution and compact refractormeter for salinity measurements[END_REF]. The resolution of this refractometer reached 10 -7 . It should be mentioned that the refractive index measurement depends on the seawater temperature varying in the range -2°C to 4 0°C. The variation of temperature impacts the measurement: Sellmeier coefficient: the variation of temperature changes the wavelength of the light source, which further impacts the refractive index of the prism. Thermo-optical coefficient: the variation of temperature changes the refractive index of the prism. Thermal expansion coefficient: the variation of temperature changes the dimension of the prism.

These facts require a refractometer design with a temperature compensation feature to minimize the refractive index change caused by temperature variation. Previously in our refractometer, we used two prisms, one had positive thermo-optical coefficients while another one had the negative thermo-optical coefficient, to compensate the refractive index change caused by the temperature variation. This compensation works in quasi-static conditions as in dynamic conditions because the two prisms are always at the same temperature as they are submitted to the same temperature variations and as their thermal mass is the same. The composition of glasses doesn't change significantly their density and then their weight. A miniature design significantly reduces both temperature expansion and thermal inertia.

Principle of turbidity measurement

Turbidity is the cloudiness or haziness of a fluid caused by individual particle. Turbiditymeter measures transmitted or scattered light at different angles, as shown in figure 1. At point C, a light sensor is placed at an angle α to the light propagation direction. When light propagates in the turbid water sample, the light is scattered in all the directions. The portion of light detected at an angle α could be used to describe the turbidity of the water sample. Among these measurement angles, two of them are the most used: 0° and 90°. At 0°, transmission of light is measured a s shown in point B, while scattering is measured at 90° are called nephelometer. The most advanced turbidi-meters combine the nephelometer with the measurement at several different angles and combine them in a ratio of nephelometric values and a weighted sum of values from all the detection angles to assess turbidity [START_REF] Ioc | The international thermodynamic equation of seawater 2010: Calculation and use of thermodynamic properties[END_REF]. We have verified the feasibility of measuring turbidity based on a refractometer from the transmission direction (0°). The turbidity T can be calculated from the following equation:

        = 0 ln 1 I I l T tr (1)
where Itr is the transmitted light intensity measured at the transmission direction, I0 is the original light source intensity and l is the light path in water. From this equation, it is clear that the reference light source intensity needs to be measured and controlled with an incertainty better than 1% to calculate turbidity. According to this equation, the turbidity derivative T' can be expressed as:

tr tr I l dI T × = ' (2)
From ( 2), to improve the resolution of this method, the light path length l in water can be increased, using a higher resolution light intensity sensor, or increasing the incident light intensity. We have verified this point by several simulations and experiments based on our refractometers [START_REF] Hou | A turbidi-meter based on a refractometer using a charge-coupled device[END_REF]. A challenge is that the light intensity measured by the detector at 0° is a combination of the transmitted light and the scattered light. Furthermore, the use of laser module as the light source causes interference by the particles in the water, which further disturbs the measurement of transmitted light intensity. The only available information from a PSD is the gravity center and total power of the light spot so that light intensity distribution would need to be measured in order to separate the transmitted light from the mixed light. This can be made with a CCD sensor. Therefore, a CCD with a low pass band filter has been preferred to eliminate the disturbance caused by interference and separate the transmitted light as shown in the figure 2.

Correlations between refractive index and turbidity measurement

Another issue related to the combination of transmitted and scattered light is that the PSD calculates the gravity center of the incident light, which includes both the gravity center of the transmitted spot and the scattered light. This results in an error when using PSD to measure the light beam deviation in a turbid medium. Due to this reason, PSD is not appropriate for a refractometer in a turbid medium. Besides the combination of transmitted light and scattered light, other factors impact the salinity measurement in a turbid medium as well. According to the principle of refractometer illustrated in section 2.1, the refractive index is measured by the light beam deviation. However, we have found that when a Gaussian beam projected to a plane with an angle of α as shown in figure 1, the projection spot is still a Gaussian spot with the same spot size but with a center shift. The shift distance is proportional to the turbidity of the medium, the squared spot size and inversely proportional to the tangent of the angle α. We established in [START_REF] Hou | A turbidi-meter based on a refractometer using a charge-coupled device[END_REF] the following relation:

(

In which

2 2 2 ) ( 0 σ p x x e I - -
is a standard gaussian beam, T is the attenuation coefficient, l is the light path of the center, u is the coefficient that describes the different attenuation inside the beam, p x is the original peak position and σ is the laser beam size. From (3), it is easy to conclude that the new laser spot is still a gaussian spot but with a peak shift of :

u T 2 σ - (4)
Because of this new gaussian center variation law, both refractive index and turbidity contribute to the light beam deviation, which results in an error to measure the refractive index from the light beam deviation. Based on this correlation, it is necessary to combine the turbidity and salinity measurements with the same seawater sample and the turbidity need to be calculated first in order to accurately indicate the salinity. Furthermore, when the gaussian beam is not collimated, the projection spot with a non-perpendicular angle becomes asymmetrical, which impacts the gravity-center-based spot location method, such as PSD and centroid algorithm.

A new method, which tracks the location that divides the mass of the spot into two equal parts, introduced in our previous study [START_REF] Hou | A turbidi-meter based on a refractometer using a charge-coupled device[END_REF] has been proved to be a good solution to this problem. Another correlation between the salinity and turbidity measurements is that the suspended particles have different refractive index from the seawater. A part of light refracted by the particles is captured by the position sensor as well, which impacts the light beam deviation. However, it should be mentioned that the light beam deviation caused by the particle-refracted light is much smaller than the required beam deviation resolution so that it is not a key issue for the refractive index measurement. Besides the impact to salinity measurement caused by turbidity, the configuration of refractometer affects the turbidity measurement as 1), the attenuation of the transmitted light is related to the geometrical light path length l in water. However, the light path length changes according to different refractive indices of seawater samples. For two different seawater samples with same turbidity, but different refractive index, the measured light intensities at 0° are different due to different light path lengths. That is to say, the design of the refracto-turbidi-meter needs to minimize the light path length variation caused by the refractive index changes.

Refracto-turbidi-meter requirements

According to the principle of salinity and turbidity measurement methods as well as the correlation between them, it is necessary to combine these two measurements with the same seawater samples. This results in some requirements to guide the design of refracto-turbidi-meter. (i) Reduce refractive index variation caused by the temperature change (ii) Minimize the sensor size to reduce the thermal expansion (iii) Maximize the light path length in water to improve the turbidity resolution (iv) Measure the original light intensity as the reference to calculate turbidity (v) Select different light source types according to the particle size in order to avoid interference (vi) Use CCD to eliminate the disturbance caused by the combination of interference and scattered light (vii) Use ratio nephelometer to improve the turbidity resolution (viii) Keep the beam projected on the position sensor at an angle of 90° in order to avoid the beam deviation caused by turbidity (ix) Minimize the light path length variation caused by refractive index change.

Therefore, any design of refracto-turbidi-meter should take the considerations of these requirements.

Design and Implementation

According to the requirements and challenges discussed in last section, a new miniature high resolution refracto-turbidi-meter is designed. Figure 3 depicts the new design of a refracto-turbidi-meter. This new model only has one piece of prism. The light beam first enters the prism from the left corner and reaches the first refraction surface AO. It is then refracted into the medium and refracted again at the refraction surface BO. This refraction directs the beam to a mirror CD, which redirects the upwards beam to the right side of the prism, where it meets the second mirror EF and redirects to a 20 mm long groove at the right side of the prism. Finally, the beam propagates through the groove and reaches the right bottom of the prism, where a CCD is used to measure the beam deviation caused by the water refractive index variation. At the inner surface of the groove, a surface treatment is used to absorb the reflected light to avoid the influence of the scattering inside the groove.

Refractive index measurement

As discussed in section 2.1, temperature has a strong impact on the in-situ refractive index measurement. In our design, we use a single prism to compensate the refractive index change caused by the temperature variation. The refractive index is not only associated to temperature but also to wavelength, illustrated by Sellmeier equation [START_REF] Ghosh | Handbook of thermo-optic coefficients of optical materials with applications[END_REF]:

3 2 2 3 2 2 2 2 1 2 2 1 2 ) ( ) ( ) ( ) ( ) ( ) ( 1 )) ( ( C T T B C T T B C T T B T n - + - + - + = λ λ λ λ λ λ λ (5)
When the temperature increases, the wavelength of the laser as well, which further results in the decrease in glass refractive index, according to [START_REF] Hou | A turbidi-meter based on a refractometer using a charge-coupled device[END_REF]. The variation caused by the wavelength change ∆nλ < 0 (when ∆T > 0) can be used to compensate the refractive index variation induced by temperature change ∆nT (when ∆T < 0). To choose a glass with this feature, we plotted two functions, dn/dT = f1(T) and dn/dλ = f2(T), for several Schott glasses at a wavelength of 635 nm in 20 °C, and found that the glass N-ZK7 best fits our requirements. Figure 4 plots the two functions and the refractive index variation from 0 _C and 40 _C for N-ZK7. As shown in figure 4, the refractive index variation of N-ZK7 from 0 °C and 40 °C is only 7.5x 10 -7 .

Figure 4. The temperature coefficient and Sellmeier coefficient of glass N-ZK7 when temperature varies from 0 °C to 40 °C Furthermore, this new design has a very small size prism (40 mmx30 mmx20 mm), resulting in a very small thermal expansion. The refractive index measurement resolution of this design can be evaluated by using the geometric method and differential method to build the relationship between the spot position and the refractive index. According to the calculation results, the relationship between the spot position P and seawater refractive index ne can be approximately obtained by a curve polynomial fitting. The equation of the fitting curve of the second order is shown in equation 4, with a standard deviation of 4.86 x10 -7 . ne = 0.000032 P 2 + 0.006 P + 1.34 [START_REF] Ioc | The international thermodynamic equation of seawater 2010: Calculation and use of thermodynamic properties[END_REF] When the refractive index changes from 1.330 to 1.350, the spot moves a distance of about 2.2 mm from -1.1 mm to +1.1 mm. According to the spot location resolution (better than 0.1 µm) checked in our previous work by using a CCD [START_REF] Hou | A turbidi-meter based on a refractometer using a charge-coupled device[END_REF], the resolution of measuring seawater refractive index reaches 0.9 x10 -6 (≈1 ppm) which is equivalent to the absolute salinity of 2mg.kg -1 . Another issue related to the refractive index measurement is the Gaussian spot center shift according to the projection angle (inversely proportional to the tangent of projection angle). To minimize this center shift, the light beam enters the groove at an angle of 90°± 0.865° with the range of refractive index from 1.33 to 1.35.

Turbidity measurement

The turbidity measurement uses the measurement in the transmission direction. A CCD is used to capture the transmitted light intensity. As discussed in section 2.2, the transmitted light is mixed with the scattered light and interference, a CCD with a low pass filter can eliminate the scattered light and interferences. The turbidity resolution based on the transmitted light measurement is inversely proportional to the light path length in water. In our previous work, we have obtained a resolution of 1:15 NTU (Nephelometric Turbidity Units) in the range from 0 to 20 NTU. The corresponding light path length in water is about 6.35 mm. The turbidity measurement resolution of this new design can be calculated as: Rt = 1.15x6.35/ (20 +MN), in which MN is the length of light path between the two refraction surfaces. This distance is directly affected by the incident light position d0 as shown in Fig. 3. The light path distance between two refraction surfaces decreases from 11 mm to 1.8 mm, when the incident beam position moves from 3 mm to 8 mm. If we use 5 mm as the incident beam position d0, the expected turbidity measurement resolution would be Rt = 1.15x6.35/(20 +MN) = 0.27 NTU in the range from 0 NTU to 20 NTU. This design also reduces the impacts caused by the light path difference with different seawater refractive indices. The variation of the light path length between two refraction surfaces is about 0.054 mm, which is only 0.23% of the entire light path, so that it can be ignored when calculating turbidity. Another improvement is the set-up of light source. In our previous refractometers, light source was directly in contact with to the prism, making it difficult to measure the reference light intensity I0. In the new model, light source is installed with the mechanics, a light splitter is used to separate the beam into two parts: one enters into the prism for the measurement, another is measured by a photo-diode as the reference light intensity I0. This design makes the use of different light sources other than laser possible. The laser module shown in Fig. 3 can be replaced by a LED at wavelength 860 nm, which is not absorbed by the natural particles and reduces interference. Another option is 530 nm, which is used for the calibration of submarine visibility model and has less absorption by the seawater. Besides these advantages brought by the light source set-up, the separation of the light source and the prism is a mandatory configuration for the temperature compensation. As mentioned in section 3.1, the refractive index variation caused by the temperature change is compensated by the positive temperature coefficient variation and negative Sellmeier coefficient variation with the same temperature change, which requires that the light source and measured sample keep the same temperature. This separation enables a good thermal isolation to avoid the different temperatures, caused by the light source heating, between the light source and the prism. The attenuation based turbidity measurement is not sensitive in ultra low and ultra high turbid medium. To improve the performance of the refracto-turbidi-meter, a new photo-diode is introduced to measure the scattered light at 90° as shown in figure 5. This new refracto-nepheloturbidi-meter equips with a mirror at 45° to the light propagation direction at the right side of the groove, which reflects the diffused light to the photo-diode B placed inside the mechanics. With this configuration, the turbidity Tu is calculated by the photodiode B (I90) and the light intensity obtained by the CCD (Iccd) as:

CCD I I I Tu 2 90 1 90 α α + = (7)
In which α1 and α2 are two constant coefficients. It should be noticed that the reference light intensity Iref given by the photo-diode A is not used in the calculation. The reason is that equation ( 7) is a ratio between the 90° scattered light and the weighted sum of transmitted light and scattered light, which is independent of the incident light intensity. However, both the nephelometer standard ISO 7027 [START_REF]Water quality -determination of turbidity[END_REF] and USEPA Method 180.1 [START_REF] Apha | Standard methods for the examination of water and wastewater[END_REF] only have a measurement range from 0 to 40 NTU. The photo-diode B introduced in this design is used to improve the performance in ultra low turbid medium, for the moderate turbidity measurement, the attenuation based measurement can be used. Based on this reason, we keep the photo-diode A to measure the reference intensity for the attenuation based turbidity measurement. This turbidi-meter is expected to have the same resolution as the one proposed in last section (1% of the measurement range), but owns better performance in ultra low turbidity medium.

Device Customization

As we have discussed in the introduction section, the measurement in different zones of the ocean has different requirements. Accordingly, the salino-turbidimeterintroduced in this paper should have different configurations to fulfill the different application scenarios. In open ocean, the seawater is clear enough so that it is not necessary to measure turbidity, and there is no impact to the salinity measurement caused by the turbidity, such as Gaussian center shift. Thus, a configuration without the groove can be used to measure salinity in open ocean. A PSD or a CCD, can be used as the light position detection device, because there is no influence caused by turbidity. In offshore area, the turbidity stays in a low level (usually smaller than 20 NTU), the refracto-turbidi-meter with the 45° mirror design can be used. The ratio calculated by ( 5) can help the turbidity measurement obtain a resolution of 1% of the measurement range. Boss and W. Scott Pegau have pointed out that 117° is a better turbidity detection angle for the ocean turbidity measurement [START_REF] Boss | Relationship of light scattering at an angle in the backward direction to the backscattering coefficient[END_REF], it is possible to change the angle of the 45° mirror to 31.5° in order to reflect the scattered light at 117°. As the turbidity impacts t he refractive index measurement, using a CCD as the light position detection device is mandatory for the measurement of light beam deviation. Different from the offshore and open ocean cases, the coastal areas have a high turbidity, which might be out of the nephelometer measurement range (usually 40 NTU). The CCD based transmissometer can be used to measure the turbidity in the high turbidity case. Furthermore, the use of CCD combined with the advanced post-processing algorithms, the influence to salinity measurement caused by turbidity can be reduced.

Conclusion

The study of physical oceanography requires the in-situ measurement of salinity and turbidity with the same seawater sample. The refractometer based salinity measurement requires that the sensor has smaller thermal expansion and less temperature influence. To obtain a higher resolution turbidity measurement, the original light intensity must be measured as the reference and the sensor should contain a long light path length in water. The turbidity impacts the combination of transmitted light and scattered light at the light propagation direction. This added to the light interference caused by the suspended particles, when using a laser, make PSD inappropriate in a turbidity medium to measure the salinity. Furthermore, the Gaussian beam will have a center shift in a turbid medium when the beam projects to a plane at a non-perpendicular angle. This proves that the turbidity disturbs the salinity measurements. According to these requirements and challenges, we have designed a new refracto-turbidi-meter equipped with a CCD to reduce the turbiditycaused impacts to the salinity measurements. To improve the turbidity measurement performance in low turbid medium a refracto-nephelo-turbidi-meter is designed as well. According to the measurement principle and previous work we designed, the new refracto-turbidi-meter with different configurations can work in different application scenarios, including open ocean, offshore, and coastal areas. The expected salinity measurement resolution reaches 2mg.kg -1 and the expected turbidity measurement resolution is 1% of the measurement range in both the low turbid medium and the high turbid one.

Figure 1 .

 1 Figure 1. Principle of refractometer and the turbidity measurement

Figure 2 .

 2 Figure 2. The captured scattered light and interference by a camera, and the corresponding images after low pass filter

  shown in equation (

Figure 3 .

 3 Figure 3. Prototype of miniature refracto-turbidi-meter

Figure 5 .

 5 Figure 5. The prototype of miniature refracto-nephelo-turbidi-meter


G N Nguyen
email: giang.nguyen@telecom-bretagne.eu

K Heggarty

P Gerard

P Meyrueis

Iterative scalar algorithm for the rapid design of wide-angle diffraction Fourier elements

We demonstrate that non-paraxial Fourier elements can be designed by an iterative Fourier transform algorithm with the help of a simple projection step. The element designed by this model shows better reconstruction in both spots position and power distribution than the one of the paraxial design and very close to the desired pattern.

Introduction

Recent development in fabrication technology has enabled manufacturing of nonparaxial Diffractive Optical Elements (DOEs), which could have higher efficiency, better uniformity and an ever broadening field of applications. However, the scalar paraxial theory fails to model these DOEs correctly, leading to the need for wideangle diffraction modeling and design. More rigorous design algorithms have been developed [START_REF] Kress | Applied Digital Optics[END_REF], but their applications are often limited due to the extensive calculation required. It has been reported that the diffraction pattern on a hemisphere in the farfield can be calculated accurately using a simple Fourier transform [START_REF] Harvey | [END_REF].

We aim to reduce computational complexity of the design model by developing a rapid algorithm using a projection of a desired pattern from the output plane to the hemisphere and iterative Fourier transform between the hemisphere and the DOE plane. The design can be used for non-paraxial Fourier elements which have many applications in practice, e.g. beam splitters, beam shapers, where the propagation distance is often much larger than the DOEs size [START_REF] Kress | Applied Digital Optics[END_REF]. The reconstruction pattern of the DOE designed by our model shows better distribution in both position and intensity than the one designed by a standard multi-stage iterative Fourier transform algorithm (IFTA) [3] and very close to the desired pattern.

Discussion

For non-paraxial scalar diffraction, a Fourier transform can be used to propagate the diffracted field to a hemisphere in the far-field [START_REF] Harvey | [END_REF], but not directly to the observation plane as in the scalar paraxial theory. The diffraction pattern on the output plane can be obtained by using a simple spherical wave projection from the hemisphere, resulting in:

ܷሺߙ, ߚ; ݖሻ = ߛ ଶ expሺ݅݇ߛ/ݖሻ ݅λݖ ݔܶ{ܷሺܨ ଵ , ݕ ଵ ; 0ሻ}
where ߙ, ߚ, ߛ are the direction cosines of the angles between the diffracted wave and the spatial coordinates. A re-sampling step is then required to interpolate the field in the output plane from angular ܷሺߙ, ߚ; ݖሻ to spatial coordinates ܷሺݔ ଶ , ݕ ଶ ; ݖሻ on a uniform sampling grid. For paraxial diffraction, the hemisphere and the output plane overlap approximately and this projection step can be neglected. However, the projection is more important for wide-angle diffraction, as shown in Fig. 1(a) where the reconstruction pattern is distorted compared to the desired pattern, which is a square grid with the diffraction angle of about 30 o at the corners. The element is designed using a paraxial IFTA with a smallest pixel size of 1μm and operating at the wavelength of 1.55μm, while the reconstruction pattern is obtained using the rigorous Rayleigh-Sommerfeld (RS) diffraction formula [START_REF] Goodman | Introduction to Fourier Optics[END_REF]. The other distortion is in intensity distribution, where the higher the diffraction angle, the less the spots power. This is due to the power estimation in the paraxial design (ܲ ∝ |ܷ| ଶ ) without taking into account the non-paraxial angle between the direction of diffracted beam with the normal of the observation plane (ܲ ∝ |ܷ| ଶ ߛ). To solve these problems for the design of non-paraxial Fourier elements, the desired pattern is back projected from the observation plane to the hemisphere before using iterative Fourier transform between the DOE plane and the hemisphere. Fig. 1(b) illustrates the projected pattern in angular coordinates on the hemisphere, where the spot intensity has been corrected according to the diffraction angle. A nearest neighbour interpolation is used for mapping the field from spatial to angular coordinates. The element is then designed by using the standard multi-stage IFTA having this image as the output pattern. Fig. 1(c) demonstrates the superposition of the desired pattern and the RS reconstruction on the observation plane of our designed DOE, which shows almost the same spatial and intensity distribution.

Conclusion

We present an efficient scalar model using the traditional iterative Fourier transform with the help of a spherical projection to design non-paraxial Fourier elements, while the calculation complexity is much simpler than the available rigorous models. The simulated results show that the DOE designed by this algorithm performs better reconstruction and uniformity than the one of the paraxial IFTA design and very close to the desired pattern. The DOEs are being fabricated for experimental verification.

Fig. 1 .

 1 (a) Result of the traditional IFTA. The red pattern is the desired output, whereas the green pattern is the reconstruction of the designed DOE. The yellow regions are where they overlap. (b) Projection of the desired pattern from the observation plane to the hemisphere. (c) Superposition of the desired pattern and the reconstruction of the DOE designed by our proposed algorithm