M Blanc

M Touratier

Modelling Elastic and Thermoelastic Thick Multilayered Composites by a New Constrained Discrete Layer Approach

A new layer by layer/discrete layer model is presented for the linearly elastic and thermoelastic responses of multilayered composite plates. The transverse variation of the displacement field is defined in terms of a one-dimensional finite element representation, using quadratic Lagrangian interpolations. The laminate thickness is subdivided into a series of one-dimensional finite elements (i.e., thickness subdivisions) whose nodes correspond to planes of constant transverse normal coordinates in the undeformed configuration of the laminate. Each of the displacements and loads is expanded in a double Fourier series in the cartesian surface coordinates. Interface and boundary conditions are exactly verified and exploited to reduce the number of independent generalized displacements, before solving the boundary value problem from the standard variational principle for displacements. This leads to a new model, whose size 3N in 3D, is less than the classic layer by layer approach requiring 3(2N + 1), also in 3D.

In addition, this model allows satisfying exactly all the interface and boundary conditions. Numerical evaluations and comparisons with exact 3D elasticity solutions confirm the efficacy of the proposed approach. Finally, the proposed modelling has been extended to the thermoelastic response of laminates and the results are compared with the exact 3D solutions. Temperature distributions required for the thermoelastic analysis are determined from the constrained discrete layer approach previously applied to heat conduction.

INTRODUCTION

Composite laminates are being increasingly used in structural applications and they require efficient analysis tools for their optimal design and reliability study. This field has been very extensively studied for several decades and is still a subject of research interest, as it is extremely difficult to find efficient modelling methods for multilayered structures with many layers and thickness zones due to the coupling between mechanics, heat transfer and moisture, and containing localized regions of complex loading or geometric or material discontinuities. Thus, Address correspondence to M. Touratier, LMSP UMR 8106 CNRS-ENSAM-Université d'Orléans, 151 bd de l'Hôpital, 75 013 Paris, France. E-mail: touratier@paris.ensam.fr several approaches have been developed, based mainly on either equivalent single layer or layer by layer descriptions, depending on the accuracy required for numerical simulations [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells: Theory and Analysis[END_REF]. Unlike the equivalent single layer theories, the layer by layer theories assume separate displacement field expansions within each material layer, thus providing a kinematically correct representation of the strains in discrete layer laminates, and allowing more accurate ply-level stresses to be determined. Therefore, only layer by layer models can lead to consistent descriptions of localized effects. Extensive recent reviews on modelling multilayered plates and shells may be found in Reddy [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells: Theory and Analysis[END_REF], Reddy and Arciniega [START_REF] Reddy | Shear deformation plate and shell theories: Stavsky to present[END_REF], and Carrera [START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF][START_REF] Carrera | Temperature profile influence on layered plates response considering classical and advanced theories[END_REF][START_REF] Carrera | Closed-form solutions to assess multilayered-plate theories for various thermal stress problems[END_REF].

In this paper, our interest concerns layer by layer descriptions along the thickness of the laminate based on quadratic Lagrangian interpolation functions. The transverse variation of the three components of the displacement field is thus defined in terms of one-dimensional finite element representations. The laminate thickness dimension is subdivided into a series of M one-dimensional quadratic finite elements (i.e., thickness subdivisions) at least equal to the number of layers. The proposed modelling incorporates the transverse normal strain, to be able to accurately determine interlaminar stresses by satisfying exactly all interface plus top and bottom boundary conditions from standard variational formulations, not using 3D equilibrium equations for transverse stress distribution determinations, as, for example, in Robbins and Reddy [START_REF] Robbins | Modelling of thick composites using a layerwise laminate theory[END_REF]. Pioneer work in layer by layer approximations such as one-dimensional finite element representation along the thickness of the laminate was published by Reddy [START_REF] Reddy | A generalization of two-dimensional theories of laminated composite plates[END_REF], and it is termed the "layerwise" theory. Applications of the layerwise theory to piezoelectric laminates can be found in Saravanos et al. [START_REF] Saravanos | Layerwise mechanics and finite element for dynamic analysis of piezoelectric composite plates[END_REF], Garcia Lage et al. [START_REF] Garcia Lage | Layerwise partial mixed finite element analysis of magneto-electro-elastic plates[END_REF], and Semedo Garção et al. [START_REF] Semedo Garção | Analysis of laminated adaptive plate structures using layerwise finite element models[END_REF], and for heat conduction see Blanc and Touratier [START_REF] Blanc | A constrained discrete layer model for heat conduction in laminated composites[END_REF]. Several different layer by layer approaches, not using one-dimensional finite element representation through thickness of laminates, are of course available in the literature (see Carrera [START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF]). A fully coupled thermoelastopiezoelectric three-dimensional solution has been presented by Xu et al. [START_REF] Xu | Three-dimensional solutions for coupled thermoelectroelastic response of multilayered plates[END_REF], using a mixed variational theorem for temperature, normal heat flux, displacements, normal transverse and transverse shear stresses, electric potential and normal electric displacements. Finally, exact three-dimensional solution for thermoelastic analysis of laminates is given in Tungikar and Rao [START_REF] Tungikar | Three-dimensional exact solution of thermal stresses in rectangular composite laminate[END_REF], while Savoia and Reddy [START_REF] Savoia | Three-dimensional thermal analysis of laminated composite plates[END_REF] carried out a three-dimensional thermomechanical analysis for laminated composite plates. In Matsunaga [START_REF] Matsunaga | A comparison between 2-D single-layer and 3-D layerwise theories for computing interlaminar stresses of laminated composite and sandwich plates subjected to thermal loadings[END_REF], comparisons between 2D single layer and 3D layerwise theories are presented for computing interlaminar stresses in laminates subjected to thermal loadings.

The objective of the present paper is to explain a new layer by layer (or discrete layer) description in which interface continuity and top and bottom boundary conditions are exploited to satisfy exactly all these conditions, and to drastically reduce model size, compared to any existing layer by layer approaches based on classic or mixed variational theorems. Thus, using quadratic Lagrangian approximation throughout thickness of laminate including N layers, the new approach of so called "Constrained Discrete Layer" modelling, requires, for 3D problems, choosing for example only one subdivision per layer, 3N independent generalized displacements, while standard layer by layer descriptions require 3(2N + 1) unknowns, and 6(2N + 1) in a mixed modified approach if normal and shear transverse stresses are not eliminated. Note that in this last case, these stresses are quadratic too, which, of course, is not the case when using a classic variational formulation and quadratic approximation throughout the layer thickness for displacement. Specifically, the objective of this paper is to present analytic three-dimensional, doubly periodic solutions for the static response of laminated composites, consisting of a number of perfectly bonded linearly elastic/thermoelastic and fiber-reinforced layers. The solutions are periodic in the surface coordinates x 1 and x 2 , and, as explained above, of quadratic variation along the thickness coordinate x 3 = z of the laminate. Periodicity for displacement and load along in-plane coordinates, allows satisfying exactly all the interface and boundary conditions within the classic displacement variational theorem. See for example Touratier [START_REF] Touratier | A refined theory for thick composite plates[END_REF], concerning difficulties solving the boundary value problem associated with thick lam-FIG. 1. The composite laminate configuration.

inates, when using a given approximate displacement field to find closed solutions. Any desired degree of displacement variation through the thickness is easily obtained by adding more one-dimensional subdivisions through the thickness. The Constrained Discrete Layer modelling is then compared with exact elasticity calculations from Pagano [START_REF] Pagano | Exact solutions for composite laminates in cylindrical bending[END_REF][START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF]. It is confirmed that this new approach gives good results for both displacement and stress distributions through the thickness of the laminate, and allows significantly reducing the size of the problem to be solved.

Responses in thermoelasticity for a rectangular laminate have been analyzed and compared with exact three-dimensional calculations from Tungikar and Rao [START_REF] Tungikar | Three-dimensional exact solution of thermal stresses in rectangular composite laminate[END_REF]. Once again, the present constrained discrete layer modelling gives accurate distributions for temperatures, displacements and thermal stresses through the thickness of a simply supported square sandwich plate, submitted to both thermal and mechanical loading.

FORMULATION OF THE BOUNDARY VALUE PROBLEM IN LINEAR ELASTICITY

Expansions for Displacement Field and Traction

Components Along Mid-Plane Coordinates A laminated plate occupying the parallelepipedic domain V of lateral edge e , constituted by N homogeneous anisotropic layers V (i) is considered, see Figure 1. The plate is referred to an orthogonal coordinate system = {O/x α , x 3 = z}, where the x α (α = 1, 2) axes lie in reference plane A of the plate and x 3 = z is in the transverse direction. The top and bottom faces of the plate and the layer interfaces, located at x 3 = x 3t , x 3b and x 3i (i = 1, N-1), are denoted by t , b and i , respectively. The plate has combinations of homogeneous cross-ply orthotropic layers of constant thickness, L 1 and L 2 are the side lengths in the x 1 and x 2 directions and h is the total thickness of the plate.

The traction components F i (i = 1 to 3) at the top and bottom surfaces of the plate, are functions of the x 1 and x 2 , i.e.:

F t i = F t i (x 1 , x 2 ) F b i = F b i (x 1 , x 2 ) (1)
where the superscripts t and b refer to the top and bottom surfaces, respectively. For laminates consisting of cross-ply layers, each of the displacement components u i is expressed as the sum of products of trigonometric functions in the Cartesian surface coordinates. The trigonometric functions are chosen so that the response quantities are periodic in x 1 and x 2 with periods 2L 1 and 2L 2 , as:

u = ∞ r=0 ∞ s=0 T (rs) (x 1 , x 2 )U (rs) (z) with U (rs) T (z) = U 1 (rs) (z), U 2 (rs) (z), U 3 (rs) (z) T (rs) =   cos a r x 1 sin b s x 2 0 0 0 sin a r x 1 cos b s x 2 0 0 0 sin a r x 1 sin b s x 2   (2)
and where a r = rπ L 1 , b s = sπ L 2 ; r and s being integers. Therefore, the following displacement conditions are satisfied by Eq. (2) as:

x 1 = 0, L 1 ; u 2 = u 3 = 0 x 2 = 0, L 2 ; u 1 = u 3 = 0 (3) FIG.
2. Thickness subdivisions of the laminate and generalized displacements.

The traction components at the top and bottom surfaces in Eq. ( 1) are also expanded in double Fourier series as follows

       F b = r s T (rs) F b rs F t = r s T (rs) F t rs F qT (rs) = F q 1(rs) F q 2(rs) F q 3(rs) ; q = b, t (4)

Discrete Layer Approach Along the Thickness of the Laminate

The laminate being subdivided into N discrete layers, see Figure 2, from Eq. (2) a continuous displacement field per layer is then written in the thickness direction by

U (rs) = N j=1 χ ]z j-1 ,z j [ U (j) (rs) χ ]z j-1 ,z j [ = 1 if z ∈]z j-1 , z j [, 0 elsewhere (5)
Introducing the non-dimensionalized thickness coordinate ξ i for each layer V (j) as in Figure 2, the generalized displacements U (j) i(rs) per layer are approximated along the thickness of the layer by quadratic Lagrange polynomia l b: bottom; m: medium; t: top

j k as      U (j) (rs) = k∈{b,m,t} l j k U (j) k (rs) = l j b U (j) b (rs) + l j m U (j) m (rs) + l j t U (j) t (rs) j = 1, . . . , N l j k = l j k (ξ j ) (6) where U
It is directly deduced from the above equations that the continuity of the displacement field at the layer interfaces is then automatically satisfied, and we obtain

U (j) (rs) = l j b U (j) b (rs) + l j m U (j) m (rs) + l j t U (j+1) b (rs) (8)
Thus, the displacement field at any point of the laminate is finally written from Eqs. (2), ( 5) and ( 8)

as u = N j=1 ∞ r=0 ∞ s=0 χ ]z j-1 ,z j [ T (rs) l j b U (j) b (rs) + l j m U (j) m (rs) + l j t U (j+1) b (rs) (9)

The New So-Called "Constrained Discrete Layer" Approach: Interface and Top and Bottom Conditions for Stresses; Constraints and Final Displacement Field

The novelty of the present approach called "constrained discrete layer" is in taking into exactly account both interface continuity, and traction on the top and bottom conditions of a laminate. Hereafter, these conditions are exactly verified to reduce the number of unknowns (number of independent generalized displacements U (j) k i(rs) ). Denoting σ (l) ij the stresses per layer, the above conditions on stresses, respectively at interfaces and at bottom and top surfaces, yield

σ (j) 3i (x 1 , x 2 , z j ) = σ (j+1) 3i (x 1 , x 2 , z j ) j = 1, . . . , N -1; i ∈ {1, 2, 3} (10)
σ (1 or N) 3i (x 1 , x 2 , 0 or h) = δ i3 F b or t 3 δ i3 Kronecker's symbol (11)
Introducing the orthotropic elastic constitutive law per jth layer, one obtains the stresses per layer as following

            σ (j) 11 σ (j) 22 σ (j) 33 σ (j) 23 σ (j) 31 σ (j) 12             =             C (j) 1111 C (j) 1122 C (j) 1133 0 0 0 C (j) 1122 C (j) 2222 C (j) 2233 0 0 0 C (j) 1133 C (j) 2233 C (j) 3333 0 0 0 0 0 0 C (j) 2323 0 0 0 0 0 0 C (j) 3131 0 0 0 0 0 0 C (j) 1212                         u (j) 1,1 u (j) 2,2 u (j) 3,3 u (j) 2,3 +u (j) 3,2 u (j) 1,3 + u (j) 3,1 u (j) 1,2 + u (j) 2,1             (12)
where C

(j) lnpq are the elastic constants of the material of the jth-layer; u (j) i are components per layer of the displacement field which, from Eq. ( 9), reduces to

u (j) (rs) = T (rs) l j b U (j) b (rs) + l j m U (j) m (rs) + l j t U (j+1) b (rs) (13)
Using expressions for stresses in jth-layer deduced from Eqs. ( 12) and ( 13), into Eq. ( 10), gives interface conditions, and bottom and top boundary conditions

-at interfaces Ĉ(j)T (rs)i3 U (j) b (rs) + Ĉ (j+1)T (rs)i3 U (j+1) b (rs) - Ĉ (j+2)T (rs)i3 U (j+2) b (rs) = - B (j)T (rs)i3 U (j) m (rs) + B (j+1)T (rs)i3 U (j+1) m (rs) i = 1, 2, 3; j = 1, 2, . . . , N -1 (14)
-at bottom and top surfaces

Ĉ (b 1 or t N )T (rs)i3 U (1 or N) b (rs) + Ĉ (b 2 or t N+1 )T (rs)i3 U (2 or N+1) b (rs) = δ i3 F b or t 3 - B(b 1 or t N )T (rs)i3 U (1 or N) m (rs) i = 1, 2, 3 (15)
Coefficients of submatrices Ĉ(l) (rs)i3 and B(n) (rs)i3 are given in Appendix. These Eqs. ( 14) and ( 15) allow us to eliminate all the generalized displacements U (rs) . Therefore, all the conditions expressed by Eqs. ( 14) and ( 15) can be written under the following algebraic system in U b (rs) :

Ĉ (rs) U b (rs) = B (rs) U m (rs) + F 3(rs) (16)
In this last algebraic system, coefficients of matrices Ĉ (rs) and B (rs) are extracted from Eqs. ( 14) and ( 15). In addition, we have established:

U bT (rs) = U (1) b T (rs) U (2) b T (rs) • • • U (N+1) b T (rs) U mT (rs) = U (1) m T (rs) U (2) m T (rs) • • • U (N) m T (rs) (17)
F T 3(rs) = 0 0F b 3(rs) 0 • • • 00 0F t 3(rs)
Then, the system given by Eq. ( 16) is solved in respect to U b (rs) . This leads to a final displacement field of the following form:

u = N j=1 ∞ r=0 ∞ s=0 χ ]z j-1 ,z j [ T (rs) l j b Ĉ -1 (rs) B (rs) U m (rs) + F 3(rs) j th row + l j m U (j) m (rs) + l j t Ĉ -1 (rs) B (rs) U m (rs) + F 3(rs) j+1 th row ( 18
)
Remark. At this stage, and adopting quadratic polynomial approximations for displacements in each layer, it is important to observe that this constrained discrete layer model based upon Eq. ( 18), leads only to 3N generalized displacements to be determined for three-dimensional laminates having N layers. In comparison, the standard discrete layer theory (based on displacement variational theorem) requires 3(2N + 1) unknowns, while the discrete layer modelling based on modified Reissner's mixed variational theorem involves 6(2N + 1) unknowns, but transverse (both normal and shear) stresses are quadratic [START_REF] Carrera | Single-layer vs multi-layers plate modelings on the basis of Reissner's mixed theorem[END_REF]. It is of course possible to improve accuracy either by discretizing layers in sub-layers, keeping the same quadratic approximation for each, or employing higher order polynomia along the thickness.

Governing Equations

The equilibrium equations and natural boundary conditions on the edge of the laminate are derived in classic manner via the standard displacement variational theorem as:

0 = V C ijkl ∂ 2 u l ∂x j ∂x k + (f i = 0) δu i + V F e i -σ ij n j δu i ;
∀δu i in V, and ∀δu i on e [START_REF] Carrera | Single-layer vs multi-layers plate modelings on the basis of Reissner's mixed theorem[END_REF] where u l are components of u in Eq.( 18), F e i are components of prescribed traction on the lateral edge e of V, σ ij are stresses computed from Eq. ( 12) in each layer, n j are components of the unitary normal outward to the edge e , δu i is the variation operator (differential operator) of u i . Developing Eq. ( 19) for each pair (r, s) of harmonics in Eq. ( 18), and performing integration throughout the thickness of the laminate, involves solving a linear algebraic system having 3N unknowns U m i(rs) . Corresponding symbolic calculations are carried out using the Mathematica system [START_REF] Wolfram | Mathematica a System for Doing Mathematics by Computer[END_REF].

EXTENSION TO THE THERMOELASTIC RESPONSE OF LAMINATES

The above formulation of the boundary value problem is now extended to thermoelastic analyses from the constrained discrete layer approach. The described method for linearly elastic composite laminates is modified as follows to reach corresponding linearly thermoelastic responses:

Step-i: Start by expanding temperature T(x 1 , x 2 , z) as the component u 3 in Eq. (2), i.e.: T = r≥0 s≥0 sin a r x 1 sinb s x 2 T (rs) (z). Then, pursue this by exchanging u 3 with T until satisfying Eq. ( 9). Replacing Eq. ( 12) by Fourier's law available for heat conduction, and inserting it in interface and top and bottom conditions for normal heat flux similarly to Eqs. [START_REF] Semedo Garção | Analysis of laminated adaptive plate structures using layerwise finite element models[END_REF] and [START_REF] Blanc | A constrained discrete layer model for heat conduction in laminated composites[END_REF], we deduce the final temperature field in the following form:

T(x 1 , x 2 , z) = N j=1 ∞ r=0 ∞ s=0 χ ]zj-1,zj[ sin a r x 1 sin b s x 2 × N i=1 L (i) m (ξ j , κ (1) , . . . κ (N) )T (i)m (rs) + L (1) b (ξ j , κ (1) , . . . κ (N) )T (1)b (rs) + L (N+1) b (ξ j , κ (1) , . . . κ (N) )T (N+1)b (rs) (20)
where coefficients L (n) k depend on Lagrange polynomial, and also, both on their first derivative and on κ

(j) = k (j) 3 /(z j+1 -z j ), k (j) 3
being thermal conductivities in the z-direction for an orthotropic material. In fact, coefficients L (n) k are nothing else than quadratic B-splines.

Step-ii: Solve the heat conduction boundary value problem by using the constrained discrete layer modelling as described above in elasticity (Section 2), or from [START_REF] Blanc | A constrained discrete layer model for heat conduction in laminated composites[END_REF]. In any case, the appropriated classic variational form to be used is expressed as:

∀ δT such that δT = 0 on T , ∀ τ: τ 0 V ρC V T ,t δTdVdt = τ 0 V q i δT ,i dVdt+ τ 0 V ρwδTdVdt + τ 0 Q (Q + q.n)δT d dt ( 21
)
where Q is the heat quantity prescribed (natural boundary conditions) on the part Q of the boundary surface = t ∪ b ∪ e of the volume V of the laminate, q is the heat flux, n the unit normal outward from the boundary surface , C v is the specific heat at constant volume, ρ is the mass density, ρw is an internal source of heat, and δT is the variation of the temperature T, at any specific time t ∈ [0, τ]. Finally, T is the complementary part of Q , where temperatures are prescribed (essential boundary conditions). The solution of the heat problem is therefore deduced from Eq. ( 19), by incorporating Eq. ( 20).

Step-iii: Modify the constitutive law Eq. ( 12) by superimposing the standard thermal contribution in stresses for an orthotropic behaviour. This requires adding the following vector to the second member of Eq. ( 12):

           -α (j) 1 C (j) 1111 + α (j) 2 C (j) 1122 + α (j) 3 C (j) 1133 T (j) -T (j) 0 -α (j) 1 C (j) 1122 + α (j) 2 C (j) 2222 + α (j) 3 C (j) 2233 T (j) -T (j) 0 -α (j) 1 C (j) 1133 + α (j) 2 C (j) 2233 + α (j) 3 C (j) 3333 T (j) -T (j) 0 0 0 0            (22)
where T (j) (x 1 , x 2 , z) is the temperature in jth layer computed at the above iith-step, T (j) 0 is the reference temperature for the jth layer, while α (j) i are coefficients of thermal expansion in the jth layer;

Step-iv: Interface, and bottom and top surfaces conditions, respectively given by Eqs. ( 14) and [START_REF] Matsunaga | A comparison between 2-D single-layer and 3-D layerwise theories for computing interlaminar stresses of laminated composite and sandwich plates subjected to thermal loadings[END_REF], must incorporate the thermal contribution into stresses by adding the following terms to the second member of these equations: * to be added to Eq. ( 14):

δ i3 Â(j) (rs)i3 - Â(j+1) (rs)i3
, at z = zj (23) * to be added to Eq. ( 15):

δ i3 Â(b1 or tN) (rs)i3
, at z = z 1 or z N (24

)
where coefficients Â(k) (rs)i3 are expressed as:

Â(k) (rs)i3 = 3 p=1 α (k) p C (k) pp33 T (k) (rs) -T (k) 0(rs) ; k ∈ {j, j+1, b 1 or t N } (25)
while temperatures per layer T (k) (rs) being computed from Eqs. ( 20)- [START_REF] Ray | Optimal control of thin circular cylindrical laminated composite shells using active constrained layer damping treatment[END_REF]; and δ i3 is Kronecker's symbol;

Step-v: This results in a new vector F 3(rs) in Eq. ( 16);

Step-vi: The final displacement field for linearly thermoelastic responses is then deduced from Eq. ( 18), using the new expression of F 3(rs) obtained at the preceding vth-step;

Step-vii: Finally, corresponding equilibrium equations and natural boundary conditions are still computed from Eq. ( 19).

NUMERICAL STUDIES

To assess the accuracy of the present model, the following problems for which an exact three-dimensional solution exists (which has been implemented here) from Pagano [START_REF] Pagano | Exact solutions for composite laminates in cylindrical bending[END_REF][START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF], have been examined: first, the cylindrical bending of a simply supported plate and, second, the bending of a simply supported rectangular plate, being submitted respectively to single and double sinusoïdal loading at the top of the plate. Hence, in Eqs. ( 2) and ( 4), integers r and s are kept as r = s = 1; and then, simply supported conditions verify Eq. (3); while in Eq. ( 4) we retain F b i(11) = 0, i = 1, 2, 3, and F t 1(11) = F t 2(11) = 0, F t 3(11) = (p sin(πx 1 /L 1 ) sin(πx 2 /L 2 ) for rectangular and square plates; and p sin(πx 1 /L 1 ) in cylindrical bending. We also present results for a square sandwich plate (Problem 3), and for the cylindrical bending of a bi-layer laminate (Problem 4). Finally, a thermoelastic analysis is produced for a bi-layer submitted to both thermal and mechanical loadings (cylindrical bending for Problem 5 and three-dimensional plate in Problem 6). A comparison with exact calculations is also given.

Otherwise, linearly elastic material properties are those given in Pagano [START_REF] Pagano | Exact solutions for composite laminates in cylindrical bending[END_REF], and are:

E L /E T = 25;
G LT /E T = 0.5; G TT /E T = 0.2; ν LT = ν TT = 0.25; E T = 6.909 GPa where L and T are the parallel and normal directions to the fibres, E is Young's modulus, G is the shear modulus, and ν is Poisson's ratio.

All the symbolic and numerical computations for both our new model and the exact solutions were performed by using the Mathematica system [START_REF] Wolfram | Mathematica a System for Doing Mathematics by Computer[END_REF].

Problem 1. Cylindrical Bending of a Simply

Supported Three-Layer 90 • /0 • /90 • Symmetric Laminate under a Sinusoidal Load Results solving Eq. ( 19) are reported in Tables 1 and2 for two thickness ratios (4 and 10), under the following dimensionalized form [START_REF] Pagano | Exact solutions for composite laminates in cylindrical bending[END_REF]:

         S = L 1 h z = z h ū1 = E T u 1 (0,z) ph u 3 = 100E T u 3 L 1 2 ,z phS 4 ( σ11 , σ33 , σ13 ) = 1 p σ 11 L 1 2 , z , σ 33 L 1 2 , z , σ 13 (0, z) (26)
The layers have equal thickness, the fibers are oriented along the x 2 -direction for the outer layers (90 • ), and along the x 1 -direction for the inner layer (0 • ). In Tables 1 and2, CDLn means with "Constrained Discrete Layer" using n sub-layers for each plate's layer (see Ray and Reddy [START_REF] Ray | Optimal control of thin circular cylindrical laminated composite shells using active constrained layer damping treatment[END_REF][START_REF] Ray | Effect of delamination on the active constrained layer damping of laminated composite beams[END_REF] for constrained layer damping). We show in Figures 3 and4, corresponding thickness distributions for displacements and stresses, using n = 3 subdivisions per layer. These are compared with the exact solution from Pagano, and a good agreement can be observed on these thickness distributions, satisfying exactly all the interface and boundary conditions. Tables 1 and2 show that convergence with respect to the number n of sub-layers for CDLn modelling is very fast between CDL1 (only one subdivision per layer) and CDL2 (two subdivisions per layer), towards exact calculations. It appears that 3 subdivisions in each layer may suffice to reach exact values at studied points in these tables. [START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF], while geometric data are such that L 2 = 3 L 1 . Two thickness ratios are again considered for thick laminates: 4 and 10. The layers have equal thickness, with fibers in the outer layers oriented in the x 1 -direction (0 • ) and those in the inner layer oriented in the x 2 -direction (90 • ). The plate is subjected to the doubly sinusoïdal distribution transverse loading expressed above. Tables 3 and4 present non-dimensionalized displacements and stresses such as:

             S = L1 h z = z h ū1 = ETu1 0, L 2 2 ,z phS 3 ū2 = ETu2 L 1 2 ,0,z phS 3 ū3 = 100ETu3 L 1 2 , L 2 2 ,z phS 4 ( σ11 , σ22 , σ12 ) = 1 pS 2 σ 11 L1 2 , L2 2 , z , σ 22 L1 2 , L2 2 , z , σ 12 (0, 0, z) ( σ13 , σ23 ) = 1 pS σ 13 0, L2 2 , z , σ 23 L1 2 , 0, z σ33 = 1 p σ 33 L1 2 , L2 2 , z (27)
Results are also presented in Figures 5 and6 for two thickness ratios (4 and 10). These three-dimensional computations show ū1 (0) 0,0139 0,0152 0,0141 0,0140 ū1 (1/3) -0,0031 -0,0048 -0,0032 -0,0032 ū1 (2/3) 0,0039 0,0056 0,0040 0,0040 ū1 (1) -0,0142 -0,0156 -0,0145 -0,0143 ū2 (0) 0,0132 0,0126 0,0132 0,0132 ū2 (1/3) 0,0049 0,0047 0,0049 0,0049 ū2 (2/3) -0,0037 -0,0036 -0,0037 -0,0038 ū2 (1) -0,0124 -0,0117 -0,0124 -0,0124 -0,0248 -0,0252 -0,0250 -0,0249 σ22 (1/3) -0,1193 -0,1147 -0,1191 -0,1193 σ22 (2/3) 0,1088 0,1041 0,1086 0,1088 σ22 (1) 0,0400 0,0404 0,0402 0,0401 σ23 (1/2) 0,0334 0,0257 0,0358 0,0330 σ13 (1/2) 0,3511 0,5444 0,3477 0,3503 σ12 (0) 0,0281 0,0277 0,0282 0,0281 σ12 (1/3) 0,0060 0,0049 0,0060 0,0060 σ12 (2/3) -0,0038 -0,0027 -0,0038 -0,0038 σ12 (1) -0,0269 -0,0437 -0,0270 -0,0269 ū1 (0) 0,00920 0,00931 0,00921 0,00920 ū1 (1/3) 0,00135 0,00120 0,00134 0,00135 ū1 (2/3) -0,00128 -0,00113 -0,00128 -0,00128 ū1 (1) -0,00917 -0,00928 -0,00918 -0,00917 ū2 (0) 0,00475 0,00473 0,00475 0,00475 ū2 (1/3) 0,00165 0,00165 0,00165 0,00165 ū2 (2/3) -0,00147 -0,00146 -0,00147 -0,00147 ū2 (1) -0,00457 -0,00455 -0,00457 -0,00457 ū3 (1/2) 0,919 0,915 0,919 0,919 σ11 (0) -0,725 -0,734 -0,726 -0,726 σ11 (1/3) -0,106 -0,094 -0,105 -0,106 σ11 (2/3) 0,104 0,092 0,103 0,103 σ11 (1) 0,726 0,735 0,727 0,726 σ22 (0) -0,0122 -0,0123 -0,0122 -0,0122 σ22 (1/3) -0,0435 -0,0435 -0,0435 -0,0435 σ22 (2/3) 0,0418 0,0418 0,0418 0,0418 σ22 (1) 0,0145 0,0146 0,0146 0,0145 σ23 (1/2) 0,0152 0,0125 0,0162 0,0151 σ13 (1/2) 0,420 0,435 0,419 0,420 σ12 (0) 0,0123 0,0123 0,0123 0,0123 σ12 (1/3) 0,0033 0,0032 0,0033 0,0033 σ12 (2/3) -0,0030 -0,0029 -0,0030 -0,0030 σ12 (1) -0,0120 -0,0120 -0,0120 -0,0120 that the trend is identical to that in the two-dimensional case: very fast convergence properties for CDLn models between n = 1 and 2, optimal subdivisions for accurate calculations along the thickness being reached as soon as n = 3, two subdivisions in each layer may suffice among points reported in Tables 3 and4.

ū3 (1/2) 2,
As for Problem 1, Problem 2 shows the efficiency of the proposed "Constrained Discrete Layer (CDL)" model to accurately compute displacements and stresses for thick laminates, keeping at most 27 unknowns in Eq. ( 19) for the above (Problem 2) 3D problem with n = 3 discretizations in each layer. Note that CDL1 and CDL2 require, in Eq. ( 19), for the above 3D problem, 9 (n = 1) and 18 (n = 2) unknowns, respectively. Using classic discrete layer modelling (i.e., not satisfying Eqs. ( 14) and ( 15)) based on quadratic Lagrangian approximation throughout thickness of the laminate as in this paper, we need 21; 39 or 57 unknowns, respectively for the above 3D problem, in adopting 1, 2 or 3 subdivisions per layer, as in Problem 2 for CDL modelling. The reader may refer to Reddy [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells: Theory and Analysis[END_REF][START_REF] Reddy | A generalization of two-dimensional theories of laminated composite plates[END_REF] about classic layerwise modelling based on Lagrangian discretization along the thickness of laminates. Symmetric Square Plate under Sinusoidal Load Material data are same as those given for Problem 2 above [START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF], while geometric data are such that L 2 = L 1 . Two thickness ratios are again considered for thick laminates: 4 and 10. The layers have equal thickness, with fibers in the outer layers oriented in the x 1 -direction (0 • ) and those in the inner layer oriented in the x 2 -direction (90 • ). The plate is subjected to the doubly sinusoïdal distribution transverse loading expressed above. Tables 5 and6 present non-dimensionalized displacements and stresses as in above Eq. ( 27).

Results are also presented in Figures 7 and8 for two thickness ratios (S = 4 and 10). These three-dimensional computations for square sandwich laminate show that the trend is identical to that in the two-dimensional case: very fast convergence properties for CDLn models between n = 1 and 2, optimal subdivisions for accurate calculations along the thickness being reached as soon as n = 3, two subdivisions in each layer may suffice at several points reported in Tables 5 and6.

As for above Problems 1 and 2, the efficiency of the proposed "Constrained Discrete Layer (CDL)" model is estab- ū1 (0) 0,00936 0,01028 0,00950 0,00941 ū1 (1/3) -0,00242 -0,00356 -0,00245 -0,00244 ū1 (2/3) 0,00317 0,00438 0,00321 0,00319 ū1 (1) -0,00969 -0,01073 -0,00985 -0,00975 ū2 (0) 0,0228 0,0210 0,0228 0,0228 ū2 (1/3) 0,0072 0,0077 0,0072 0,0072 ū2 (2/3) -0,0066 -0,0071 -0,0067 -0,0067 ū2 ( 1 [START_REF] Pagano | Exact solutions for rectangular bidirectional composites and sandwich plates[END_REF]), a bi-directional two layer laminate is considered with the T = 90 • and L = 0 • directions aligned parallel to the x 1 -axis in the top and bottom layers, respectively, the layers being of equal thickness. Material properties and loading conditions are the same as those in Problem 1. Only the thickness ratio L 1 /h = 4 is tested in this Problem. Results are presented in Figure 9. Comparisons are shown between exact calculations and our CDL solution, by using either one, three, or four subdivisions in each layer. It is established for this difficult problem as well, that the proposed constrained discrete layer model converges towards the exact calculations, for thick laminates. But this problem as is known appears very difficult, because the convergence of stresses is slower than in the above problems.

Problem 5. Thermoelastic Response of a

Bi-Directional Two-Layer 90 • /0 • Laminate Subjected to Both Sinusoidal Thermal and Mechanical Loadings To evaluate our approach in thermoelasticity, we again consider the bi-layer analyzed in elasticity in the above Problem 4, but boundary conditions are now expressed as follows:

-mechanical boundary conditions:

σ 13 (x 1 , 0) = 0 u 3 (x 1 , 0) = 0 σ 13 (x 1 , h) = 0 σ 33 (x 1 , h) = p Sin πx 1 L 1 (28)
-thermal boundary conditions :

T (x 1 , 0) = -T 0 Sin πx 1 L 1 T (x 1 , h) = T 0 Sin πx 1 L 1 (29)
where T(L 1 /2, h)/T 0 = 1.

Mechanical properties have been given at the top of Section 4. From [START_REF] Tungikar | Three-dimensional exact solution of thermal stresses in rectangular composite laminate[END_REF], thermal properties are: k

L = 36.42 Wm -1• C -1 ; k T = 0.96 Wm -1• C -1 ; α L = 0.57 10 -6• C -1 ; α T = 35.6 10 -6• C -1 .
For this kind of supported antisymmetric bi-layer (90 • /0 • ) laminated composite subjected to the above thermal and mechanical loadings, we determine an exact thermoelastic solution from [START_REF] Tungikar | Three-dimensional exact solution of thermal stresses in rectangular composite laminate[END_REF] by using the Mathematica package, and a constrained discrete layer solution from the calculations summarized above in Eqs. ( 20) and ( 21) (also see [START_REF] Blanc | A constrained discrete layer model for heat conduction in laminated composites[END_REF]). Corresponding temperature distributions computed from Eqs. ( 20) and ( 21), are then shown in Figure 10. Results from the present constrained discrete layer modelling for linearly thermoelastic analysis are now compared with exact calculations for displacement and thermal stress distributions through the thickness of the bi-layer laminate as depicted in Figure 11. It is concluded again on good convergence properties for the constrained discrete layer approach, as soon as the number of subdivisions in each layer increases up to 6, having only one independent generalized temperature per FIG. 10. Temperature distribution along the thickness of the bi-directional bi-layer laminate in Problem 5, at its mid-length. subdivision, and two independent generalized displacements per subdivision.

Problem 6. Thermoelastic Response of a

Three-dimensional Two-Layer 90 • /0 • Laminate Subjected to Both Doubly Sinusoïdal Thermal and Mechanical Loads In this section, we consider the above laminate as a threedimensional medium. Boundary conditions are as follows.

-mechanical boundary conditions:

       σ 13 (x 1 , x 2 , 0) = σ 13 (x 1 , x 2 , h) = σ 23 (x 1 , x 2 , 0) = σ 23 (x 1 , x 2 , h) = 0 u 1 (x 1 , 0, 0) = u 1 (x 1 , L 2 , 0) = u 2 (0, x 2 , 0) = u 2 (L 1 , x 2 , 0) = 0 u 3 (x 1 , x 2 , 0) = 0 and σ 33 (x 1 , x 2 , h) = pSin πx1 L1 Sin πx2 L2 (30)
thermal boundary conditions:

     T(0, x 2 , x 3 ) = T(L 1 , x 2 , x 3 ) = T(x 1 , 0, x 3 ) = T(x 1 , L 2 , x 3 ) = 0 T(x 1 , x 2 , 0) = -T 0 Sin πx 1 L 1 Sin πx 2 L 2 T(x 1 , x 2 , h) = T 0 Sin πx 1 L 1 Sin πx 2 L 2 (31)
The temperature distribution for this problem has been obtained from Eqs. ( 20) and ( 21), and is shown in Figure 12. The thermoelastic analysis corresponding to this present problem is shown in Figure 13. Results from our CDL modelling are compared both with exact calculations and the standard layerwise model, called here the DL model (see [START_REF] Reddy | Mechanics of Laminated Composite Plates and Shells: Theory and Analysis[END_REF][START_REF] Robbins | Modelling of thick composites using a layerwise laminate theory[END_REF]). It is concluded that the CDL model is better than the DL on stress distributions, for a reduced computational model size. Finally, the proposed constrained discrete layer modelling allows reducing the size problem for laminates without any loss of accuracy, and gives very good thickness distributions for both displacements, stresses, temperatures, as both interface and top and bottom boundary conditions are exactly satisfied. The capability of this modelling in capturing particular transverse shear stress distributions for very thick laminates must be remarked (here, thickness ratio 4), see Figures 3 to 13.

CONCLUSIONS

Discrete layer solutions are presented for static responses of linearly elastic and thermoelastic multilayered composite plates, and have been found in good agreement with exact threedimensional solutions.

The proposed constrained discrete layer modelling is new, and allows satisfying exactly all the interface and boundary conditions for temperature, heat flux, displacements and stresses. Using thickness subdivisions along the thickness of the laminate from quadratic Lagrangian interpolations in the manner of quadratic one-dimensional finite elements through the thickness, interface and boundary conditions are exploited to reduce the number of independent generalized temperatures and independent generalized displacements. The result is that independent generalized temperatures and displacements are only those at the medium node of the subdivision. The boundary value problem is then solved from standard variational theorems (for thermal and mechanical analyses), assuming double Fourier series expansions for both temperature, displacement and loads, in respect to the surface coordinates. Adopting for example one subdivision per layer in a laminate having N layers, and quadratic Lagrangian interpolation along the thickness of this laminate, the constrained discrete layer modelling leads to N independent generalized temperatures, 3N independent generalized displacements, while the standard discrete-layer modelling requires 2N -1 independent generalized temperatures and 3(2N + 1) independent generalized displacements. This is the main advantage of the constrained discrete layer approach, in addition to the exact satisfaction of interface and boundary conditions. Finally, the constrained discrete layer modelling has been successfully extended to thermoelastic analysis of laminates, by computing temperature distributions followed by a thermomechanical response, using the same method of discretisation, and comparing it with exact thermoelastic calculations.

- 1 ≤

 1 ξ j ≤ 1 ξ j = (2z-(z j-1 +z j )) (z j -z j-1 )

FIG. 5 . 4 . 3 .

 543 FIG. 5. Thickness distributions of non-dimensionalized displacements and stresses extracted from Eq. (27) for the three-layer rectangular plate in Problem 2. (S = 4; L 2 /L 1 = 3; solid lines: elasticity;

FIG. 8 .FIG. 9 .

 89 FIG. 8. Thickness distributions of non-dimensionalized displacements and stresses extracted from Eq. (27) for the three-layer square plate in Problem 3. (S = 10; solid lines: elasticity; dashed lines:

FIG. 11 .

 11 FIG. 11. Transverse distributions for non-dimensionalized thermo-elastic displacements and stresses extracted from Eq. (26), for the laminate associated with Problem 5; S = 4.

TABLE 1

 1 Convergence of the CDLn model with respect to the number n = 1, 2, 3 of subdivisions in each layer. Calculations are extracted from Eq. (26). Problem 1: S = 4; cylindrical bending of a three-layer laminate

		Elasticity	CDL1	CDL2	CDL3
	ū1 (0)	0,9212	1,0010	0,9339	0,9259
	ū1 (1/3)	-0,2018	-0,3086	-0,2054	-0,2039
	ū1 (2/3)	0,2509	0,3621	0,2547	0,2531
	ū1 (1)	-0,9457	-1,0332	-0,9594	-0,9508
	ū3 (1/2)	2,8901	2,7380	2,8860	2,8893
	σ11 (0)	-18,087	-19,655	-18,336	-18,180
	σ11 (1/3)	4,0400	6,1219	4,1119	4,0806
	σ11 (2/3)	-4,7555	-6,9234	-4,8316	-4,7984
	σ11 (1)	18,819	20,536	19,087	18,919
	σ13 (1/2)	1,4319	1,5034	1,4435	1,4289

TABLE 2

 2

	Convergence of the CDLn model with respect to the number
	n = 1, 2, 3 of subdivisions in each layer. Calculations are
	extracted from Eq. (26). Problem 1: S = 10; cylindrical
		bending of a three-layer laminate
		Elasticity	CDL1	CDL2	CDL3
	ū1 (0)	9,37	9,49	9,39	9,38
	ū1 (1/3)	1,389	1,240	1,382	1,386
	ū1 (2/3)	-1,326	-1,178	-1,320	-1,324
	ū1 (1)	-9,35	-9,46	-9,36	-9,35
	ū3 (1/2)	0,933	0,929	0,933	0,933
	σ11 (0)	-73,61	-74,52	-73,71	-73,65
	σ11 (1/3)	-10,84	-9,69	-10,79	-10,82
	σ11 (2/3)	10,60	9,45	10,55	10,58
	σ11 (1)	73,66	74,57	73,76	73,69
	σ13 (1/2)	4,239	4,386	4,245	4,235

FIG. 3. Thickness distributions of non-dimensionalized displacements and stresses extracted from Eq. (26) for the laminate under cylindrical bending in Problem 1. (S = 4; solid lines: elasticity; dashed lines: present CDL3 model.) FIG. 4. Thickness distributions of non-dimensionalized displacements and stresses extracted from Eq. (26) for the laminate under cylindrical bending in Problem 1. (S = 10; solid lines: elasticity; dashed lines: present CDL3 model.) 4.

2. Problem 2. Simply Supported Laminated 0 • /90 • /0 • Symmetric Rectangular Plate under Doubly Sinusoidal Load Material

  data are same as those given for Problem 1

TABLE 3

 3 Convergence of the CDLn model with respect to the number

	Elasticity	CDL1	CDL2	CDL3

n = 1, 2, 3 of subdivisions in each layer. Calculations are extracted from Eq. (27). Problem 2: S = 4; sandwich rectangular plate L 2 /L 1 = 3

TABLE 4

 4 Convergence of the CDLn model with respect to the number n = 1, 2, 3 of subdivisions in each layer. Calculations are extracted from Eq. (27). Problem 2: S = 10; sandwich rectangular plate L 2 /L 1 = 3

	Elasticity	CDL1	CDL2	CDL3

TABLE 5

 5 Convergence of the CDLn model with respect to the number n = 1, 2, 3 of subdivisions in each layer.

	Calculations are
	extracted from Eq. (27). Problem 3: S = 4; sandwich
	square plate

. Problem 4. Cylindrical Bending of a Simply Supported Two-Layers 90 • /0 • Laminate under a Sinusoidal Load In

  this problem (see

	)	-0,0228	-0,0210	-0,0228	-0,0228
	ū3 (1/2)	2,006	1,900	2,003	2,005
	σ11 (0)	-0,755	-0,826	-0,766	-0,759
	σ11 (1/3)	0,190	0,279	0,193	0,192
	σ11 (2/3)	-0,230	-0,324	-0,233	-0,232
	σ11 (1)	0,801	0,881	0,813	0,805
	σ22 (0)	-0,079	-0,074	-0,079	-0,079
	σ22 (1/3)	-0,556	-0,597	-0,562	-0,559
	σ22 (2/3)	0,534	0,574	0,540	0,537
	σ22 (1)	0,095	0,091	0,095	0,095
	σ23 (1/2)	0,217	0,122	0,249	0,215
	σ13 (1/2)	0,256	0,267	0,254	0,255
	σ12 (0)	0,051	0,049	0,051	0,051
	σ12 (1/3)	0,0074	0,0064	0,0075	0,0075
	σ12 (2/3)	-0,0055	-0,0043	-0,0055	-0,0055
	σ12 (1)	-0,0511	-0,0499	-0,0513	-0,0512
	lished to accurately compute displacements and stresses for thick
	laminates.
	4.4

Elasticity

CDL1 CDL2 CDL3 ū1 (0) 0,00936 0,01028 0,00950 0,00941 ū1 (1/3) -0,00242 -0,00356 -0,00245 -0,00244 ū1 (2/3) 0,00317 0,00438 0,00321 0,00319 ū1 (1) -0,00969 -0,01073 -0,00985 -0,00975 ū2 (0) 0,0228 0,0210 0,0228 0,0228 ū2 (1/3) 0,0072 0,0077 0,0072 0,0072 ū2 (2/3) -0,0066 -0,0071 -0,0067 -0,0067 ū2 (1) -0,0228 -0,0210 -0,0228 -0,0228 ū3 (1/2) 2,006 1,900 2,003 2,005 σ11 (0) -0,755 -0,826 -0,766 -0,759 σ11 (1/3) 0,190 0,279 0,193 0,192 σ11 (2/3) -0,230 -0,324 -0,233 -0,232 σ11 (1) 0,801 0,881 0,813 0,805 σ22 (0) -0,079 -0,074 -0,079 -0,079 σ22 (1/3) -0,556 -0,597 -0,562 -0,559 σ22 (2/3) 0,534 0,574 0,540 0,537 σ22 (1) 0,095 0,091 0,095 0,095 σ23 (1/2) 0,217 0,122 0,249 0,215 σ13 (1/2) 0,256 0,267 0,254 0,255 σ12 (0) 0,051 0,049 0,051 0,051 σ12 (1/3) 0,0074 0,0064 0,0075 0,0075 σ12 (2/3) -0,0055 -0,0043 -0,0055 -0,0055 σ12 (1) -0,0511 -0,0499 -0,0513 -0,0512

APPENDIX T T

(rs)i is the line i (i = 1, 2, 3) in T (rs) Ĉ (j)T (rs)i3 =


Clément Pellegrini
email: pelleg@math.univ-lyon1.fr

Existence, uniqueness and approximation of stochastic Schrödinger equation: the diffusive case

Recent developments in quantum physics make heavy use of so-called "quantum trajectories". Mathematically, this theory gives rise to "stochastic Schrödinger equations", that is, pertubations of Schrödinger-type equations under the form of stochastic differential equations. But such equations are in general not of the usual type as considered in the litterature. They pose a serious problem in terms of: justifying the existence and uniqueness of a solution, justifying the physical pertinence of the equations.

In this article we concentrate on a particular case: the diffusive case, for a two-level system. We prove existence and uniqueness of the associated stochastic Schrödinger equation. We physically justify the equations by proving that they are continuous time limit of a concrete physical procedure for obtainig quantum trajectory.

Introduction

Belavkin equations (also called stochastic Schrödinger equations) are classical stochastic differential equations describing the evolution of an open quantum system undergoing a continuous quantum measurement. The solutions of such equations are called quantum trajectories and describe the time evolution of the state of the system. The random nature of the result of quantum measurement is at the origin of the stochastic character of the evolution.

The first rigorous description of a state undergoing a continuous measurement is due to Davies in [START_REF] Davies | Quantum theory of open systems[END_REF]. It describes in quantum optics, the behaviour of an atom from which we observe the photon emission. This is the so-called Resonnance Fluorescence experiment (see [START_REF] Kümmerer | An ergodic theorem for quantum counting processes[END_REF] and [START_REF] Bouten | Stochastic Schrödinger equations[END_REF]) .

In the literature, essentially two kinds of Belavkin equations are considered: they are either driven by a Brownian motion or by a counting process. But the kind of equations which are obtained this way are of non-usual type compared to the usual theory of stochastic differential equations. In particular there is no reference in physics nor in mathematics, where the existence and the uniqueness of the solution of such equations is discussed. Furthemore, the physical justification of the apparition of these equations requires in general quite heavy mathematical framework (Von-Neumann algebra, conditional expectation, filtering...). the high technology of such tools contrasts with the simplicity and the intuition of the physical model.

An approach to such equations, which is physically very intuitive, is the one of repeated quantum interactions. The setup is the following. The continuous measurement model is obtained as a limit of discrete models. This discrete model is a naive approach to the interaction of a simple system interacting with a field. The field is represented as a chain of independent copies of small pieces of environment. The simple system interacts, for a time interval h, with one piece of the environment. After that interaction an observable of the piece of environment is measured. The random result of the measurement induces a random new state for the small system. The small system then interacts again with another piece of the environment for a time interval h. A measurement of the same observable of this second copy is performed. And so on.

This experiment gives rise to a discrete evolution of the state of the small system, which is a Markov chain. The continuous time limit (h → 0) of this evolution should give rise to the quantum trajectories.

Repeated quantum interactions have been considered by Attal-Pautrat in [START_REF] Attal | From repeated to continuous quantum interactions[END_REF] and by Gough in [START_REF] Gough | Stochastic Schrödinger equations as limit of discrete filtering[END_REF]. The continuous limit of repeated quantum interactions is rigorously shown to converge to a quantum stochastic differential equation in [START_REF] Attal | From repeated to continuous quantum interactions[END_REF]. The setup of measuring an observable of the chain after each interaction is considered in [START_REF] Gough | Stochastic Schrödinger equations as limit of discrete filtering[END_REF], but the continuous limit, the existence and the uniqueness of the solutions are not all treated rigorously in this reference.

The aim of this article is to study the diffusive Belavkin equation, to show existence and uniqueness of the solution, to show its approximation by repeated quantum interaction models. The same results for the equation concerning the counting process are developed in another article [START_REF] Pellegrini | Existence, uniqueness and approximation of stochastic Schrödinger equation: the Poissonian case[END_REF]. This article is structured as follow: In section (1) we present the model of quantum measurement interaction. In particular we precise the mathematical way to construct the principle of quantum repeated interactions. We define then the sequence of states resulting to the measurement and we show the Markov property. To finish we establish a stochastic finite difference equation which appears as a discrete stochastic differential equation.

The section (2) is devoted to the continuous stochastic Belavkin equation. We deal with the problems of existence and uniqueness of solutions of such equations.

In section (3) we combine the two previous section. Thanks to the discrete model of section [START_REF]Open quantum systems[END_REF] we obtain an approximation sequence of states. We define the candidate to be the good approximation of the diffusive solution. To conclude we present the convergence theorem which justifies the diffusive model of Belavkin equation.

The discrete quantum trajectory

The main goal of this section is to describe the mathematical model of repeated quantum measurement. We present the general theory of indirect measurement which is necessary to obtain significant information. We describe the evolution of the small system undergoing successive measurement through the discrete quantum trajectory. We focus on probabilistic properties of the sequence representing the discrete quantum trajectory, in particular the Markov character.

Repeated quantum measurement

This section is devoted to lay out the mathematical framework of the naive physical model presented in the introduction.

To obtain information about the evolution of the state of a system, in quantum optics, the principle of indirect measurement is used. A small system interacts with a field (a photon stream for example) on which we perform a measurement and we are interested in the evolution of the small system. We are going to point out this method in the case of quantum repeated interaction.

The field is represented as a chain of independent copies of pieces of environment. Each copy is represented by a Hilbert space H and the small system is described by H 0 . The copies interact with the small system one after the other. The mathematical description of one interaction is the following.

The compound system is described by the tensor product H 0 ⊗ H and the interaction is characterized by a unitary evolution U acting on this Hilbert space. In the Schrödinger picture, if ρ denotes any state on the tensor product the evolution is given by:

ρ → U ρ U ⋆ .
Here the time of interaction is not specified. It will be important when considering approximations, all these questions are treated in [START_REF] Attal | From repeated to continuous quantum interactions[END_REF] and this is the subject of the section 3. For instance we consider that this time is fixed (equal to one for example) but we keep in mind that the unitary operator U depends on the interaction length time. After the above interaction, we consider a second copy of H which interacts with H 0 in the same way.

This sequence of interaction is described by the state space:

Γ = H 0 ⊗ k≥1 H k (1)
The countable tensor product k≥1 H k means the following. We consider that H 0 and H are finite dimensional Hilbert spaces. Let {X 0 , X 1 , . . . , X n } be a fixed orthonormal basis of H, the projector on X 0 : |X 0 X 0 | being the ground state (or vacuum state) of H (this is the bra-ket notation in mathematical physics see the remark below). The tensor product is taken with respect to X 0 (for all details concerning countable tensor product see [START_REF] Attal | From repeated to continuous quantum interactions[END_REF]).

Remark: A vector Y in a Hilbert space H is represented by the application |Y from C to H which acts with the following way |Y (λ) = |λY . The linear form on H are represented by the operators Z| which acts on the vector |Y by Z||Y = Z, Y where , denotes the scalar product of H. The unitary evolution describing the k-th interaction is given by U k which acts non trivially like U on H 0 ⊗ H k whereas it acts like the identity operator on the other copies. If ρ is a state on Γ, the effect of the k-th interaction is:

ρ → U k ρ U ⋆ k
Hence the result of the k first interactions is described by the operator V k on B(Γ) defined by the recursive formula:

V k+1 = U k+1 V k V 0 = I (2)
and its effect on a state is given by:

ρ → V k ρ V ⋆ k .
The above equation gives us the description of quantum repeated interaction. Let us move on to the principle of quantum repeated indirect measurement. The idea is to perform a measurement with respect of an observable of the field at each interaction. Let A = p j=1 λ j P j be any observable on H, then we consider its natural ampliation which defines an observable on Γ by:

A k := k j=0 I ⊗ p j=1 λ j P j ⊗ j≥k+1 I (3)
As a consequence if ρ is any state on Γ the probabilistic theory of a quantum measurement with respect to A k gives:

P [to observe λ j ] = T r[ ρ P k j ], j ∈ {1, . . . , p}
If we have observed the eigenvalue λ j the "projection" postulate called "wave packet reduction" imposes the new state to be

ρ j = P k j ρ P k j T r[ ρ P k j ]
.

Remark: This new state is then the reference state of our system. If we want to perform another measurement of the observable A k we obtain P [to observe λ j ] = 1. As a consequence a repeated interaction is necessary to obtain more than one information on H 0 . This above description gives the principle of measurement on the k-th copy. The quantum repeated measurement principle is the combination of the measurement principle and the repeated quantum interactions. Physically it means that each photons interacts with the atom and we perform a measurement after each interaction. After each procedure we have a new state given by the projection postulate: this is our discrete quantum trajectory.

The initial state on Γ is chosen to be

µ = ρ ⊗ j≥1 β j
where ρ is any state on H 0 and each β i = β is the reference state on H. We denote by µ k the state representing the new state after the k first interactions, that is

µ k = V k µV ⋆ k .
Let us now define the probabilistic framework and describe the effect of the successive measurements. We put Ω = {1, . . . , p} and on Ω N we define the cylinders of size k:

Λ i 1 ,...,i k = {ω ∈ Ω N /ω 1 = i 1 , . . . , ω k = i k }.
We endow Ω N with the σ-algebra generated by all these sets. This is the cylinder σ-algebra.

Remarking that for all j, the unitary operator U j commutes with all P k for all k < j, for {i 1 , . . . , i k } (corresponding to the index of eigenvalues). Thus we can define the following non normalized state :

μ(i 1 , . . . , i k ) = I ⊗ P i 1 ⊗ . . . ⊗ P i k ⊗ I . . . µ k I ⊗ P i 1 ⊗ . . . ⊗ P i k ⊗ I . . . = P k i k . . . P 1 i 1 µ k P 1 i 1 . . . P k i k .
So we can define a probability on Ω N , defined on the cylinders:

P [Λ i 1 ,...,i k ] = T r[μ(i 1 , . . . , i k )].
This probability satisfies the Kolmogorov consistency criterion, it defines then a probability on Ω N . Hence we define the following random sequence of states:

ρk (.) Ω N -→ B(Γ) ω -→ ρk (ω 1 . . . ω k ) = μ(ω 1 ...ω k ) T r[μ(ω 1 ...ω k )]
This random sequence of states is our discrete quantum trajectory and the operator ρk (i 1 , . . . , i k ) represents the state, if we have observed the results (i 1 , . . . , i k ) during the k first measurement. This fact is precised in the following proposition.

Proposition 1 Let (ρ k ) be the above random sequence of states we have for all ω ∈ Ω N :

ρk+1 (ω) = P k+1 ω k+1 U k+1 ρk (ω) U ⋆ k+1 P k+1 ω k+1 T r ρk (ω) U ⋆ k+1 P k+1 ω k+1 U k+1 .
This proposition is obvious but summarizes the quantum repeated measurement principle. The sequence ρk is the quantum trajectory rendering the effect of the successive measurements on Γ. The following theorem is an easy consequence of the previous proposition.

Theorem 1 The sequence (ρ n ) n is a Markov chain valued on the set of states of H 0 i≥1 H i It is described as follows:

P ρn+1 = µ/ρ n = θ n , . . . , ρ0 = θ 0 = P ρn+1 = µ/ρ n = θ n
If ρn = θ n then ρn+1 takes one of the values:

P n+1 i (U n+1 (θ n ⊗ β)U ⋆ n+1 ) P n+1 i T r (U n+1 θ n U ⋆ n+1 ) P n+1 i i = 1, . . . , p with probability T r (U n+1 θ n U ⋆ n+1 ) P n+1 i .
In quantum theory it was assumed that we do not have access to the field (because it is more complicated), we just have access to the small system. So the mathematical tools rendering this phenomenon is the partial trace operation given by the following theorem.

Definition-Theorem 1 If we have a state α on a tensor product H ⊗ K. There exists a unique state η on H which is characterized by the property:

∀X ∈ B(H) T r H [ η X ] = T r H⊗K [ α(X ⊗ I) ].
Hence to obtain the trajectory concerning the small system we have to take the partial trace on H 0 . Let E 0 denotes the partial trace on H 0 with respect to k≥1 H k . We then define a random sequence of states on H 0 . For all ω in Ω N we put:

ρ n (ω) = E 0 [ρ n (ω)]. (4)
This defines a sequence of state on H 0 which contains the "partial" information given by the measurement and we have the following theorem which is a consequence of theorem (1).

Theorem 2

The random sequence defined by formula (4) is a Markov chain with values in the set of states on H 0 . If ρ n = χ n then ρ n+1 takes one of the values:

E 0 I ⊗ P i U(χ n ⊗ β)U ⋆ I ⊗ P i T r[ U(χ n ⊗ β)U ⋆ I ⊗ P i ] i = 1 . . . p with probability T r [U(χ n ⊗ β)U ⋆ P i ].
Remark: Let us stress that

I⊗P i U (χn⊗β)U ⋆ I⊗P i T r[ U (χn⊗β)U ⋆ I⊗P i ]
is a state on H 0 ⊗ H, we have kept the notation E 0 to denote the partial trace on H 0 .

The next section is devoted to the case H 0 = H = C 2 which represents a two-level atom in contact with a photon stream. Because of physical consideration this particular case is often the central case in the literature concerning continuous measurement. The different results will be establish in this setting.

A two-level atom

In this section we want to establish a discrete quantum evolution equation for (ρ n ) which is a discrete equivalent of the Belavkin equation. In the previous section we have treated the general case, here we are going to work in 2 dimensionnal Hilbert spaces.

The main goal of this section is to obtain a formula of the following form:

ρ k+1 = f (ρ k , X k+1 ). ( 5
)
where (X k ) k is a sequence of random variables. In order to obtain such a formula we study how to obtain ρ k+1 through the measurement after the (k + 1)-th interaction when the initial state after k procedures is ρ k . The state ρ k can be namely considered as a initial state (according to the Markov property cf theorem ( 2)). Thus we consider a single interaction with a system (H, β) (actually this is the k + 1-th copy). Remember that each Hilbert space are C 2 . We consider an observable of the form A = λ 0 P 0 + λ 1 P 1 and the unitary operator describing the interaction is a unitary 4 × 4 matrix. We consider it as an operator on H 0 :

U = L 00 L 01 L 10 L 11
where each L ij are operators on H 0 . In order to compute the state given by the projection postulate we choose a suitable basis.

If (X 0 = Ω, X 1 = X) is an orthonormal basis of C 2 , for H 0 ⊗ H we consider the following basis Ω ⊗ Ω, X ⊗ Ω, Ω ⊗ X, X ⊗ X.
This basis allows us to consider the above way of writing for U. For β we choose:

β = |Ω Ω|
As a consequence, the state after the interaction is:

µ k+1 = U(ρ k ⊗ β)U ⋆ = L 00 ρ k L ⋆ 00 L 00 ρ k L ⋆ 10 L 10 ρ k L ⋆ 00 L 10 ρ k L ⋆ 10 . (6)
We apply the indirect quantum measurement principle. For the two possible results of the measurement we put:

L 0 (ρ k ) = E 0 [I ⊗ P 0 (µ k+1 )I ⊗ P 0 ] (7) L 1 (ρ k ) = E 0 [I ⊗ P 1 (µ k+1 )I ⊗ P 1 ]. ( 8
)
Thanks to the partial trace these are operators on H 0 . We denote the two probability by

p k+1 = T r[L 0 (ρ k )] and q k+1 = T r[L 1 (ρ k )].
The non normalized state: L 0 (ρ k ) appears with probability p k+1 and L 1 (ρ k ) with probability q k+1 .

Thanks to this two probabilities we can define a random variable ν k+1 on {0, 1} by: ν k+1 (0) = 0 with probability p k+1 ν k+1 (1) = 1 with probability q k+1

As a consequence we can describe the state on H 0 with the following equation. We have for all ω ∈ Ω N :

ρ k+1 (ω) = L 0 (ρ k (ω)) p k+1 (ω) (1 -ν k+1 (ω)) + L 1 (ρ k (ω)) q k+1 (ω) ν k+1 (ω) (9)
In order to obtain the final discrete quantum evolution equation we consider the centered and normalized random variable:

X k+1 = ν k+1 -q k+1 √ q k+1 p k+1 .
We define the associated filtration on {0, 1} N :

F k = σ(X i , i ≤ k).
So by construction we have E[X k+1 /F k ] = 0 and E[X 2 k+1 /F k ] = 1. Thus we can write the discrete evolution equation for our quantum trajectory.

ρ k+1 = L 0 (ρ k ) + L 1 (ρ k ) + [- q k+1 p k+1 L 0 (ρ k ) + p k+1 q k+1 L 1 (ρ k )]X k+1 . (10)
The above equation can be considered in a general way and the unique solution starting from ρ 0 is our quantum trajectory. Here the time interaction is chosen arbitrarly to be one. In Section 3 we are going to consider this equation with a interaction time h which is supposed later to go to zero. From the physical point of view, we perform a sequence of measurements spaced by a time h. In section 3 we present the link between the different models and the convergence results. Before treating these questions we shall focus on the equations which are supposed to model the continuous measurement, this is the content of the next section.

Belavkin equation

It is commonly assumed that the evolution of a system undergoing a continuous measurement is described by stochastic differential equation. High technical tools like continuous quantum filtering are usually necessary to obtain rigorous results about existence and uniqueness. This heavy machinery is refered to fine properties of Von-Neumann algebra and needs important background. Otherwise heuristic but more intuitive rules can be used to obtain such description in a non rigorous way (see [START_REF] Bouten | Stochastic Schrödinger equations[END_REF]).

The framework of our subject is the following. Consider a two-level atom (the small system) describing by C 2 and any state ρ in interaction with an environment (classically described by a Fock space endowed with a reference state). The time evolution is given by a unitary-process (U t ) which satisfies a quantum Langevin equation (cf [START_REF] Parthasarathy | An introduction to quantum stochastic calculus[END_REF]). Without measurement the evolution of the small system is given by a norm continuous semigroup {T t } t≥0 i.e ρ t = T t (ρ). The Linblad generator of (T t ) is denoted by L and we have the master equation:

dρ t dt = L(ρ t ) = -i[H, ρ t ] - 1 2 {CC ⋆ , ρ t } + Cρ t C ⋆
where C is any operator and H is the Hamiltonian of the atom.

In the theory of time continuous measurement L is decomposed as the sum of L + J where J represents the instantaneous state change taking place when detecting a photon, and L describes the smooth state variation in between these instants. These operators are defined by

L(ρ) = -i[H, ρ] - 1 2 {CC ⋆ , ρ} J (ρ) = CρC ⋆ .
Thanks to the works of Davies in [START_REF] Davies | Quantum theory of open systems[END_REF] which describes the evolution of a state during a continuous measurement. We can obtain with more or less rigorous argument two different equations whose solutions are called quantum trajectories:

The diffusive equation (Homodyne detection experiment) is:

dρ t = L(ρ t )dt + [ρ t C ⋆ + Cρ t -T r (ρ t (C + C ⋆ )) ρ t ]dW t
where W t design a one-dimensional brownian motion.

The jump equation (Resonnance fluorescence experiment) is:

dρ t = L(ρ t )dt + J (ρ t ) T r[J (ρ t )] -ρ t (d Ñt -T r[J (ρ t )]dt)
where Ñt is assumed to be a counting process with stochastic intensity

t 0 T r[J (ρ s )]ds.
The main goal of this article is the justification of the Belavkin model through the previous discrete description. In the last section the continuous equation appears as a limit of the discrete equation. Such convergence theorem need a deep study of the continuous stochastic differential equation. For instance we don't speak about validity, we are going to study such equations in a general way.

In this paper we consider only the diffusive case, the jump-equation and all convergence theorems refering to this case are treated in details in [START_REF] Pellegrini | Existence, uniqueness and approximation of stochastic Schrödinger equation: the Poissonian case[END_REF] with different techniques.

Let ρ 0 be any state, the aim is to show existence and uniqueness for the stochastic differential equation:

ρ t = ρ 0 + t 0 L(ρ s )ds + t 0 [ρ s C ⋆ + Cρ s -T r[(ρ s (C + C ⋆ )) ρ s ]dW s . (11)
Classical theorems concerning such equations can not be applied because the coefficients are not Lipschitz. A finer study is necessary to come to the result. From a mathematicalphysical point of view if there exists a solution we must check that the solution-process is valued on set of states. The problem of existence and the fact that the solution must be a state valued process are actually not independent. Before to answer the general problem we investigate the question of purity.

An important feature of the differential equation is that it preserves the property to be a pure state. In quantum theory a pure state is a one dimensionnal projector. Indeed if the initial state is pure and if there is a solution, the solution-process is valued on the set of pure state. This idea is mentioned in the following proposition and will be resumed in the final theorem: Proposition 2 Let (W t ) be a standard Brownian motion on (Ω, F , F t , P ) and let |ψ 0 be any vector in C 2 of norm one. Let ν t = 1 2 ψ t , (C + C ⋆ )ψ t , if the following stochastic equation:

d|ψ t = (C -ν t I)|ψ t dW t + -iH - 1 2 C ⋆ C -2ν t C + ν 2 t I |ψ t dt (12)
admits a solution (|ψ t ) then almost surely for all t ψ t = 1 and the process (|ψ t ψ t |) is a solution of the diffusive Belavkin equation [START_REF] Kurtz | Weak limit theorems for stochastic integrals and stochastic differential equations[END_REF].

Proof: Let |ψ 0 be any vector in C 2 and let (|ψ t ) be the solution of [START_REF] Kurtz | Wong-Zakai corrections, random evolutions, and simulation schemes for SDEs[END_REF]. We can compute d ψ t 2 . Thanks to Ito formulas: dW t dW t = dt and dW t dt = 0 and the fact that H is self-adjoint a straightforward computation shows that:

d ψ t 2 = d ψ t , ψ t = dψ t , ψ t + ψ t , dψ t + dψ t , dψ t = (C -ν t I)ψ t , ψ t dW t + (-iH - 1 2 (C ⋆ C -2ν t C + ν 2 t I)ψ t , ψ t dt + ψ t , (C -ν t I)ψ t dW t + ψ t , (-iH - 1 2 (C ⋆ C -2ν t C + ν 2 t I)ψ t dt + (C -ν t I)ψ t , (C -ν t I)ψ t dt = (2ν t -2ν t ψ t , ψ t )dW t
If ψ 0 2 = 1 it implies that almost surely for all t > 0:

ψ t 2 = ψ 0 2 = 1.
We define for all t ≥ 0 ρ t = |ψ t ψ t |, thanks to the fact that ψ t = 1 we then have for all y ∈ C 2 : ρ t |y = ψ t , y |ψ t .

So we can compute d ρ t |y by the Ito formula:

d ρ t |y = dψ t , y |ψ t + ψ t , y d|ψ t + dψ t , y d|ψ t = (C -ν t )ψ t , y |ψ t dW t + (-iH - 1 2 (C ⋆ C -2ν t C + ν 2 t )ψ t , y |ψ t dt + ψ t , y (C -ν t )|ψ t dW t + ψ t , y (-iH - 1 2 (C ⋆ C -2ν t C + ν 2 t )|ψ t dt + (C -ν t )ψ t , y (C -ν t )|ψ t dt
Let us recognize the equation [START_REF] Kurtz | Weak limit theorems for stochastic integrals and stochastic differential equations[END_REF]. It was clear that

ν t = 1 2 T r[|ψ t ψ t |(C + C ⋆ )].
As a consequence the term in front of the Brownian motion becomes:

(C -ν t )ψ t , y |ψ t + ψ t , y (C -ν t )|ψ t = (C|ψ t ψ t | + |ψ t ψ t |C ⋆ -T r[|ψ t ψ t |(C + C ⋆ )]|ψ t ψ t |) |y
A similar computation show that the term in front of dt is:

L(|ψ t ψ t |)|y
Hence we recognize the Belavkin equation ( 9) and the proposition is proved.

We can formulate the theorem concerning the existence and the uniqueness of a solution of ( 11):

Theorem 3 Let (Ω, F , F t , P ) be a probabilistic space which supports a standard Brownian motion (W t ) and let ρ 0 be any state, the stochastic differential equation:

ρ t = ρ 0 + t 0 L(ρ s )ds + t 0 [ρ s C ⋆ + Cρ s -T r[(ρ s (C + C ⋆ )) ρ s ]dW s
admits a unique solution (ρ t ) valued on the set of states and defined for all t ∈ [0, ∞[. Furthemore if the initial condition is a pure state the solution takes value in the set of pure states. The stochastic differential equation for a wave function is given by:

d|ψ t = (C -ν t )|ψ t dW t + -iH - 1 2 C ⋆ C -2ν t C + ν 2 t |ψ t dt where ν t = 1 2 ψ t , (C + C ⋆ )ψ t .
Proof: As the coefficients are not Lipschitz we can not apply directly the usual theorem for SDE's. However the coefficients are C ∞ , so locally Lipschitzian and we can use a truncature method. Our equation is of the following form :

dρ t = L(ρ t )dt + Θ(ρ t )dW t ( 13
)
where Θ is

C ∞ and Θ(A) = AC ⋆ + CA -T r [A(C + C ⋆ )] A. Let k ∈ N, k > 1, we define the truncation function ϕ k from R to R defined by ϕ k (x) =    -k if x ≤ -k x if -k ≤ x ≤ k k if -k ≤ x ≤ k
For a matrix A = (a ij ) we define by extension φk (A) = (ϕ k (Re(a ij )) + iϕ k (Im(a ij ))). Thus Θ • φk is a Lipschitz function. Now we consider the truncated equation:

dρ k,t = L • φk (ρ k,t )dt + Θ • φk (ρ k,t )dW t .
The theorem concerning stochastic differential equation driven by a drift term and a Brownian motion can be applied and there exists a unique solution t → ρ k,t defined for all t. Besides the solution is continuous.

We define the stopping random time

T k = inf{t, ∃(ij)/|Re(a ij (ρ k,t ))| = k or |Im(a ij (ρ k,t ))| = k}.
As ρ 0 is a state, for a suitable norm we have ρ 0 ≤ 1. Thanks to continuity, if k is choosed large enough we have T k > 0 and for all t ≤ T k we have φk (ρ k,t ) = ρ k,t . Thus t → ρ k,t is the unique solution of equation ( 11) (whithout truncation) on [0, T k ]. The classical method to solve an equation with non Lipschitz coefficients is to put T = lim k T k and to show that T = ∞.

In addition to the proof of existence of a solution we must prove that the process is valued on set of states. So when we deal with any state ν we have ν ≤ 1 and we have |ν(ij)| ≤ 1. Hence if we prove that on [0, T 2 ], the process (ρ 2,t ) is valued on set of states it proves that T 2 = ∞ a.s and then we have proved that there exists a unique solution valued on set of states. Let us prove this fact.

In the proof of the existence and uniqueness of a solution in the case of Cauchy-Lipschitz coefficients we create a sequence which converges to the solution:

ρ n+1 (t) = ρ n (0) + t 0 L • φk (ρ n (s))ds + t 0 Θ • φk (ρ n (s))dW s ρ 0 (t) = ρ (14)
With the right definition of Θ and L if ρ 0 is a state it is easy to see that this sequence is self-adjoint with trace one. This conditions are closed and at the limit the process is self-adjoint with trace one. The condition of positivity does not result of this sequence. We introduce the random time:

T 0 = inf{t ≤ T 2 /∃X ∈ C 2 / X, ρ 2,t X = 0} (15)
We have X, ρ 0 X ≥ 0 for all X, so by continuity we have X, ρ 2,t X ≥ 0 on [0, T 0 ] which implies that ρ 2,t is a state for all t ≤ T 0 . If T 0 = T 2 a.s the result is proved. Otherwise we have T 0 < T 2 , hence by continuity there exists X such that X, ρ 2,T 0 X = 0 and for all Y Y, ρ 2,T 0 Y ≥ 0. It implies that ρ 2,T 0 is a pure state because we work in dimension 2. Let us denote by ψ T 0 a vector of norm one such that ρ 2,T 0 = |ψ T 0 ψ T 0 |. We consider the equation:

d|ψ t = (C -ν t )|ψ t dW t + -iH - 1 2 C ⋆ C -2ν t C + ν 2 t |ψ t dt
with ψ T 0 as initial condition. The problem of existence and uniqueness for this equation is solved by a truncation method too. The fact that if we have a solution it is of norm one 1 shows that the solution obtained by truncation (defined for all t) is actually the solution of [START_REF] Kurtz | Wong-Zakai corrections, random evolutions, and simulation schemes for SDEs[END_REF]. The proposition (2) and the uniqueness of ρ 2,t on [T 0 , T 2 ] shows that the solution:

|ψ t = |ψ T 0 + t T 0 (C -ν s )|ψ s dW s + -iH - 1 2 C ⋆ C -2ν s C + ν 2 s |ψ s ds defines a process (|ψ t ψ t |) equals to ρ 2,t on [T 0 , T 2 ]
. Hence the process obtained by truncation is valued on set of states and the result is proved.

If the initial condition is a state the Belavkin equation ( 11) admits a unique solution valued on set of states. Moreover if there exist t 0 such that ρ t 0 is a pure state, for all t ≥ t 0 , the state ρ t is pure too.

Change of measure

Usually the stochastic equation appearing in the litterature is of the following form:

ρ t = ρ 0 + t 0 L(ρ s )ds + t 0 [ρ s C ⋆ + Cρ s -T r[ρ s (C + C ⋆ )]] d Ws ( 16
)
where:

Wt = W t - t 0 T r[ρ s (C + C ⋆ )]ds (17)
In the last section we will show a convergence of a discrete process to the solution of equation [START_REF] Kurtz | Weak limit theorems for stochastic integrals and stochastic differential equations[END_REF], the link between the two different equation is given by the Girsanov's theorem (see [START_REF] Protter | Stochastic integration and differential equations[END_REF] for a good introduction):

Theorem 4 Let (W t ) be a standard brownian motion on (Ω, F , F t , P ) and let H be a càdlàg process. Let:

X t = t 0 H s ds + W t (18)
and define a new probability by dQ dP = exp -

T 0 H s dW s -1 2
T 0 H 2 s ds for some T > 0. Hence under Q, the process (X t ) is a standard brownian motion for 0 ≤ t ≤ T .

The link between the two equation is then obvious. Let (ρ t ) be the solution of equation [START_REF] Kurtz | Weak limit theorems for stochastic integrals and stochastic differential equations[END_REF] given by the theorem (3) on (Ω, F , F t , P ). For some T > 0 we define the probability Q by:

dQ dt = exp T 0 T r[ρ t (C + C ⋆ )]dW s - 1 2 T 0 T r[ρ t (C + C ⋆ )] 2 ds (19)
The above theorem claim that Wt is the a standard Brownian motion under Q for 0 ≤ t ≤ T . As a consequence in the next section (concerning approximations and convergence theorems), we are going to consider (W t ) when we deal with a standard Brownian motion but we keep in mind that it is equivalent with the process which appears in physical literature.

Approximation and convergence

In this section we will construct a discrete quantum trajectory which converges to the solution of the diffusive Belavkin equation. The guiding idea is the fact that the discrete quantum trajectory satisfies a discrete stochastic equation which is of the same form as the continuous one.

The discrete approximation

At the beginning of the section 2 we have announced that there are essentially two kinds of Belavkin equation. We are going to see how the diffusive case appears from discrete process.

In section 1 we have obtained the following equation:

ρ k+1 = L 0 (ρ k ) + L 1 (ρ k ) + - q k+1 p k+1 L 0 (ρ k ) + p k+1 q k+1 L 1 (ρ k ) X k+1 (20)
Hence we have:

ρ k+1 -ρ 0 = k i=0 [ρ i+1 -ρ i ] = k i=0 [L 0 (ρ i )) + L 1 (ρ i ) -ρ i ] + k i=0 - q i+1 p i+1 L 0 (ρ i ) + p i+1 q i+1 L 1 (ρ i ) X i+1 (21)
The discrete process ρ k appears as the solution of a kind of a discrete time stochastic differential equation. This idea is going to be developed in order to obtain an approximation of the solution of [START_REF] Kurtz | Weak limit theorems for stochastic integrals and stochastic differential equations[END_REF].

We want now to introduce a discretisation of time. Consider a partition of [0, T ] in subintervals of equal size 1 n . The time of interaction is supposed now to be 1 n , the dynamic laws concerning the evolution of an open quantum system imposed that the unitary operator of evolution depends on the time interaction. We then have:

U(n) = L 00 (n) L 01 (n) L 10 (n) L 11 (n) .
The work of Attal-Pautrat [START_REF] Attal | From repeated to continuous quantum interactions[END_REF] has shown that the asymptotics of the coefficients L ij (n) must be properly rescaled in order to obtain a non-trivial limit. Indeed they have shown that V [nt] = U [nt] (n) . . . U 1 (n) which represents the discrete dynamic of quantum repeated interaction converges to an operator V t representing the continuous dynamic. Through the measurement theory we find a part of this result again. We choose:

L 00 (n) = I + 1 n -iH - 1 2 CC ⋆ + • 1 n (22) L 10 (n) = 1 √ n C + • 1 n (23)
The corresponding Hamiltonian H tot (n) describing the interaction during a time interval 1 n on H 0 ⊗ H is of the following form:

H tot (n) = H ⊗ I + I ⊗ 1 0 0 0 + 1 √ n C ⊗ 0 0 1 0 + C ⋆ ⊗ 0 1 0 0 + • 1 n
where H is the Hamiltonian of the small system and C is any operator.

With the time discretization we obtain:

ρ k+1 (n) = L 0 (n)(ρ k (n)) + L 1 (n)(ρ k (n)) + - q k+1 (n) p k+1 (n) L 0 (n)(ρ k (n)) + p k+1 (n) q k+1 (n) L 1 (n)(ρ k (n)) X k+1 (n)
The sequence of random variables (X k (n)) is defined through the two probabilities:

p k+1 = T r[L 0 (ρ k )] q k+1 = T r[L 1 (ρ k )].
Each L i depends on the measured observable: A = λ 0 P 0 + λ 1 P 1 . The aim of this section is to show the convergence of ρ [nt] (n) to the solution of diffusive Belavkin equation [START_REF] Kurtz | Weak limit theorems for stochastic integrals and stochastic differential equations[END_REF]. The Brownian noise will appear thanks to the sequence (X k ). By definition we have:

X k (n)(i) =    -q k+1 (n) p k+1 (n) with probability p k+1 (n) if i = 0 p k+1 (n) q k+1 (n) with probability q k+1 (n) if i = 1 (24)
As the probability and the operators L i depends on the observable, we are going to classify the observables in order to determine which ones give the diffusive nature.

If the observable is of the form A = λ 0 1 0 0 0 + λ 1 0 0 0 1 , we obtain the asymptotics for the probabilities:

p k+1 (n) = 1 - 1 n T r J (ρ k (n)) + • 1 n q k+1 (n) = 1 n T r J (ρ k (n)) + • 1 n
The discrete equation becomes:

ρ k+1 (n) -ρ k (n) = 1 n L(ρ k (n)) + •( 1 n ) + J (ρ k (n)) T r [J (ρ k (n)))] -ρ k (n) + •(1) q k+1 (n)p k+1 (n) X k+1 (n)
If the observable is non diagonal in the basis (Ω, X), we consider P 0 = p 00 p 01 p 10 p 11

and P 1 = q 00 q 01 q 10 q 11 we have:

p k+1 = p 00 + 1 √ n T r ρ k (p 01 C + p 10 C ⋆ ) + 1 n T r ρ k p 00 (C + C ⋆ ) + • 1 n q k+1 = q 00 + 1 √ n T r ρ k (q 01 C + q 10 C ⋆ ) + 1 n T r ρ k q 00 (C + C ⋆ ) + • 1 n
The discrete equation becomes then

ρ k+1 -ρ k = 1 n L(ρ k ) + • 1 n + e iθ Cρ k + e -iθ ρ k C ⋆ -T r[ρ k (e iθ C + e -iθ C ⋆ )] ρ k + •(1) 1 √ n X k+1
There appears a deterministic parameter θ which represents a kind of phase. This real parameter can be express thanks to the coefficients of the projector P i and can be introduced in the continuous equation. The results are the same. We then consider θ = 0 for the rest of the paper.

In [START_REF] Pellegrini | Existence, uniqueness and approximation of stochastic Schrödinger equation: the Poissonian case[END_REF] it is shown that the first case (when A is diagonal) gives the jump-Belavkin equation. In this section we shall show that the second case give rise to the diffusive Belavkin equation.

Convergence theorems

Before presenting the theorem concerning the convergence of the quantum trajectories, we expose a theorem concerning the average of the processes.

Theorem 5 Let (Ω, X) be any orthonormal basis of C 2 . For all non diagonal observable A, the deterministic function t → E[ρ [nt] (n)] converges in L ∞ ([0, T ]) to the function t → E[ρ t ] when n goes to infinity. That is:

sup 0<s<T E[ρ [ns] (n)] -E[ρ s ] -→ n→∞ 0.
Furthemore the function t → E[ρ t ] is the solution of the master equation:

d ν t = L(ν t )dt.
Proof: First of all we show the second part of theorem. We can consider the function t → E[ρ t ] because we have existence and uniqueness of the solution which is integrable (because ρ t is a state for all t). It was obvious that this deterministic function is process valued on set of states. As ρ 0 is a deterministic state we must show:

E[ρ t ] = ρ 0 + t 0 L(E[ρ s ])ds (25)
The state valued process (ρ t ) satisfies: According to the fact that the process (W t ) is a martingale and that the process (ρ t ) is predictable because continuous, the properties of stochastic integral with respect to a martingale give: We then have the integral form of the solution of the master equation and the second part is proved.

We will use a theorem of Kurtz and Protter (cf [START_REF] Kurtz | Wong-Zakai corrections, random evolutions, and simulation schemes for SDEs[END_REF]) to prove the convergence. If we consider a process define by: X n (t) = ρ 0 + ε n (t) + t 0 L(X n (s-)dV n (s) + t 0 Θ(X n (s-))dW n (s) Theorem 6 Suppose that W n is a martingale and V n is a finite variation process. Assume that for each t ≥ 0:

sup n E[[W n , W n ] t ] < ∞ sup n E[T t (V n )] < ∞
and that (W n , V n , ε n ) converges in distribution to (W, V, 0) where W is a standard brownian motion and V (t) = t for all t.

Suppose that X satisfies:

X t = X 0 + t 0 L(X s )ds + t 0 Θ(X s )dW s
and that the solution of this stochastic differential equation is unique. Then X n converges in distribution to X.

Recall that [X, X] is defined for a semi-martingale by the formula [X, X] t = X 2 tt 0 X s-dX s . The operation T t denotes the total variation of a finite variation process. For more details see [START_REF] Protter | Stochastic integration and differential equations[END_REF].

The convergence in distribution means the convergence in distribution for stochastic processes. We need another theorem to prove this type of convergence: Then M n converges in distribution to a standard Brownian motion. The conclusion is the same if we have:

lim n→∞ E [|[M n , M n ] t -t|] = 0.
Let us verify the different hypothesese to apply these theorems. We define a filtration for the process (W n (.)):

F n t = σ(X i , i ≤ [nt]).

ρ t = ρ 0 + t 0 L(ρ s )ds + t 0 [

 00 ρ s C ⋆ + Cρ s -T r (ρ s (C + C ⋆ )) ρ s ]dW s

E t 0 [

 0 ρ s C ⋆ + Cρ s -T r (ρ s (C + C ⋆ )) ρ s ]dW s = 0.Hence we have by linearity of L:E[ρ t ] = ρ 0 +

Theorem 7

 7 Let (M n ) be a sequence of martingales. Suppose thatlim n→∞ E[sup s≤t |M n (s) -M n (s-)|] = 0 and [M n , M n ] t -→ n→∞ t.

We shall now compare E[ρ [nt] (n)] with E[ρ t ] in order to obtain the convergence result. Like in the continuous case, the martingale argument is replaced by the fact that the process (X k ) is centered. Remember that we have:

As a consequence, considering k = [nt] and taking expectation in the discrete equation we have:

This is a kind of Euler scheme and we can conclude with a discrete Gronwall lemma to have the convergence:

The average of the discrete process is then an approximation of the average of ρ t . In [START_REF] Attal | From repeated to continuous quantum interactions[END_REF] this result was shown in the case of repeated interaction whithout measurement, it is a consequence of the asymptotics of the unitary-operator coefficients, so it justifies the choice of the coefficients of U(n). We are going to prove a similar result for the processes.

The discrete process which is the candidate to converge to the diffusive quantum trajectory satisfies for k = [nt]

Thanks to this equation we can define the processes:

By observing that this four processes are piecewise constant we can write the process (ρ n (t)) t≥0 like a solution of a stochastic differential equation in the following way:

) is a martingale. The process (W n (.)) converges to a standard Brownian motion W . when n goes to infinity and

Furthemore, we have the convergence in distribution for the process (W n , V n , ε n ) to (W, V, 0) when n goes to infinity.

Proof: Thanks to the definition of the random variable X k , we have

s ] = 0 for t > s. Thus if t s we have the martingale property:

By definition of [Y, Y ] for a stochastic process we have

Hence we have sup

Let us prove the convergence of (W n ) to (W . ). According to the theorem [START_REF] Jacod | The approximate Euler method for Lévy driven stochastic differential equations[END_REF] we must prove that:

It is a convergence in L 1 . We are going to prove a convergence in L 2 . In order to show the convergence in L 2 we will use the following property:

)] = 0. Thus we have:

Thanks to the fact that p 00 and q 00 are not equal to zero (because the observable is not diagonal) E[(X 2 i -1) 2 ] is bounded uniformly in i so we have:

Like [nt] n -→ t in L 2 we have the desired convergence. The convergence in distribution of (W n ) and (V n ) implies the convergence of (ε n ) to 0. This property is the essential point in the theorem of Kurtz an Protter. Thus we can express the final theorem.

Theorem 8 Let (Ω, X) be any orthonormal basis of C 2 and A be any observable non diagonal (in this basis). Let ρ be any state on C 2 .

Let (ρ [nt] (n)) be the discrete quantum trajectory obtained from the quantum repeated measurement principle with respect to A. The process (ρ [nt] (n)) then satisfies:

Let (ρ t ) be the solution of the diffusive Belavkin equation (11) which satisfies:

Thus we have the following convergence in distribution for stochastic process:

Proof: It is a simply compilation of the two theorems and the property for (W n (.)). To conclude we use the existence and uniqueness property proved in the second section (see theorem (3)).