François Hild

Stéphane Roux

Comparison of local and global approaches to digital image correlation

Keywords: 4-noded quadrilaterals, Displacement, Full-eld measurements, Resolution, Standard uncertainty, Strain

Local and global approaches to digital image correlation are compared when the displacement interpolation is based upon bilinear shape functions (i.e., with fournode quadrilaterals). The resolution in terms of displacements and strains associated with both techniques are evaluated a priori and validated a posteriori by using series of images of real experiments. It is shown that global approaches generally out-perform a local approach.

Introduction

Digital image correlation (DIC) is a popular technique to measure 2D and 3D displacement elds [START_REF] Bergonnier | Digital image correlation used for mechanical tests on crimped glass wool samples[END_REF]. In its 2D version, DIC consists in registering two images by measuring the displacement eld that enables for the best match. Generally, planar surfaces normal to the optical axis are observed to minimize perspective distortions, and non-rigid image registration is sought for. A random pattern is often applied on the observed Laboratoire de Mécanique et Technologie (LMT-Cachan) ENS Cachan / CNRS / UPMC / PRES UniverSud Paris 61 Avenue du Président Wilson, F-94235 Cachan Cedex, France ⋆ Corresponding author: hild@lmt.ens-cachan.fr surface, or the natural surface texture and markings are exploited to measure the displacement eld based on the assumption that the surface texture passively follows the displacement of the analyzed solid. DIC, which is used in solid mechanics, is equivalent to Particle Image Velocimetry (PIV) utilized in uid mechanics. DIC [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF] was initiated a few years later than PIV [812]. However, because the sought displacement and strain resolutions were smaller than the velocity resolutions, its development was slower; algorithmic challenges being more dicult to address. However, the technique has reached centipixel standard resolutions and uncertainties [1316], and millipixel biases or less [START_REF] Schreier | Systematic errors in digital image correlation caused by intensity interpolation[END_REF][START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF][START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF][START_REF] Wang | Quantitative Error Assessment in Pattern Matching: Eects of Intensity Pattern Noise, Interpolation, Strain and Image Contrast on Motion Measurements[END_REF] that are compatible with a large variety of cases encountered in solid mechanics [START_REF] Bergonnier | Digital image correlation used for mechanical tests on crimped glass wool samples[END_REF].

In its original developments, the registration was performed on zones of interest (ZOIs) that are small windows of the considered region of interest. This type of approach will be referred to as local. The early developments were based on pure translations that were sought either in the physical space [START_REF] Leclerc | Digital Volume Correlation: What are the limits to the spatial resolution?[END_REF] or in Fourier space [1820].

The kinematic hypotheses were changed later on to account for warping of the ZOI and rst order interpolations were implemented [START_REF] Sutton | Application of an optimized digital correlation method to planar deformation analysis[END_REF][START_REF] Vacher | Bidimensional strain measurement using digital images[END_REF]. Higher degree interpolations are also used [START_REF] Bornert | Assessment of Digital Image Correlation measurement errors: Methodology and results[END_REF] but in all the cases, the only information that is kept is the mean displacement of each ZOI.

In computer vision, very early on, alternative techniques were implemented. They are based upon variational formulations (as one would use in solid mechanics when resorting to, say, nite element formulations [START_REF] Zienkievicz | The Finite Element Method[END_REF]). In the work of Horn and Schunck, a spatial regularization was introduced to reduce the displacement uctuations [START_REF] Horn | Determining optical ow[END_REF]. However, in cases in which discontinuities are expected, the previous approach is not appropriate FAH IA [START_REF] Mitiche | Computation and analysis of image motion: A synopsis of current problems and methods[END_REF]. Smoother penalizations based on robust statistics were implemented [START_REF] Black | Robust Incremental Optical Flow[END_REF][START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. Another way of regularizing the registration is to enforce directly, say, continuity or even more so that no choice is left concerning the weighting associated with the regularization term. Various displacements elds were studied, namely, linear [START_REF] Roux | Correlation Image Velocimetry: A Spectral Approach[END_REF], Fourier expansions [START_REF] Roux | Correlation Image Velocimetry: A Spectral Approach[END_REF][START_REF] Wagne | Spectral Approach to Displacement Evaluation From Image Analysis[END_REF], B-spline [START_REF] Cheng | Full-eld Speckle Pattern Image Correlation with B-Spline Deformation Function[END_REF][START_REF] Hild | Displacement Measurement Technique for Beam Kinematics[END_REF], elastic solutions (closed-form [31, ! 32] and numerical [START_REF] Leclerc | Integrated digital image correlation for the identication of mechanical properties[END_REF][START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis fractured samples[END_REF]). In all these cases, because the registration is performed over the whole region of interest (ROI), it will be referred to as global approach.

It is also possible to use the same interpolations as those used in local approaches.

In that case, there is a need for discretization schemes and the nite element method is the most suited one. In computer vision, this type of approach was developed in the 1990s [START_REF] Simoncelli | Bayesian Multi-Scale Dierential Optical Flow[END_REF]. Finite-element based global approaches were introduced in the mid-2000s in solid mechanics. Either the correlation product [START_REF] Sun | Finite-element formulation for a digital image correlation method[END_REF] or the sum of squared dierences [START_REF] Broggiato | Adaptive image correlation technique for full-eld strain measurement[END_REF][START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF] were used as the minimization quantity. 4-noded bilinear elements are the most logical choice for comparison purposes. They will be used herein in addition to mesh adaption to allow for better resolutions of nodal displacements in the whole ROI. It is worth noting that unstructured meshes could also be used [START_REF] Leclerc | Integrated digital image correlation for the identication of mechanical properties[END_REF]. However, there are very few, if any, local approaches that use such types of discretizations.

The aim of the paper is to compare local and global approaches. For local as for global approaches, dierent matching algorithms are used [START_REF] Bornert | Assessment of Digital Image Correlation measurement errors: Methodology and results[END_REF][START_REF] Broggiato | Adaptive image correlation technique for full-eld strain measurement[END_REF][START_REF] Sun | Finite-element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]. To make the comparison objective, the same algorithms should be used in terms of gray level interpolation, kinematic interpolation, and minimization procedure. It will be the case hereafter. The paper is organized as follows. First, the principles of local and global approaches will be presented when the sum of squared dierences is minimized to evaluate the sought displacements. The sensitivity to acquisition noise is discussed in terms of the uncertainty levels of nodal displacements, mean displacements and mean strains per ZOI or element. Theoretical predictions are compared with an a priori resolution analysis. A rst practical example deals with a series of pictures taken with a CCD camera during a biaxial test on a thin lm. A couple of pictures acquired in a scanning electron microscope (SEM) are analyzed as a second example.

Correlation procedure

There are various ways of registering 2D pictures. Among them, the most popular in local approaches is related to the correlation product [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF]. An alternative approach is given by the sum of squared dierences. The latter was used by the authors when implementing global approaches in the context of nite element discretizations [START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]3234]. It will be used for both techniques in the sequel.

Let us consider f the picture of the reference conguration, and g that in the deformed conguration. Those two images capture the speckle pattern of the solid surface that is assumed to be passively advected by the supporting solid. The gray level conservation at any pixel location x then reads

f (x) = g(x + u(x)) ( 1
)
where u is the displacement vector. Because the conservation law is never strictly satised due to acquisition noise and the number of unknowns exceed the number of equations, it is not possible to determine the displacement vector u from the sole pictures f and g. The correlation procedure has therefore to be written on a given domain that contains more than 1 pixel. Thus a correlation procedure aims, for instance, at minimizing the sum of squared dierences

T = Ω (f (x) -g[x + u(x)]) 2 dx (2)
over the considered domain Ω (i.e., a ZOI in a local approach, and a ROI in a global approach) in which the displacement eld is interpolated as

u(x) = n unψn(x) (3)
where ψn are (chosen) vector functions, and un the associated degrees of freedom. The measurement problem then consists in minimizing T with respect to the unknowns un.

Let us note that choosing nodal displacements is a convenient way to emphasize the similarities and dierences between local and global approaches. However, the traditional parametrization used in local approaches is equivalent but dierent. The center point displacement being the key quantity of interest it is singled out, and other degrees of freedom correspond to a rst order polynomial expansion about the center.

A Newton iterative procedure is followed to circumvent the nonlinear aspect of the minimization problem. Let u i denote the displacement at iteration i, and {u} i the vector containing all the unknown degrees of freedom. By assuming small increments du = u i+1u i of the solution, a Taylor expansion is used to linearize g(x

+ v(x)) ≈ g(x) + v(x) • ∇g(x) ≈ g(x) + v(x) • ∇f (x)
and then, ∂T /∂{u} i is recast in a matrixvector product as

∂T ∂{u} i = [M]{du} -{b} i = {0} (4) with Mmn = Ω [∇f (x) • ψm(x)][∇f (x) • ψn(x)]dx (5) b i m = Ω [f (x) -gi (x)][∇f (x) • ψm(x)]dx (6)
and

gi (x) = g(x + u i (x)) (7)
It is worth noting that matrix [M] is computed once for all, and only vector {b} is updated from one iteration to the next. This choice (other alternatives are possible [START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]) is only dictated by computational eciency. The iterations stop when the displacement corrections {du} reach a small level that is chosen by the user. In many situations, the initial value is {u} i = {0}, but there are other propositions [START_REF] Hild | Displacement Measurement Technique for Beam Kinematics[END_REF][START_REF] Besnard | Characterization of necking phenomena in high speed experiments by using a single camera[END_REF].

To evaluate the resolution of the developed technique, since the un-noised reference is not known, the noise is assigned to the deformed picture with a variance 2σ 2 , where σ 2 is the variance of Gaussian white noise [START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF]. Under this assumption, the covariance matrix [C] simply reduces to [START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF][START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF] [C] = 2σ 2 p 2 [K] [START_REF] Barker | Measuring uid velocities with speckle patterns[END_REF] The resolution of a measuring system is the `smallest change in a quantity being measured that causes a perceptible change in the corresponding indication' [START_REF]2007, International Vocabulary of Metrology -Basic and General Concepts and Associated Terms, VIM[END_REF].

$ with [K] = [M]

-1 , and p the physical size of one pixel. This result is very general.

Similarly, it can be shown that when the noise level remains small, the estimate of u is unbiased.

Local approach

In the present setting, a local approach to DIC consists in minimizing the sum of squared dierences

T l = ZOI (f (x) -g[x + u(x)]) 2 dx (9)
over the considered ZOI. At this level of generality, many choices can be made for the displacement in each ZOI. In the following, shape functions used in nite element formulations are chosen. It is worth noting that the connectivity between elements is not yet enforced in this part. Since pictures are sampled in pixels it is logical to use square or rectangular ZOIs. Only regular square ZOIs will be considered herein. Their size is denoted by ℓ. The separation between ZOIs can be chosen, namely, it can be equal to ℓ (i.e., contiguous ZOIs), less than ℓ (i.e., overlapping ZOIs), or greater than ℓ (i.e., separate ZOIs). To compare with nite element calculations, contiguous ZOIs will be used in the sequel. The simplest interpolation to consider is bilinear (i.e., 1, x, y, xy functions, where x and y denote the local coordinates of any point M (x, y) belonging to a given ZOI). Each component of the displacement eld is treated the same way (i.e., only scalar functions Nn(x) are considered) so that the displacement u e (x) in each ZOI Ωe reads u e (x) = ne n=1 α a e αn Nn(x)eα [START_REF] Grousson | Study of ow pattern in a uid by scattered laser light[END_REF] where ne denotes the number of nodes (here ne = 4), eα the unit vector associated with direction α = 1, 2, and a e αn the unknown nodal displacements. Elementary matrices [M e ] whose components read

M e αmβn = Ωe ∂αf (x)Nm(x)∂ β f (x)Nn(x)dx (11)
and elementary vectors {b e } b e αm = [f (x) -g(x)] ∂αf (x)Nm(x)dx [START_REF] Pickering | Speckle Laser in Fluid Flows: Signal Recovery with Two-Step Processing[END_REF] where ∂αf = ∇f • eα, and g the corrected image, are computed for each ZOI. Matrix

[M e ] is symmetric, positive, and denite (when invertible). Contrary to classical FE procedures, quadrature formulas are not used because of the irregularity of the picture texture. Instead, a pixel summation is implemented [START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]. For each ZOI, an 8 × 8 linear system is solved iteratively

[M e ]{da e } = {b e } i (13)
where vector {a e } collects all sought degrees of freedom a e αn . The resolution analysis of Section 2.1 is applied to the kinematics chosen in this part. For each degree of freedom a e αm , the standard variation σ a e αm reads

σ a e αm = √ 2σp K e αmαm (14)
with [K e ] = [M e ] -1 . This result shows that depending on the underlying texture, and more precisely its gradient, the resolution of the correlation technique is not necessarily uniform over the whole ROI.

Let us assume that the correlation length of the texture is greater than a few pixels and less than the ZOI size ℓ. In such a case, a mean-eld assumption is used, namely, the squared image gradient that appears in the expression of [M e ] (Equation ( 11)) is extracted from the integral and changed into its expectation value

∂αf ∂ β f = G 2 f δ αβ (15)
where G f is the square mean of the image gradient projected along any direction. Note that here the image texture is assumed to be isotropic. This mean-eld assumption is expected to hold for large ZOI sizes, but it may break down for small ones. Under such circumstances, sub-matrix M e αβ reads

[M e αβ ] ∝

G 2 f ℓ 2 36 δ αβ          4 2 2 1 2 4 1 2 2 1 4 2 1 2 2 4          (16)
& Consequently, the covariance matrix associated with the measured degrees of freedom reads

[C αβ ] ∝ 8σ 2 p 2 G 2 f ℓ 2 δ αβ          4 -2 -2 1 -2 4 1 -2 -2 1 4 -2 1 -2 -2 4          (17)
The covariance matrix for the degrees of freedom is simply related to the inverse of The largest eigenvalue corresponds to a uniform translation. The second eigenvalue, which is twice degenerate and three times smaller than the rst one, corresponds to a uniform gradient along x and y. Finally, the fourth eigenvalue is a eld whose nodal value is proportional to {1, -1, -1, 1}, and as it is orthogonal to the three other eigenmodes, it corresponds to a deformation mode with zero translation and zero mean strain. However, it is nine times smaller than the major one. Thus noise has a much more pronounced eect on the latter one.

An estimate σa of the standard displacement resolution of each degree of freedom is given by

σa = 4 √ 2σp G f ℓ (18)
For any quantity of interest, which is linear in the displacement ϕ = i κ i a i , or ϕ = {κ} t {a}, its standard uncertainty reads

σϕ = {κ} t [C]{κ} 1/2 (19)
where κ i are known coecients collected in vector {κ}. This general formula can be used to evaluate the mean displacement uncertainty, the strain uncertainty, but also holds for more complex linear forms such as stress intensity factors (based on interaction integrals [START_REF] Réthoré | Noise-robust Stress Intensity Factor Determination from Kinematic Field Measurements[END_REF][START_REF] Réthoré | Optimal and noise-robust extraction of Fracture Mechanics parameters from kinematic measurements[END_REF]), or simply mean strain evaluated over arbitrary shaped and sized domains.

' For any point x within a ZOI, the displacement can be evaluated from the dierent shape functions, see Equation [START_REF] Grousson | Study of ow pattern in a uid by scattered laser light[END_REF]. Consequently, Equation ( 19) can be used to map out the uncertainty of σu(x). The shape of the uncertainty function has a minimum at the centre and a maximum value at the mesh nodes. The map is shown in Figure 1(a).

Consequences of this observation are that if one single information is to be kept from the local analysis then the center value is the most reliable one. The displacement at the center of ZOI reads

u = 1 4 m=1,4 a e αm ( 20
)
and the corresponding standard deviation

σ u = √ 2σp m=1,4 n=1,4 K e αmαn = σa 4 (21)
To evaluate the strain eld, the components of the displacement gradient ∇ ⊗ u e is computed in each ZOI Ωe by using the gradient of the shape functions

(∇ ⊗ u e )(x) = ne n=1 α a e αn (∇Nn(x) ⊗ eα) (22)
The mean component per ZOI is calculated by using the divergence theorem

p 2 ℓ 2 ∇ ⊗ u e = ∂Ωe (u e ⊗ n)ds ( 23
)
where n is the outward normal of any point belonging to the boundary ∂Ωe characterized by its curvilinear abscissa s. Let uα,α and u α,β denote the two components of the displacement gradient associated with uα.

so that their corresponding standard deviations reads

σ uα,α = √ 2σ 2ℓ m=1,4 n=1,4 (-1) m+n M -1 mn σ u α,β = √ 2σ 2ℓ m=1,4 n=1,4 (-1) m+n M -1 nm ( 25
)
As a rst approximation, the standard uncertainty level for the mean gradient compo-

nents σ ǫ becomes σ ǫ = 2 √ 6σ G f ℓ 2 (26)
so that it is related to σ u and σa by

σ ǫ = 2 √ 3 σ u ℓp = √ 3 2 σa ℓp (27)
Usually, the mean displacement gradients are evaluated by interpolating the mean displacements of neighboring ZOIs. In the following analysis, a set of 2 × 2 ZOIs is considered. The mean displacement gradients become

uα,α = u α4 -u α3 + u α2 -u α1 2ℓp , u α,β = u α4 + u α3 -u α2 -u α1 2ℓp (28)
and the corresponding uncertainty level σ ǫ becomes

σ ǫ = √ 2σ G f ℓ 2 = σ u ℓp (29)
where • denotes averages taken over 2 × 2 ZOIs.

All these results will be analyzed by using a real picture in Section 3. Figure 2 summarizes the results when a local approach is used and 10 × 10 ZOIs are considered.

By construction, the local approach yields displacement and strain uncertainties (i.e., σa, σ u , σ ǫ ) that are independent of the location when the texture content (i.e., (∂αf ) 2 is uniformly distributed over the whole ROI) and provided the ZOIs do not overlap.

Global approach

In a global approach to DIC the sum of squared dierences is minimized when dened over the whole ROI

Tg = ROI (f (x) -g[x + u(x)]) 2 dx (30)
In the following, the same shape functions as those chosen in the local approach will be considered. A regular square mesh made of 4-noded elements with a bilinear interpolation of the displacements is constructed. The main dierence with the previous approach is that the continuity of the displacement eld between elements is satised. Each component of the displacement eld satises Equation [START_REF] Grousson | Study of ow pattern in a uid by scattered laser light[END_REF]. Matrix [M] (see Section 2.1) is assembled by considering elementary matrices [M e ] whose components are dened in Equation ( 11) and vector {b} is assembled by calculating elementary vectors {b e } expressed in Equation [START_REF] Pickering | Speckle Laser in Fluid Flows: Signal Recovery with Two-Step Processing[END_REF]. Matrix [M] is symmetric, positive, denite (when invertible), and sparse as in any FE computation. A single global problem is solved iteratively

[M]{du} = {b} i (31)
where vector {b} i is updated thanks to the picture corrections g, and vector {u} collects all sought degrees of freedom. Because of connectivities between elements, the maximum component of matrix [M] is multiplied by 4 for inner nodes when compared with a local approach. Further, edge nodes and corner nodes have smaller values compared with inner nodes since the connectivities are dierent.

As shown above, the spectrum of eigenvalues and eigenmodes allows for a simple interpretation of noise sensitivity. In a global approach, and with a torus topology (in order to avoid edge eects), matrix

with kx = nxπ/L, ky = nyπ/L, and 0 ≤ nx ≤ L 0 ≤ ny ≤ L. Figure 3 shows the magnitude of the eigenvalues as a function of the wavevector k. The largest eigenvalue,

ξmax = G 2 f ℓ 2 , corresponds to k = 0.
The corresponding eigenmode is thus a uniform translation. As the wavelength decreases, so does the eigenvalue, and for a large system size, a series expansion of ξ as a function of k provides the following expression

ξ = G 2 f ℓ 2 1 -(1/6)|k| 2 (33)
where only the modulus of the wavenumber matters. At the other end of the spectrum, the smallest eigenvalue corresponds to a checkerboard nodal eld. In that case, kx = ky = π, and hence

ξ min = (1/9)G 2 f ℓ 2 (34)
The condition number of matrix [M] is thus equal to 9, provided the mean-eld assumption is valid. A consequence of this analysis is that when the noise level is large (or the element size small), the nodal displacement eld exhibits a checkerboard type component that can easily be seen. Hopefully, those high frequency modes have little impact on estimates of strain or mean displacement.

The resolution analysis of Section 2.1 is applied to the nite element implementation [START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF][START_REF] Besnard | Finite-element displacement elds analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]. Each inner node belongs to four elements so that the shape function associated with a given inner node spans over 4 elements. Consequently, the variance is divided by 4, and the standard displacement resolution σa is divided by 2 when compared with the value of a single Q4 element. However, this result does not account for the additional requirement of continuity in the displacement eld (i.e., the global approach deals with matrix-vector system that is written on a level of the ROI). To estimate the additional gain associated with the continuity constraint, a global matrix is formed by assembling elementary matrices Mαα for, say, a mesh made of 10 × 10 Q4 elements (Figure 2). For inner nodes, it is found that

σa = √ 6σp G f ℓ (35)
The multiplicative factor induced by continuity is therefore equal to √ 3/2. Edge nodes are shared by two elements so that the standard displacement resolution σa is divided by √ 2 when compared with its value for a single Q4 element. Again, this result does not account for the continuity of the displacement eld. When the latter is accounted for, the following result is obtained

σa = 2 √ 3σp G f ℓ (36)
The multiplicative factor due to continuity is equal √ 3/2 as for the previous case. Last, corner nodes belong to one element only, and the standard displacement resolution σa ! is identical to that of a single Q4 element. By following the same procedure as above, it is found that

σa = 2 √ 6σp G f ℓ (37)
The multiplicative factor is equal to Because the number of connectivities of each node of the mesh is dierent (i.e., from 1 to 4), it leads to a standard uncertainty σa that is no longer uniform over the whole ROI. This is also the case of σ u and σ ǫ (Figure 2). For inner nodes, by using the same approach as above, the following approximations are obtained for σ u

√ 3/2 in
σ u = 0.63 √ 2σp G f ℓ (38) and σ ǫ σ ǫ = 1.67 √ 2σ G f ℓ 2 (39)
In both cases, when normalized by the corresponding value or σa, they lead to upper bounds.

In the case of σ u , it can be noted that the factor 2 when compared with a local approach is not found. An upper bound to the ratio is equal to 0.36 (it is reached for inner nodes). This eect can be understood by comparing the uncertainty maps of σu(x) for both approaches (Figure 1). The shape of the uncertainty map is similar " with a minimum at the center and a maximum value at the mesh nodes. However, it is observed that consistently with the value reported for the nodal and center values, the entire function for the global case lies below that of the local case. It is also observed that the gain is much more important at the nodes, and much smaller for the average value. However, for the global approach the same result holds albeit it may appear as less trivial. In particular, when the displacement eld is to be compared with another such eld only at discrete points it is much more advantageous to select the central values rather than the nodal ones, while preserving the same density of points, the uncertainty is reduced by more than a factor of 2.

Conversely, for σ ǫ there is a gain of 0.48. In that case, the bonus of continuity is equal to 4 %. However, let us note that the mean strain is never calculated in such a way for a local approach. To compare with a local approach, the displacement gradients are evaluated by interpolating the nodal displacements

and the corresponding uncertainty level σ ǫ reads

σ ǫ = √ 2σ √ 3G f ℓ 2 = σa 3ℓp (41)
for inner nodes. Even though nodal displacements are used, a multiplication factor of

1/ √

3 is obtained when the local and global approaches are compared. Had the mean displacements of each elements been used instead, the result of Equation ( 29) still applies (σ ǫ = σ u /ℓp), so that

σ ǫ = 0.63 √ 2σ G f ℓ 2 = σa 2.75ℓp (42)
which is about 8 % higher compared to the previous case.

The use of a uniform mesh leads to a heterogeneous distribution of standard displacement and strain uncertainties due to the heterogeneity of matrix [M] (see Figure 2). It is possible to make the distribution more uniform. Depending on the type # of studied quantity, the result will be dierent. In the following analysis, the standard uncertainty of the nodal displacements is chosen. As a rst approximation, it was shown that σa was 2 times larger for corner nodes and √ 2 larger for edge nodes when compared to the value of inner nodes. Since the standard displacement uncertainty is inversely proportional to the element size, it suces to multiply the size of the corner elements by a factor of 2, and edge elements then become rectangular. Figure 2 shows the eect on the global matrix for which corner and edge nodes have components that have the same order of magnitude than inner nodes.

Under this assumption, the uctuation of the G f ℓσa/ √ 2σ eld varies between 1.73 (≈ √ 3) and 1.82 (Figure 2). The eld of G f ℓσ u / √ 2σ then varies between 0.87 and 0.63, and the eld G f ℓ 2 σ ǫ / √ 2σ between 0.7 and 1.67. These uctuations are smaller than those observed with a regular mesh.

Analysis of a synthetic case

In this rst analysis, an a priori resolution analysis is performed. It consists in considering a given picture f and adding Gaussian white noise with dierent standard deviations to form picture g. These two pictures are correlated by DIC. In a rst part, a local analysis with Q4 ZOIs is considered. The aim of the study is to evaluate the measurement uncertainties associated with the measured degrees of freedom, the average displacement per ZOI. The same analysis is run for the rst displacement gradients. It was obtained by spraying a black and white paint. It was used to identify damage parameters by using a biaxial test on a vinylester matrix reinforced by a mat of Eglass bers [START_REF] Claire | A nite element formulation to identify damage elds: The equilibrium gap method[END_REF][START_REF] Roux | Digital Image Mechanical Identication (DIMI)[END_REF]. In the present case, a 512 × 512-pixel denition is considered. The histogram (Figure 4(b)) does not cover completely the whole 8-bit dynamic range. The correlation radius of the texture is of the order of 3 pixels. Values of σ ranging from 1 to 32 gray levels are considered in the sequel. In many practical situations, σ is of the order of 1 percent of the dynamic range for CCD or CMOS cameras. However, other imaging systems may be used (e.g., scanning electron microscopes or atomic $ force microscopes) for which the noise level is signicantly larger than that observed for conventional cameras (see Section 5).

The rst output of a correlation code to be checked is the residual eld η(x) = f (x) -g[x + u(x)] and its RMS η over the whole ROI. In the present case, the ZOI separation is equal to the ZOI size so that this measure counts only once each pixel belonging the ROI. Figure 5 shows the change of the RMS residual η normalized by the dynamic range ∆f = max ROI f -min ROI f as a function of the dimensionless noise level √ 2σ/∆f . A linear trend is shown with a slope of 1. When the noise level reaches very high values, there is a deviation as the number of ZOIs increases, thereby indicating a higher sensitivity to noise, which is to be expected. Otherwise, there is an almost 1 to 1 correspondence between the two quantities. The fact that the correlation residual is virtually identical to the noise level is an indication that the registration was successful. Furthermore, when the noise level is unknown, as will be shown below, it is a simple way of assessing it by acquiring a series of pictures, performing correlations to account for slight kinematic deviations, and assessing the RMS value of the residual, which is equal to √ 2σ as a rst approximation.

The result of Equation ( 18) is analyzed next. Figure 6(a) shows the change of σa with the ZOI size ℓ and for dierent values of σ. The uncertainty level is increasing as σ increases and ℓ decreases. This trend is described by Equation [START_REF] Chen | Digital Speckle-Displacement Measurement Using a Complex Spectrum Method[END_REF]. On a more quantitative way, when σa is divided by the noise level σ, it should fall onto a single curve that varies as 1/ℓ. The results of Figure 6(b) show that this dependence is captured. A deviation is observed when only one ZOI is considered (i.e., ℓ = 512 pixels) or for very small ZOIs. For large ZOIs, very few points are available, making the statistical content poor. Conversely, for small ZOI sizes, the mean eld approximation is less likely to apply as the ZOI size becomes closer to the correlation length of the texture. On a quantitative way, it is found that ℓσa∆f / √ 2σp ≈ 47 to be compared with 4∆f /G f = 50. This last result validates Equation [START_REF] Chen | Digital Speckle-Displacement Measurement Using a Complex Spectrum Method[END_REF]

since σa ≈ 4 √ 2σp/G f ℓ.
When the mean displacement per ZOI is computed, it is expected that its uncertainty level σ u decreases when compared with that of the measured degrees of freedom % σa (see Equation ( 21)). Figure 7 shows a plot of the ratio σ u /σa as a function of σa.

The average value of the ratio is found to be equal to 0.24 (to be compared with the theoretical value of 1/4). The dierent data are scattered because, for small ZOIs in particular, there are uctuations of the mean picture gradient, and the hypothesis of scale separation is less satised.

Equation ( 26) is now studied. The change of σ ǫ with the ZOI size ℓ and for dierent values of σ is shown in Figure 8(a). The uncertainty level again increases with σ and decreases with ℓ as expected from Equation [START_REF] Odobez | Robust multiresolution estimation of parametric motion models[END_REF]. When σ ǫ is normalized by the noise level σ it falls onto a single curve that varies as 1/ℓ 2 (see Figure 8(b)). A deviation is also present when only one ZOI is considered for the same reason as above. The interpolation of the previous results leads to ℓ 2 σ ǫ ∆f / √ 2σ ≈ 39 to be compared with

2 √ 3∆f /G f = 43. The value 2 √
3 is therefore a very good estimate of the actual ratio

σ ǫ / √ 2σℓ 2 G f .
When the mean strain per ZOI is computed, it is expected that its uncertainty level σ ǫ depends linearly on the standard uncertainty of the measured degrees of freedom σa (see Equation ( 27)) and is inversely proportional to ℓ. Figure 9 shows a plot of the ratio σ ǫ /σa as a function of σa. The average value of the ratio is found to be equal to 0.84, which is very close to the theoretical value ( √ 3/2 ≈ 0.87). All the results found in the present section validate the closed-form solutions derived in Section 2.2.

The performances of a global approach is compared with a local one in Figure 10 when the mean values of the standard uncertainties are considered for all the noise levels. As the ZOI or element sizes decrease, the standard uncertainties σa, σ u and σ ǫ of the global approach tend to the theoretical limits when compared with the values of the local approach. This is due to the fact that for small elements, the proportion of inner nodes becomes dominant. These results validate the a priori predictions of the standard resolution in terms of nodal displacement, mean displacement and mean strain per ZOI or element for local and global approaches.

The series of pictures analyzed in this section were taken before starting a biaxial test on supported thin W lms. The sample was positioned in the goniometer center of the DiAbs experimental station at SOLEIL synchrotron source [START_REF] Geandier | Development of a synchrotron biaxial tensile device for in-situ characterization of thin lms mechanical response[END_REF]. Strains using

x-ray diraction and digital image correlation can be measured in this new setup (Fig-

ure 11(a)) The two techniques led to very similar results not only in terms of mean strain levels but also concerning the standard strain uncertainties [START_REF] Djaziri | Combined synchrotron X-ray and image-correlation analyses of biaxially deformed W/Cu nanocomposite thin lms on Kapton[END_REF]. The optical microscope is composed of a telecentric lens mounted on a CCD camera (12-bit dynamic range, denition: 1392 × 1024 pixels). This setup allows for the minimization of out-of-plane eects since a stereoscopic equipment was not possible to implement in the experimental environment. Values of σ/∆ ranging from 1/256 to 1/32 ∆f are considered in the sequel. Figure 12 shows the result of the resolution analysis for a local approach. The same trends as those observed in Figure 5 are obtained. The predictions of the a priori resolution analysis will be compared to actual measurements. A series of ten pictures, in addition to the reference one, was acquired. The stepping motors of the testing machine were motionless. This is typical of an analysis that can be performed experimentally to assess the actual resolution of the measuring system [START_REF]2007, International Vocabulary of Metrology -Basic and General Concepts and Associated Terms, VIM[END_REF]. The direct DIC analysis of the 10 pictures provides the values of the RMS residual η in addition to the measurement uncertainties. In the present case, η/∆f varies between 0.008 and 0.025 (see shaded zone in Figure 12). From this information, it is concluded that σ varies between 0.5 % and 1.8 % of the dynamic range of the camera, with an average value of 1.0 %. Figure 13 shows the change of the standard displacement and strain uncertainties as a function of the ZOI size. The latter ones are compared with the levels predicted by the resolution ' analysis when using the bounds of η and the results of Figure 12. A good agreement is obtained, thereby validating the a priori resolution analysis.

Let us now compare the measurement uncertainties of local and global approaches.

Lower uncertainty levels are to be expected with a global approach. Figure 14(a) shows the change of the standard displacement uncertainty of the measured degrees of freedom for both approaches. In both cases, as expected, the uncertainties increase when the ZOI or element sizes decrease. All approaches converge toward the same value when only one ZOI or element is considered. When 2 × 2 elements are chosen, the two global approaches coincide since no mesh adaption was performed. As the number of elements increases, the eect of mesh adaption is visible and a gain is observed in terms of σa level for the same element size. The gain is even more important when compared with the results of the local approach.

Figure 14(b) shows the change of the standard uncertainty of the average displacement per ZOI (local approach) or element (global approach). Let us note that the mean displacement is generally not an information that is given when using a global approach. Generally, the values that are shown are the nodal displacements. Yet, for comparison purposes, the levels of these values are compared herein. As expected from the a priori analysis, the gain between local and global approaches is less important.

However, the global approaches still out-perform the local approach as the ZOI or element sizes decrease.

Figure 14(c) shows the change of the standard uncertainty of the mean strain per ZOI (local approach) or element (global approach). It is worth noting that this information is generally not given as an output of a local approach. However, for comparison purposes it is reported herein. As the number of ZOIs or elements increases, the ratio of 2 between σ ǫ of the global approach compared with that of a local approach is found. This is faster when mesh adaption is considered.

A couple of pictures taken 30 minutes apart in an SEM is nally analyzed. They have an 8-bit digitization. A ROI of 750 × 750 pixels is extracted. This case was chosen because of the high acquisition noise level generally associated with SEM pictures as will be shown hereafter. In the present case, the secondary electron mode is selected with an Everhart-Thornley detector (ETD). Because the sample surface was not textured enough, a grid was deposited by resorting to microphotolithography (Figure 15(a)).

This type of texture leads to a bimodal histogram (Figure 15(b)) corresponding to the grids (high gray levels) and the raw surface (low gray levels). Moreover, the whole dynamic range of the detector is not used in the present case, but only 80 %.

When the couple of pictures is analyzed, the RMS residual η/∆f is found to be of the order of 7 % for local and global approaches, which is signicantly higher than what was observed above. The a priori analysis is run on the reference picture of Figure 15(a) for values of σ/∆ ranging from 1/64 to 1/8. By using the results of Figure 16, it is concluded that the noise level associated with the SEM pictures is signicantly higher than for a conventional CCD or CMOS camera. This case is therefore interesting to check whether the theoretical predictions survive when the noise level of pictures becomes higher.

In Figure 17 For very large ZOIs, the local approach out-performs the global approach. This may occur since it was shown that the ratio σ ǫ /σa is lower for a local approach (i.e., 1/4) compared with a global approach (i.e., 0.36 at most for inner nodes). However, the trend is inverted for smaller values of ℓ for which the global approaches out-perform again the local one.

Figure 17(c) shows the change of the standard uncertainty of the mean strain per ZOI or element as a function of ℓ. In the present case, there is a clear benet from the global approach. Very early, the ratio of 2 between the global and local approaches is reached. Because of the high noise level, there is an additional bonus in using global approaches.

Summary

A comparison between local and global approaches to DIC was performed by using 4noded zones of interest (local approach) and 4-noded elements (global approach). The same minimization procedure was used in both approaches so that algorithmic issues are not the cause of dierences. In both cases, the sum of squared dierences was minimized by using a modied Newton scheme. It was shown that the global approach out-performs the local one thanks to the continuity requirements and the fact that shape functions span over four elements for inner nodes. It is worth noting that the results derived herein with Q4 shape functions can be generalized to 8-noded cubes with trilinear shape functions [START_REF] Leclerc | Digital Volume Correlation: What are the limits to the spatial resolution?[END_REF] when performing Digital Volume Correlation [START_REF] Roux | Three dimensional image correlation from X-Ray computed tomography of solid foam[END_REF].

General results concerning the resolution in terms of displacement and strains were derived. Under simplifying assumptions, closed-form expressions were obtained. They were checked against articial and practical cases, and validated in all the examples analyzed herein. In particular, a mesh adaption strategy was proposed to yield nodal displacement uncertainties that are uniform for a global approach, irrespective of the fact that nodes are located inside the region of interest, or on its edges or at corners.

The gain in terms of measurement resolution allows for ner meshes, namely, a ratio of two of the element size compared with the ZOI size yields, as a rst approximation, the same uncertainty level for the displacement and strain quantities investigated herein. This result shows that more complex displacement elds can be captured by a global approach when compared with a local approach. Furthermore, there is an ad-ditional bonus induced by the continuity constraint on the displacement eld, though less signicant (i.e., of the order of 10 %) as compared with the previous eect. It was shown that in case of very noisy images, this requirement made the displacement measurements more robust with a global approach.

The next advantage of FE-based global approaches is their direct link with numerical simulations. The measured displacement elds are very often used for comparison purposes with numerical simulations. The latter ones are, for instance, used for identication and / or validation purposes of constitutive laws and numerical models. If these issues are addressed, it is desirable to have the same kinematic basis whenever possible. If it is not the case, the measured eld can be interpolated by using the shape functions of the underlying discretization, and not some exogenous interpolation that was not used during the measurement stage. Standard strain uncertainty σ ǫ as a function of the ZOI size ℓ when σ = 1, 2, 4, 8, 16, 32 gray levels (a). The higher σ the higher σ ǫ for a given ZOI size. Normalized standard strain uncertainty σ ǫ ∆f / √ 2σ as a function of the ZOI size ℓ. The solid line shows an interpolation by the inverse of square of the ZOI size.

Figure 4 (

 4 Figure 4(a) shows the considered reference picture of the random speckle pattern.

  Figure 11(b) shows the considered reference picture. It is the raw surface with a 512 × 512 pixel denition. Even though the quality of the texture seems poor at rst sight, the histogram (Figure11(c)) indicates that most of the 12-bit dynamic range is used with a slight saturation so that image acquisition is close to optimal. The correlation radius of the texture is of the order of 3 pixels and isotropic. Hence, the images are considered as suitable for DIC.

  (a) the change of the standard uncertainty of the measured degrees of freedom for the three approaches is shown. All the general trends observed in Figure 14(b) are also present in this dicult case. In particular, the global approaches become even more interesting as they allow for a signicant decrease of the uncertainty level.

Figure 17 (

 17 Figure 17(b) shows the change of the standard uncertainty of the average displacement per ZOI (local approach) or element (global approach) as a function of ℓ. For

1 29 2 30 3 31 4 32 5 33 6 34 7 35 8Fig. 1 Fig. 2 Standard

 2930313233343512 Fig. 1 Map of the standard displacement uncertainty σu (normalized by √ 2σp/G f ℓ) for local (a) and global (b) approaches. Cuts along the diagonal x = y (c).

Fig. 3 Fig. 4 Fig. 5

 345 Fig. 3 Map of the eigenvalues in Fourier space.

Fig. 6 StandardFig. 7

 67 Fig. 6 Standard displacement uncertainty σa as a function of the ZOI size ℓ when σ = 1, 2, 4, 8, 16, 32 gray levels (a). The higher σ the higher σa for a given ZOI size. Normalized standard displacement uncertainty σa∆f / √ 2σ as a function of the ZOI size ℓ. The straight line shows an interpolation by the inverse of the ZOI size.

  Fig.8

Fig. 9 Normalized

 9 Fig.9Normalized standard strain uncertainty pℓσ ǫ /σa as a function of σa. The solid line shows the result of Equation (27).

Fig. 10 Fig. 12 Normalized

 1012 Fig. 10 Standard displacement uncertainty of the measured degrees of freedom (a), of the mean displacement per ZOI (b), of the mean strain per element for a global and optimized approach as functions of the same quantity for a local approach.

Fig. 13

 13 Fig. 13 Standard displacement uncertainty σa (a) and strain uncertainty σ ǫ (b) as functions of the ZOI size ℓ for the ten analyzed pictures. The gray shaded areas show the predictions based on the a priori resolution analysis.

Fig. 14 Fig. 15 ReferenceFig. 16 Normalized

 141516 Fig. 14 Standard displacement uncertainty of the measured degrees of freedom (a), of the mean displacement per ZOI or element (b), of the mean strain per ZOI or element as functions of ZOI or element size ℓ for the three approaches used herein.

Fig. 17

 17 Fig. 17 Standard displacement uncertainty of the measured degrees of freedom (a), of the mean displacement per ZOI or element (b), of the mean strain per ZOI or element as functions of ZOI or element size ℓ for the three approaches used herein.

  [M] is easily diagonalized. Translational invariance makes it a convolution. It is thus diagonal in Fourier space, and hence all harmonic functions e ik•x are eigenmodes, where k is the wave vector. Matrix [M]/(G 2 f ℓ 2 ) extends over second neighbors with three dierent values for inner nodes, 16/36, 4/36, and 1/36 respectively for diagonal, rst and second neighbors.

	Eigenvalues, ξ, are determined by applying [M] onto such Fourier modes. They
	obey
	ξ =	(G 2 f ℓ 2 ) 9	(4 + 2 cos(kx) + 2 cos(ky) + cos(kx) cos(ky))

  that case as well. All these results show that there is an additional bonus associated with the continuity requirement on the uncertainty level of the nodal displacements. If this eect is neglected, it is shown that the same uncertainty level is reached for a global approach as for a local one, but with a mesh that is two times ner in both directions. This can be understood by comparing the number of unknown degrees a freedom per direction eα, which is equal to 4NxNy when Nx × Ny Q4 ZOIs are analyzed in a local approach, and (Nx +

1)(Ny + 1) when Nx × Ny Q4 elements are considered in a global approach. The variance level is estimated as the ratio of the total number of pixels in the ROI Lx × Ly divided by the number of unknowns. Consequently, the gain in uncertainty level is (Nx + 1)(Ny + 1)/4NxNy ≈ 1/2 when Nx, Ny ≫ 1 when the same Q4 discretization is compared for a local and a global approach.

  2 × 2 elements. The

	mean displacement gradients become
	uα,α =	a e α3 -a e α1 + a e α6 -a e α4 + a e α9 -a e α7 6ℓp
	u α,β =	a e α7 -a e α1 + a e α8 -a e α2 + a e α9 -a e α3 6ℓp

Acknowledgments

The authors wish to thank the members of the Micromechanical Testing Laboratory at EDF R & D Les Renardières for kindly providing the two SEM pictures analyzed herein.

The denition of the picture is equal to 256 kpixels with a 12-bit digitization.


F Jendrzejewski

K Müller

J Richard

A Date

T Plisson

P Bouyer

A Aspect

V Josse

Coherent Backscattering of Ultracold Atoms

Keywords: numbers: 03.75.-b, 67.85.-d, 05.60.Gg, 42.25.Dd, 72.15.Rn

We report on the direct observation of coherent backscattering (CBS) of ultracold atoms in a quasi-two-dimensional configuration. Launching atoms with a well-defined momentum in a laser speckle disordered potential, we follow the progressive build up of the momentum scattering pattern, consisting of a ring associated with multiple elastic scattering, and the CBS peak in the backward direction. Monitoring the depletion of the initial momentum component and the formation of the angular ring profile allows us to determine microscopic transport quantities. We also study the time evolution of the CBS peak and find it in fair agreement with predictions, at long times as well as at short times. The observation of CBS can be considered a direct signature of coherence in quantum transport of particles in disordered media. It is responsible for the so called weak localization phenomenon, which is the precursor of Anderson localization.

Quantum transport differs from classical transport by the crucial role of coherence effects. In the case of transport in a disordered medium, it can lead to the complete cancelling of transport when the disorder is strong enough: this is the celebrated Anderson localization (AL) [1]. For weak disorder, a first order manifestation of coherence is the phenomenon of coherent backscattering (CBS), i.e., the enhancement of the scattering probability in the backward direction, due to a quantum interference of amplitudes associated with two opposite multiple scattering paths [2][3][4] (see inset of Fig. 1). Direct observation of such a peak is a smoking gun of the role of quantum coherence in quantum transport in disordered media. CBS has been observed with classical waves in optics [5][6][7][8], acoustics [9,10], and even seismology [11]. In condensed matter physics, CBS is the basis of the weak localization phenomenon (see e.g. [START_REF] Altshuler | Quantum Theory of Solids[END_REF]), which is responsible for the anomalous resistance of thin metallic films and its variation with an applied magnetic field [START_REF] Altshuler | [END_REF]14]. In recent years, it has been possible to directly observe Anderson localization with ultracold atoms in one dimension [15,16] and three dimensions [17,18]. Convincing as they are, none of these experiments includes a direct evidence of the role of coherence.

In this Letter, we report on the direct observation of CBS with ultracold atoms, in a quasi-two-dimensional (2D) configuration [19]. Our scheme is based on the proposal of Ref. [21] that suggested observing CBS in the momentum space. A cloud of noninteracting ultracold atoms is launched with a narrow velocity distribution in a laser speckle disordered potential (Fig. 1). Time of flight imaging, after propagation time t in the disorder, directly yields the momentum distribution, as shown in Fig. 2. As expected for elastic scattering of particles, we observe a ring that corresponds to a redistribution of the momentum directions over 2π while the momentum magnitude remains almost constant. The evolution of the initial momentum peak and of the angular ring profile yields the elastic scattering time and the transport time. But the most remarkable feature is the large visibility peak, FIG. 1. Experimental set-up. A cloud of noninteracting ultracold atoms, released from an optical trap (beams along y and z axis, red) and suspended against gravity by magnetic levitation (horizontal coils, yellow), is launched with a welldefined momentum pi along the z axis. It is submitted to an anisotropic laser speckle disordered potential (blue beam), propagating along the x axis and elongated along that direction, leading to a quasi-2D diffusive motion in the y -z plane (see text). The atomic momentum distribution in this plane is monitored by fluorescence imaging after a time of flight of 150 ms. Inset: physical origin of CBS. The coherent enhancement of scattering in the backward direction originates from the interference between each multiple scattering path (solid line) and its reversed counterpart (dashed line). which builds up in the backward direction. The height and width of that peak, and their evolution with time, are an indisputable signature of CBS, intimately linked to the role of coherence.

To understand the origin of that CBS peak, let us consider an input plane matter wave with initial momentum p i that experiences multiple scattering towards a final momentum p f (inset of Fig. 1). For each multiple scattering path, we can consider the reversed path with the same input p i and output p f . Since the initial and final atomic states are the same, we must add the two corresponding complex quantum amplitudes, whose phase difference is δφ = (p i + p f ) • R/ (R is the spatial separation between the initial and final scattering events and = h/2π the reduced Planck constant). For the exact backward momentum p f = -p i , the interference is always perfectly constructive, whatever the considered multiple scattering path. This coherent effect survives the ensemble averaging over the disorder, so that the total scattering probability is twice as large as it would be in the incoherent case. For an increasing difference between p f and -p i , the interference pattern is progressively washed out as we sum over all interference patterns associated with all multiple scattering paths. It results in a CBS peak of width inversely proportional to the spread ∆R in the separations [START_REF]More precisely[END_REF]. For diffusive scattering paths, the distribution of R is a Gaussian whose widths increase with time as t 1/2 , and the CBS widths decrease according to ∆p CBS,α = / √ 2D α t for each direction of space (α = y, z), D α being the diffusion constant along that direction. This time resolved dynamics of the CBS peak has been observed in acoustics [9,10] and optics [START_REF] Vreeker | [END_REF]24].

The crux of the experiment is a sample of noninteracting paramagnetic atoms, suspended against gravity by a magnetic gradient (as in [18]), and launched along the z axis with a very well-defined initial momentum p i (see Fig. 1). This is realized in four steps. First, evaporative cooling of an atomic cloud of 87 Rb atoms in a quasiisotropic optical dipole trap (trapping frequency 5 Hz) yields a Bose-Einstein condensate of 9 × 10 4 atoms in the F = 2, m F = -2 ground sublevel. Second, we suppress the interatomic interactions by releasing the atomic cloud and letting it expand during 50 ms. At this stage, the atomic cloud has a size (standard half-width along each direction) of ∆r α = 30 µm, and the residual interaction energy (E int /h ∼ 1 Hz) is negligible compared to all relevant energies of the problem. Since the atomic cloud is expanding radially with velocities proportional to the distance from the origin, we can use the "deltakick cooling" technique [26], by switching on a harmonic potential for a well chosen amount of time. This almost freezes the motion of the atoms, and the resulting velocity spread ∆v α = 0.12±0.03 mm/s is just one magnitude above the Heisenberg limit (∆r α m∆v α ∼ 5 , with m the atom mass). Last, we give the atoms a finite momentum p i along the z direction, without changing the momentum spread, by applying an additional magnetic gradient during 12 ms. The first image of Fig. 2 shows the resulting 2D momentum distribution. The average velocity is

v i = 3.3 ± 0.2 mm/s (k i = p i / 4.5 µm -1
), corresponding to a kinetic energy E K = p 2 i /2m (E K /h 1190 Hz). This momentum distribution is obtained with a standard time of flight technique that converts the velocity distribution into a position distribution. Because of the magnetic levitation, we can let the atomic cloud expand ballistically for as long as 150 ms before performing fluorescence imaging along the x axis. The overall velocity resolution of our experiment that takes into account the initial momentum spread, which writes ∆v res = [∆v α 2 + (∆r α /t tof ) 2 ] 1/2 = 0.23 mm/s, is nevertheless mainly limited by the size ∆r α of the atomic cloud.

To study CBS, we suddenly switch on an optical disordered potential in less than 0.1 ms, let the atoms scatter for a time t, then switch off the disorder and monitor the momentum distribution at time t. The disordered potential is the dipole potential associated with a laser speckle field [27,[START_REF] Goodman | Speckle Phenomena in Optics: Theory and Applications[END_REF], obtained by passing a laser beam through a rough plate, and focusing it on the atoms (Fig. 1). It has an average value V R (the disorder "amplitude") equal to its standard deviation. Its autocorrelation function is anisotropic, with a transverse shape well represented by a Gaussian of standard half-widths σ y = σ z = σ ⊥ 0.2 µm, and a longitudinal Lorentzian profile of half-width σ x 1 µm (HWHM) [START_REF] Piraud | [END_REF]. The laser (wavelength 532 nm) is detuned far off-resonance (wavelength 780 nm), yielding a purely conservative and repulsive potential. The disorder amplitude V R is homogenous to better than 1% over the atom cloud (profile of half-widths 1.2 mm along y,z, 1 mm along x).

The anisotropy of the speckle autocorrelation function (elongated along x) allows us to operate in a quasi-two-dimensional configuration by launching the atoms perpendicularly to the x axis (along the z axis). In the y-z plane, the atoms are scattered by a potential with a correlation length shorter than the reduced atomic de Broglie wavelength (k i σ ⊥ 0.9), so that the scattering probability is quasi-isotropic, and we will replace the subscript α = y, z by ⊥ in the rest of this Letter. The dynamics within this plane develops on the typical time scale of a single scattering event, that is, the elastic scattering time τ s . In contrast, the correlation length along the x axis is larger than the reduced atomic de Broglie wavelength (k i σ x 4.5), so that each scattering event produces but a small deviation out of the y-z plane [30]. The diffusive dynamics along x is then slower than in the yz plane, and for short times the dynamics is quasi 2D.

Figure 2 shows the time evolution of the momentum distribution for a disorder amplitude V R /h = 975±80 Hz. In order to analyze these data quantitatively, we perform a radial integration of the 2D momentum distribution on a thin stripe between p i -∆p res and p i + ∆p res (inset of Fig. 3) (∆p res = m∆v res is the momentum resolution). This yields the angular profile n(θ, t), displayed in Fig. 3 for increasing diffusion times.

We first extract the elastic scattering time τ s from the exponential decay of the initial peak n(θ = 0, t) ∝ e -t/τs [START_REF] Akkermans | Montambaux Mesoscopic Physics of Electrons and Photons[END_REF]. We find τ s = 0.33 ± 0.02 ms (mean free path l s = v i τ s = 1.1 µm). Note that this quantity also plays a role in the radial width of the ring (inset of Fig. 3), which is associated to a Lorentzian energy spread ∆E dis = /τ s (HWHM) acquired by the atoms when the disorder is suddenly switched on. Combining the corresponding momentum spread with the resolution of our measurement, we find a width in agreement with the observed ring width. The measured value of τ s is in accordance with numerical simulations adapted to our configuration, but is about 2 times the value predicted by a perturbative calculation [START_REF] Piraud | [END_REF]. This is consistent with the fact that we are not fully in the weak disorder regime defined by ∆E dis /E K = 2/k i l s

1. Here we have ∆E dis /E K ∼ 0.4 (k i l s ∼ 5) [32].

Monitoring the isotropization of the momentum distribution, we obtain another important quantity: the transport time τ , after which, in the absence of coherence, the information about the initial direction would be lost. It is determined from the exponential decay e -t/τ of the first component of the Fourier series expansion of n(θ, t) [21,34]. We find τ = 0.4 ± 0.05 ms, also in agreement with numerics. Note that τ takes into account the CBS phenomenon and its calculation must include weak localization corrections. It is only slightly larger than τ s , as expected for a nearly isotropic scattering probability (k i σ ⊥ = 0.9). The transport time sets the time scale of the onset of the diffusive dynamics, which is well established only after several τ . In the experiment, we observe that the momentum distribution has become fully isotropic after t ∼ 2.5 ms (i.e. ∼ 6 τ ), with a steady and flat angular profile of mean value n(θ, t) = n, except around θ = π where the CBS peak is still present.

To analyze the evolution of the CBS signal, we fit (see Fig. 3) the angular profiles by the function n incoh (θ, t) + n coh (t) exp[-(θπ) 2 /2∆θ(t) 2 ]. In this formula, n incoh is (within higher order terms) the multiple scattering background that would be obtained in the absence of coherence [START_REF] Akkermans | Montambaux Mesoscopic Physics of Electrons and Photons[END_REF] and we assume that n incoh (θ, t) has a parabola shape around θ = π. The fit allows us to measure the contrast C(t) = n coh (t)/n incoh (π, t) and the width ∆θ(t) of the CBS peak, and to compare them to the results of theoretical predictions. Their evolutions are shown in Fig. 4. A CBS peak appears as soon as scattering in the backward direction is significant, but the contrast starts decreasing before reaching the ideal value of 1. For the width, we observe the predicted monotonic decrease, but it tends asymptotically towards a non-null value rather than zero.

To compare these observations with theoretical predictions, we must take into account our finite resolution. The black line in Fig. 4(b) represents the calculated CBS width that results from the convolution of our resolution ∆θ res = ∆p res /p i = 0.07 with the expected CBS width ∆θ CBS = ∆p CBS,⊥ /p i = /p i √ 2D ⊥ t in the diffusive regime (the widths are added quadratically). The diffusion constant is evaluated using the standard relation D ⊥ = v 2 i τ /2, so that the solid line does not involve any adjustable parameter. We see that the agreement with the data is good when we enter the multiple scattering regime [for t 1.5 ms in Fig. 4(b)], but not at The theoretical predictions for the multiple scattering regime are represented by the solid black lines at long times (t 4τ ; τ = 0.4 ms), and by dotted black lines at short times. The dashed-dotted red lines correspond to the calculation of [START_REF] Gorodnichev | [END_REF] at short times (i.e. for t ∼ τs; τs = 0.33 ms), where single scattering events cannot be neglected.

short times.

The broadening of the CBS peak by the finite resolution is also responsible for a decrease of the contrast, as represented in Fig. 4(a). Here also, we observe that the prediction (which again involves no adjustable parameter) is very different from the observed values at short times, but is in reasonable agreement with the measurements around t ∼ 1.5 ms. On the other hand, the measured contrast is definitely smaller than the theoretical prediction when t increases yet more. We relate this observation at long times to the onset of the dynamics in the third direction x. Using a second imaging system yielding the momentum distribution in the xy plane, we estimate a typical time of 4 ms for this out-of-plane dynamics to become significant, and render the 2D approach wanting.

Deviations at short times were to be expected. Indeed, CBS demands multiple scattering, or at least double scattering, to happen (see inset of Fig. 1), whereas single scattering events do not participate to the CBS peak. At short times (t ∼ τ s ), the contribution of single scattering to backscattering is not negligible compared to multiple scattering. This entails a reduction of the contrast (see e.g. [3]), and a modification of the shape (no longer Gaussian), whose width decreases at this stage as 1/t (ballistic motion between the first two scatterers). In the case of light, a calculation for isotropic scattering [START_REF] Gorodnichev | [END_REF] predicts a short time evolution of the contrast C = (2t/πτ s )/(1 + 2t/πτ s ) and width ∆θ CBS ∼ 3/k i l s (τ s /t). This prediction is plotted in Fig. 4 and is found in fair agreement with the observations in this time domain. Finally, note that the width around t ∼ τ s is linked to the disorder strength quantified by k i l s . Here we find a maximum value of ∆θ max ∼ 0.3 rad, that is ∆θ max ∼ 1.5/k i l s (k i l s ∼ 5, see above).

Similar measurements and analysis have been repeated for weaker disorder (V R /h = 525, 750 Hz) and smaller initial momentum p i (E K /h = 160, 220, 620 Hz), and we have found a similar agreement between data and theory. In contrast to the observed moderate changes in the maximum peak contrast (C max ∼ 0.5 -0.7), the maximum peak width ∆θ max increases significantly with the amplitude of the disorder and the inverse of p i . The highest observed width of 1.2 rad (from which we infer k i l s ∼ 1.25) suggests that we are very close to the strong disorder regime, where AL is expected to be experimentally observable in 2D systems. Such an observation, however, would demand a longer 2D evolution in the disorder, which is limited in the present experiment because of the cross over to the 3D regime. Increasing this time, as for instance in [37], will then constitute the next step towards AL, with the possibility to observe the coherent forward scattering peak predicted in [38].

In conclusion, we have demonstrated experimentally that the time resolved study of the momentum distribution of ultracold atoms in a random potential is a powerful tool to study quantum transport properties in disordered media. We have been able to extract the elastic scattering time, the transport time, and to observe and study the evolution of the CBS peak. Let us emphasize that the theoretical analysis as well as numerical simulations render an account of the observations not only in the multiple scattering regime but also at short time, during the onset of multiple scattering. Such agreement gives a strong evidence of the fundamental role of coherence in that phenomenon. Further evidences of the role of coherence could be sought in the predicted suppression of the CBS peak [39] when scrambling the disorder, or when dephasing the counter-propagating multiple scattering paths using artificial gauge fields [40], in the spirit of pioneering works in condensed matter physics [14] or optics [41]. Finally, this work also opens promising prospects to study the effect of interactions on CBS (see e.g. [42,43]).

FIG. 2 .

 2 FIG.2. Observed momentum distributions after different propagation times t in the disorder. The images correspond to an averaging over 20 experimental runs. Note that the vertical scale is different in the three first (t =0, 0.5, and 1 ms), whereas it is the same in the three last images (t =1.5, 2, and 2.5 ms).

FIG. 3 .

 3 FIG. 3. Normalized angular profiles n(θ, t)/n corresponding to the distributions shown in Fig. 2. The red solid lines correspond to the double structure fit around the backscattering direction [a parabola for the incoherent background (dashed line) and a Gaussian for the CBS peak, see text]. Inset: false color representation of the momentum distribution (t = 2 ms). The angle θ refers to the scattering direction with respect to the initial direction (pi).

FIG. 4 .

 4 FIG.4. Time resolved dynamics of the CBS peak: evolution of the contrast C (a) and the width ∆θ (b) versus the scattering time t. The blue points are experimental data and the error bars correspond to the 95% confidence intervals for the fitted values (±2 standard deviations). The theoretical predictions for the multiple scattering regime are represented by the solid black lines at long times (t 4τ ; τ = 0.4 ms), and by dotted black lines at short times. The dashed-dotted red lines correspond to the calculation of[START_REF] Gorodnichev | [END_REF] at short times (i.e. for t ∼ τs; τs = 0.33 ms), where single scattering events cannot be neglected.

We thank T. Bourdel, C. Müller and B. van Tiggelen for fruitful discussions and comments. This research was supported by ERC (Advanced Grant "Quantatop"), ANR (ANR-08-blan-0016-01), MESR, DGA, RTRA Triangle de la Physique, IXBLUE and IFRAF.


Caroline Lafond
email: c.lafond@rennes.unicancer.fr

Frédéric Gassa

Christophe Odin

Gaël Dréan

Justine Even

Renaud De Crevoisier

Pascal Pommier

Jean-Pierre Manens

Marie-Claude Biston

Comparison between two treatment planning systems

Keywords: VMAT, Pinnacle, Monaco, Prostate Cancer

teaching and research institutions in France or abroad, or from public or private research centers.

INTRODUCTION

Volumetric modulated arc therapy (VMAT) is a new radiotherapy technique which allows to achieve treatment plans of similar or improved quality compared to fixed-field intensitymodulated radiation therapy (IMRT) while reducing the treatment time per fraction [START_REF] Otto | Volumetric modulated arc therapy: IMRT in a single gantry arc[END_REF]. In practice, to obtain highly modulated dose distributions delivered efficiently, a treatment planning system (TPS) with a powerful optimization and segmentation algorithm is required.

While a lot of users are in the process of replacing fixed-field IMRT by VMAT, or directly implementing VMAT in their radiotherapy department, there is a lack of information concerning the relative performances of the mainly used TPS for VMAT planning. To our knowledge, only three studies deal with this topic [START_REF] Rao | Comparison of Elekta VMAT with helical tomotherapy and fixed field IMRT: plan quality, delivery efficiency and accuracy[END_REF][START_REF] Masi | Comparison of 3 TPS for VMAT Plans Optimization[END_REF][START_REF] Wiezorek | Rotational IMRT techniques compared to fixed gantry IMRT and tomotherapy: multi-institutional planning study for head-and-neck cases[END_REF]. In Rao et al, ERGO++ (Elekta, Crawley, UK) was compared to Pinnacle (Philips Medical Systems, Madison, WI) direct machine parameter optimization (DMPO) combined with a home-made arc-sequencer and Pinnacle SmartArc inverse planning module [START_REF] Rao | Comparison of Elekta VMAT with helical tomotherapy and fixed field IMRT: plan quality, delivery efficiency and accuracy[END_REF]. In Masi et al, the performances of Monaco (CMS-Elekta, Crawley, UK) were compared to ERGO++ and Oncentra (Nucletron-Elekta) [START_REF] Masi | Comparison of 3 TPS for VMAT Plans Optimization[END_REF]. Finally, in Wiezorek et al, VMAT plans obtained with Monaco and Eclipse (Varian Medical System, Palo Alto, CA) were evaluated [START_REF] Wiezorek | Rotational IMRT techniques compared to fixed gantry IMRT and tomotherapy: multi-institutional planning study for head-and-neck cases[END_REF]. In these studies, the comparisons were made by fixing common planning objectives on PTVs and OARs and comparing the dosimetric results and treatment delivery efficiency (number of monitor units and treatment time).

The aim of the present study was to investigate the performances of two TPS that have not been compared yet in VMAT mode, both using different approaches for VMAT plan optimization: Monaco based on a two-stage constrained optimization [START_REF] Semenenko | Evaluation of a commercial biologically based IMRT treatment planning system[END_REF] and Pinnacle SmartArc [START_REF] Bzdusek | Development and evaluation of an efficient approach to volumetric arc therapy planning[END_REF]. This work was performed by two institutions. The aim was to compare VMAT plans performed by Monaco and Pinnacle regarding to dosimetric performances and treatment delivery specificities. We therefore fully put in evidence the differences observed in terms of dose distributions, delivery efficiency, treatment delivery parameters (mean dose rate, mean segment area) and quality control results on 10 prostate cancer cases.

MATERIALS AND METHODS

Patients

Ten prostate adenocarcinoma patients referred to our institutions for a radical external beam irradiation to the prostate and seminal vesicles (SV) were considered for this dosimetric comparative analysis.

Anatomic data acquisition, volumes definition and dose

Organs at risk [rectal wall (5 mm thickness), bladder wall (7 mm thickness), femoral heads (FH)] and target volumes (prostate, SV) were delineated on dedicated 2 mm-thick CT slices.

The first clinical target volume (CTV1) comprised the prostate and SV. The CTV2 was limited to the prostate only. Planning target volumes (PTVs) were automatically generated adding a 3D 1 cm uniform margin around the CTVs, except in the posterior direction, where a 0.5 cm margin was added to protect the rectum.

The total prescription dose was 46 Gy to the PTV1 and an additional 32 Gy to the PTV2 using a standard fractionation (2 Gy per fraction, 5 days a week) for a total dose of 78 Gy using a sequential technique.

A dose objectives set was fixed for PTVs and OARs : for PTV1 : 95% of the PTV covered by 97% of the prescribed dose, and less than 5 % of the PTV receiving more than 107% of the prescribed dose; PTV2: 95% of the PTV covered by 95% of the prescribed dose; Rectum: maximum dose (into 1.8 cc) < 76 Gy, V72 1 25% , V60 1 50% ; Bladder V70 1 25%, V60 1 50%; Femoral heads: V50 1 5%.

Treatment planning

For each patient, two VMAT plans were generated: one with Monaco 3.0 (CMS-Elekta Ltd, Crawley, UK) and one with Pinnacle 9.0 (Philips Medical Systems, Madison, WI).

The irradiation was delivered, using 6-MV photons with an Elekta Synergy ® machine equipped with a Cone-Beam Computed Tomography (CBCT) device (XVI ® ) and with a multi-leaf collimator (MLCi2 ® ) consisting of 40 paired leaves, each measuring 1 cm in width at the isocenter. The possible dose rate values were 25 MUs/min, 50 MUs/min, 100 MUs/min, 200 MUs/min and 400 MUs/min. For each treatment plan a single 360° arc was used.

MONACO PLANNING

For Monaco planning, the optimization constraints were established on the basis of biological cost functions (i.e. Serial or parallel complication model for OARs and Poisson cell kill function for the PTVs). The prescription template applied to all patients is given in Table 1. The optimisation was first performed in a constrained mode, meaning that all constraints to the OARs are treated as hard constraints and all optimization criteria must be met. Conversely, the constraints to the targets are considered as objectives. The pareto mode which gives priority to PTV coverage was used secondarily to achieve the PTV coverage detailed above.

Sequencing parameters used for PTV1 and PTV2 irradiation were: 124 control points (CP) to achieve in practice 120 CP; target dose rate 300 MUs/min; minimum segment width 0.5 cm; fluence smoothing: low.

For final Monte Carlo dose calculations, a calculation grid of 3 mm and a 3 % variance were used. With these parameters, the time needed for final dose calculation was about 10 min on an Intel Xeon CPU 3 GHz and 12 GB RAM platform. The time for optimisation stage and adjusting the prescription parameters was about 20 min.

PINNACLE PLANNING

For Pinnacle planning, inverse optimization was performed using the SmartArc algorithm [START_REF] Bzdusek | Development and evaluation of an efficient approach to volumetric arc therapy planning[END_REF].

The optimization objectives were defined with physical dose points. The template is shown in Table 1. The arc sampling parameter was fixed at 3 degrees to obtain 120 CP for the full arc. The delivery time parameter was fixed at 180 s firstly; then was eventually increased to 240 s to allow more dose modulation for the most complex cases. Final dose was computed with a collapsed cone algorithm using a dose grid resolution of 3 mm. With these parameters, the time needed for optimization and final dose calculation was about 13 min on an Intel quadruple-Core (Xeon) 2.8

GHz and 16 GB RAM platform. Time for parameters adjustment was 10 min.

Preliminary work

Although this study was performed by two institutions, an important number of constraints were set to limit the influence of the planners and planning philosophy of the two hospitals. First, a preliminary comparison study was performed on a water-equivalent cylindrical phantom with a C-Shape target surrounding a central avoidance structure (data not shown) as described by the AAPM task group 119 [START_REF] Ezzell | IMRT commissioning: multiple institution planning and dosimetry comparisons, a report from AAPM Task Group 119[END_REF]. This preliminary work allowed to harmonize both planning methods and to verify that for a simple geometry both institutions were able to produce plans of similar quality regarding dose distribution and delivery efficiency.

Treatment plans comparisons

DOSE DISTRIBUTION

In order to limit the uncertainties on DVHs calculations between both TPS, the results were evaluated in the ARTiView 1.12 software (Aquilab, Lille, France) by comparing DVHs for targets and OARs (mean dose and doses at selected points of the DHVs). Patient-averaged DVHs were compared. In addition, several quality indexes for PTV1 and total plans were assessed: homogeneity index (HI) was calculated as (D 5% -D 95% )/D mean within the PTV; D 5% and D 95% being the dose received by 5 and 95 % of the PTV [START_REF] Claus | Report of a study on IMRT planning strategies for ethmoid sinus cancer[END_REF] ; conformity index (CI) was calculated as the ratio between the volume of the reference isodose (V 95% ) and the PTV volume (V PTV ) [V 95% /V PTV ] [START_REF] Shaw | Radiosurgery for the treatment of previously irradiated recurrent primary brain tumors and brain metastases: initial report of radiation therapy oncology group protocol (90-05)[END_REF]; healthy tissue coverage index (HCO) evaluates the percentage of reference isodose which is outside the PTV volume. HCO was calculated as [100*(1-(V PTV, 95% /V 95% ))]; V PTV, 95% was the volume of PTV covered by the reference isodose.

Statistical analysis used two-sided Wilcoxon-signed rank test, a nonparametric test, calculated with PASW Version 18.0.0 (SPSS Inc., Chicago, IL). A value of p<0.05 was considered statistically significant.

To underscore the spatial localization differences between the two TPS, a patient-averaged dose distribution was performed. To obtain the averaged dose distribution, the anatomy of all the patients was registered to a typical patient anatomy. After an initialization done by an affine registration, we used an organ-driven non-rigid registration method using the demons algorithm between the CTs and between each considered organ (prostate, bladder, rectum). The final deformation field was computed by merging the different deformation fields by weighting them according to the distance between voxels and organs [START_REF] Thirion | Image matching as a diffusion process: an analogy with Maxwell's demons[END_REF][START_REF] Drean | Inter-individual organ-driven CT registration for dose mapping in prostate cancer radiotherapy[END_REF][START_REF] Ourselin | Reconstructing a 3D structure from series histological section[END_REF]. The transformation was then applied to the dose distribution.

ROBUSTNESS

In case of daily CBCT, we assume to have a geometric accuracy better than 2 mm as mentioned in the AAPM recommendations [START_REF] Bissonnette | Quality assurance for image-guided radiation therapy utilizing CT-based technologies: a report of the AAPM TG-179[END_REF]. To investigate the robustness of Monaco and Pinnacle dose distributions, an isocenter shift of 2 mm was applied on one representative prostate case in unfavorable directions for the main OARs (i.e. in anterior and posterior direction). The impact on dose distribution was judged with the quality indexes described previously for PTV2, rectum wall and bladder wall.

DELIVERY EFFICIENCY

The delivery time, MU per fraction, mean dose rate and time-average MLC aperture area were used to evaluate the VMAT delivery efficiency.

Quality assurance

In order to evaluate the dose algorithm accuracy of Monaco and Pinnacle, the consistency between calculated and delivered dose was verified for one typical plan with EBT3 radiochromic films (Ashland, Covington, USA) in a pelvis anthropomorphic phantom (CIRS, Norfolk, USA).

A dosimetric validation was then performed for all plans with the Octavius phantom (PTW, Freiburg, Germany). The dose delivery was measured using the PTW 2D-array Seven29 (PTW) ionization chamber matrix. The dose was measured in coronal and sagittal planes.

Comparisons were performed with Verisoft software (PTW). The 3D gamma method was used to compare the measured dose distributions with the calculated 3D dose distribution. The dose criterion was 3% of the local dose and the distance criterion was 3 mm. The evaluated areas were areas with doses higher than 30% of the maximum dose.

RESULTS

Dose distribution

The patient-averaged differences in the dose distributions for the two TPS are shown in Figure 1.

Large differences are observed. First, Monaco favoured the dose delivery on gantry angles that allowed crossing a minimum volume of healthy tissue to reach the PTV. Therefore, the volume of healthy tissue receiving higher doses was more important with Pinnacle than with Monaco. Furthermore, we observed that Pinnacle solutions result in more dose delivered on the left-right direction (i.e. around the FH), whereas Monaco solutions result in more doses delivered on the antero-posterior direction (i.e. around the bladder and the rectum).

The averaged dose-volume histograms (DVHs) of total plans are shown in Figure 2. The continuous lines represent the mean values and are surrounded by two dashed curves representing the 2.5 th percentiles of the data for the lower dashed curves and the 97.5th percentiles for the upper dashed curves. DVHs data for PTVs and OARs and conformal indexes are reported in Table 2 as averages for the investigated patients. Figure 2 shows that Pinnacle plans resulted in more dose into the PTVs than Monaco but with a more important scattering. While the mean dose to PTV1 was significantly different between both modalities (p=0.008), the difference was not statistically significant on PTV2 (p=0.241) (Table 2). Furthermore, Monaco plans had higher conformality with significantly better CI and HCO at the expense of significantly lower dose homogeneity into the PTV than Pinnacle. Regarding the dose distribution to the OARs, Monaco provided a significantly lower mean dose to the bladder wall (p=0.019), to the FH (p=0.017) and to the healthy tissue (p=0.005) while Pinnacle provided a significantly lower mean dose (p=0.047) and a lower D 50% value (p=0.047) for the rectal wall (Table 2 and Figure 2).

Robustness

Details of dose with and without isocentre shift for the main organs are shown Table 3. For the rectal wall, the 2-mm posterior shift resulted in an increase of the D max of 1.3 Gy for both TPS, and an increase of the D mean of 1.9 and 1.7 Gy for Monaco and Pinnacle plans, respectively. The D 95% to PTV2 was the same for Monaco plans and resulted in an increase of 0.3 Gy for Pinnacle plans. For the bladder wall, the 2-mm anterior shift resulted in an increase of the D max of 0.3 Gy and 0.2 Gy and an increase of the D mean of 0.8 and 0.7 Gy, for Monaco and Pinnacle plans, respectively. The D 95% to PTV2 was decreased of 0.7 and 1.3 Gy for Monaco and Pinnacle plans, respectively.

Treatment delivery evaluation

Details of efficiency parameters are shown in Table 4; Figure 3 shows the variations of MLC aperture area and dose rate as a function of gantry angle during typical VMAT delivery. Monaco plans needed one and a half more MUs than Pinnacle plans. Nevertheless we noticed similar delivery times because the time-averaged dose rates were higher with Monaco (230 MUs/min) than with Pinnacle (160 MUs/min).

The better efficiency of Pinnacle plans was due to a larger MLC aperture area (Fig. 3, Table 4).

On average, total MUs to MLC aperture area ratio was 265 MUs/cm² for Monaco plans and 110 MUs/cm² for Pinnacle plans. The variation of the area according to the gantry angle (Fig. 3) shows that the segmentation is based on an alternative Sliding Window pattern for Monaco VMAT (i.e. all the leaves moves first on one way and then on the other way alternatively for the full arc) .

Quality assurance

The analysis of the films placed inside an anthropomorphic phantom showed that 99.1% and 98.4% of the points passed the 3%/3 mm criterion with mean gamma values of 0.34 and 0.31 for Monaco and Pinnacle respectively. Central area encompassing the PTV showed high consistency between calculated and delivered dose for Pinnacle and Monaco; the mean gamma values were 0.29 for both TPS.

The results obtained with the 2D matrix ionisation chamber showed a high passing rate (> 97.4%) for the gamma index for both TPS (Table 4). Nevertheless, Monaco provided better dosimetric agreement than Pinnacle. To explain the QA results differences, we investigated the dosegradient values [START_REF] Moran | A dose gradient analysis tool for IMRT QA[END_REF] in the two measurement plans of the Octavius phantom (Table 4). For coronal planes, mean dose-gradients were similar for both TPS but maximum dose-gradient was higher with Pinnacle. For sagittal planes, all reported dose-gradient values were higher with Pinnacle.

DISCUSSION

The aim of the present study was to investigate the performances of Monaco and Pinnacle TPS for VMAT plan optimization. To our knowledge, this is the first study yielding a global comparison of two TPS for VMAT planning, from the prescription phase and dose distribution evaluation to the delivery efficiency. We also included beam geometry, treatment delivery parameters and quality control results.

For prostate cancer, VMAT solutions proposed by both TPS offered good PTV coverage and OARs sparing, with similar delivery time. Note that the lowest doses to the OARs were achieved with Monaco, except for the rectal wall. Regarding dose distribution and delivery parameters large differences were observed. First on dose distribution, we noticed that Monaco favoured the dose delivery from gantry angles that allowed crossing a minimum volume of healthy tissue to reach the PTV. As a consequence, much lower doses were delivered to the FH in comparison with Pinnacle, despite the fact that contrarily to Pinnacle, no constraint was assigned to the FH on Monaco prescription. Likewise, healthy tissue contained in the CT area received less dose with Monaco than with Pinnacle. However, this result must be taken with care since the uncertainties in TPS dose calculation are generally larger in the lower dose regions. Indeed, the uncertainty in Monte Carlo calculation is influenced by the variance reduction techniques and inversely proportional to the square root of the number of histories used for dose calculation. The uncertainty of the entire plan is always less than the variance value for dose calculation since that value is per segment (i.e. 3 % in our case, for Monaco TPS). Therefore the noise which can be observed is substantially reduced in high dose area where a more important number of segments contribute to the dose (superposition effect) [START_REF] Fotina | Advanced kernel methods vs. Monte Carlo-based dose calculation for high energy photon beams[END_REF] . Furthermore, it is well known that the main part of out-of-field doses are due to the linac-head scatter and leakage radiation, proportional to MUs [16;17]. In this case Monaco could provide a higher out-of-field dose. Hence, both TPS provide different irradiation patterns regarding low doses without superiority of one TPS was clearly demonstrated.

Finally, higher doses were delivered in the antero-posterior direction, above the bladder or below the rectum with Monaco compared to Pinnacle. Therefore, with Monaco TPS, an anatomical deformation of the rectum or bladder will probably have more consequences on the delivered dose to these organs. In particular, it has been shown that the relative volume of the rectum can vary by more than 150% from one day to another [START_REF] Cazoulat | From imageguided radiotherapy to dose-guided radiotherapy[END_REF]. Thus in this case, the risk of rectal toxicity is increased and it is of importance to make a daily imaging control. The robustness investigation of dose distributions showed that the impact of 2 mm isocentre shifts provided very low and similar dose deviations for Monaco and Pinnacle plans : Table 3 showed that the dose deviations were always less than 5% of the prescribed dose and the dose deviations differences were less than 1% of the prescribed dose between the both TPS. Therefore both TPS provided robust plans.

To conclude on this part, most of the differences observed may be closely linked to the objective functions used for optimizing the dose distribution. Whereas on both TPS biological and physical cost functions are available, Monaco planner chose to use biological in combination with physical cost functions while Pinnacle planner used only physical cost functions. In one study IMRT plans performed with Pinnacle and Monaco biological-dose based prescriptions were compared to Pinnacle physical-dose based prescriptions [START_REF] Qi | Improved critical structure sparing with biologically based IMRT optimization[END_REF]. Biological prescription led to improve OARs sparing compared to physical dose based prescription with a similar performance for Monaco and Pinnacle. Contrarily to Qi et al, another study comparing VMAT plans performed with Pinnacle biological, physical and mixed prescriptions reported that, in some particular cases, the physical prescriptions result in superior treatment plans [START_REF] Mihaylov | Biological optimization in volumetric modulated arc radiotherapy for prostate carcinoma[END_REF]. Therefore, there is still an open debate on the question of the superiority of the biological against physical cost functions.

In practice, it is possible to get the same results using physical or biological cost functions even though physical cost functions are not as easy to control since one constraint only controls one dose point. Therefore one single biological cost function will always lead to a better dosimetric solution than one single physical cost function, but is this true for one biological cost function against 3 physical cost functions? There will always be a high proportion of subjectivity concerning all the studies performed on this subject since the results highly depend:

1/ on the planner and his degree of experience with either physical or biological cost functions 2/ on how many physical and biological cost functions are used in the prescription The above papers did not give any details about the prescription templates the authors used to get their results. In the present study, the planners had to use the prescription template that they had found to be best suitable for treating prostate cancer. Monaco planner chose to combine biological and physical cost functions, since from their practices, the best results were obtained with mixed prescriptions, whereas Pinnacle planner preferred physical cost functions, since they did not get better results with biological prescriptions.

Regarding the treatment delivery efficiency, the treatment control system of the Elekta accelerator adjusts the dose rate, the gantry and leaves speeds to offer the lowest delivery time possible for VMAT plans. This adjustment is efficient since despite different numbers of MUs, Monaco and Pinnacle plans were delivered with similar times. We noticed this adjustment particularly with the dose rate variation; Monaco plans were delivered with time-average dose rate 44% higher than Pinnacle plans. Previous studies compared VMAT performed with Monaco or other TPS [3;4]. Similarly to our results they found that Monaco provided more MUs than other TPS.

CONCLUSION

For prostate cancer patients, VMAT planned with Monaco and Pinnacle TPS offered clinically acceptable dose distributions. Monaco plans showed enhanced OAR sparing but lower doses into the PTV compared to Pinnacle plans. Similar delivery times were found for both TPS but Pinnacle solution required less MUs. Finally, a good dosimetric agreement with measured doses was achieved with both TPS, but Monaco offered a slightly higher passing rate in the gamma index analysis. Further investigations are in progress to confirm the performances of both TPS on more complex volumes (head and neck cancer or prostate with pelvic node involvement).

Conflicts of interest:

Caroline Lafond' contribution is part of a PhD thesis supported by a grant from Elekta.

Figures and Tables legends

  Figures and Tables legends:

Figure 1 .

 1 Figure 1. Patient-averaged dose distribution differences for axial (A), coronal (B) and sagittal (C)

Figure 2 .

 2 Figure 2. Average composite DVHs for VMAT performed either with Monaco (continuous grey

Figure 3 .

 3 Figure 3. Area and dose rate variation as a function of gantry angle during VMAT delivery for a

Figure 3

 3 Figure 3

Table 2 .

 2 Average dosimetric indices of bladder wall, rectum wall, femoral heads, healthy tissue, PTV1 and PTV2 from VMAT with Monaco or Pinnacle. P<0.05 was considered significant for

	Wilcoxon test.

Table 2

 2

		Monaco	Pinnacle	Wilcoxon test		Monaco	Pinnacle	Wilcoxon test
	Rectal wall				PTV1
	D max (Gy)	75.1 (0.5)	74.7 (0.8)	0.385	D mean (Gy)	73.2 (2.8)	74.9 (3)	0.008
	D mean (Gy)	46.1 (5.7)	43.4 (7.2)	0.047
	D 25% (Gy)	61.7 (7.6)	60.5 (8.5)	0.093
	D 50% (Gy)	45.4 (8.4)	41.5 (9.9)	0.047
	Bladder wall				PTV2
	D max (Gy)	77.9 (0.7)	78.1 (1.1)	0.541	D mean (Gy)	78.4 (0.4)	78.9 (0.8)	0.241
	D mean (Gy)	31.9 (7.3)	33.9 (8.3)	0.019	D max (Gy)	81.7 (0.7)	81.1 (1.1)	0.102
	D 25% (Gy)	52.1 (12.6)	54.5 (11.5)	0.221	D 5% (Gy)	81.1 (0.6)	80.8 (1)	0.260
	D 50% (Gy)	22.4 (9.4)	24.7 (10.9)	0.202	D 95% (Gy)	74.7 (0.2)	75.4 (0.6)	0.005
	Femoral heads
	D mean (Gy)	19.6 (3.5)	25.3 (2.7)	0.017
	D 5% (Gy)	33.1 (5.3)	37.1 (2.5)	0.059
	Patient				Quality indexes
	D mean (Gy)	8.3 (1.2)	9.2 (1.4)	0.005	HI	0.08 (0.01)	0.07 (0.01)	0.027
	D max (Gy)	82.1 (0.7)	81.2 (1.1)	0.047	CI	1.19 (0.08)	1.30 (0.08)	0.036
	V 5% (cc)	34.8 (4.3)	36.4 (4.7)	0.008	HCO	18.10 (4.40)	25.37 (6.04)	0.028
	V 10% (cc)	29.4 (3.7)	30.5 (4.1)	0.053
	V 95% (cc)	1.1 (0.3)	1.2 (0.4)	0.012

Table 3

 3

			Without shift	2 mm shift in anterior direction	2 mm shift in posterior direction
		Monaco	Pinnacle	Monaco	Pinnacle	Monaco	Pinnacle
	PTV2
	D mean (Gy)	78.8	78.5	78.6	78.3	78.8	78.7
	D max (Gy)	83.0	81.1	82.8	81.0	82.9	81.2
	D 5% (Gy)	81.9	80.6	81.8	80.6	81.9	80.8
	D 95% (Gy)	74.8	74.9	74.1	73.6	74.8	75.2
	Rectum wall
	D max (Gy)	74.0	74.4	72.7	73.2	75.3	75.7
	D mean (Gy)	46.4	47.4	44.7	45.8	48.3	49.1
	D 25% (Gy)	65.3	65.3	61.9	61.5	68.3	68.5
	D 50% (Gy)	46.9	47.9	45.3	46.7	48.8	49.3
	Bladder wall
	D max (Gy)	77.2	77.8	77.5	78.0	76.9	77.6
	D mean (Gy)	27.6	26.8	28.4	27.5	27.2	26.2
	D 25% (Gy)	44.8	41.8	46.1	43.1	43.5	40.6
	D 50% (Gy)	20.7	19.2	21.5	19.6	20.3	18.7

Table 4 Monaco Pinnacle Coronal Sagittal Coronal Sagittal

 4

	γ γ γ γ index test: criterion 3% / 3mm, dose level: 30%
	% accepted point	100	98.5	98.2	97.4
	γ γ γ γ mean	0.24	0.31	0.38	0.40
	γ γ γ γ max	0.78	1.22	1.37	1.33
	Dose gradient
	Mean gradient (cGy/mm)	4.6	4.6	4.3	5.1
	Maximum gradient (cGy/mm)	25.6	28.3	29.4	32.8
	Efficiency
	Delivery time (s)	169 [130-237]	165 [139-203]
	MU	688 [477-902]	452 [322-614]
	Mean Dose Rate (MU/min)	230		160
	Mean Area (cm²)	2.6 [2.0-3.9]	4.1 [2.6-6.7]
	Minimum Area (cm²)	0.5 [0.2-1.0]	1.3 [0.7-3.0]
	Maximum Area (cm²)	4.9 [3.4-7.7]	6.3 [4.1-8.6]

Concerning QA analysis, many precautions have been taken. First, to limit the impact of the beam-modeling, one set of beam data was used for the commissioning of both TPS [15;21].

Previous QA results showed that TPS commissioning are consistent for Monaco and Pinnacle.

We found similar dosimetric agreement for static fields for both TPS: for 3×3 cm² and 20×20 cm² field sizes the mean gamma values were 0.208 and 0.238, and 0.163 and 0.265 in a homogeneous phantom for Monaco and Pinnacle, respectively.

Then to avoid uncertainties with QA implementation, the same operator performed the measurements simultaneously for both TPS on the linear accelerator. Our results demonstrated a high passing rate for VMAT plan QA for both TPS. Nevertheless, Monaco provided better dosimetric agreement than Pinnacle. However, the lowest dosimetric agreement results were obtained for the planes where the highest dose-gradients were observed (i.e. sagittal Pinnacle planes).

In addition to this data , film analysis in an anthropomorphic phantom showed that for both TPS more than 98 % of points passed the (3%, 3 mm) gamma index criterion. Therefore, the dose algorithm accuracy of both TPS was verified in inhomogeneous conditions for pelvis cases thus proving that that there is no advantage to use Monaco at the expense of Pinnacle regarding to dosimetric accuracy for the studied case. Shrink Margin = 0 4920 -3420 cGy