R Ruffilli Ac

M Berkani

P Dupuy

S Lefebvre

Y Weber

B Warot-Fonrose

C Marcelot

M Legros

Aluminum metallization and wire bonding aging in power MOSFET modules

Keywords: Power MOSFET-based device, Metallization microstructure aging, Scanning electron microscopy (SEM), Ion milling, Ion microscopy, Focused ion beam (FIB) tomography, Transmission electron microscopy (TEM), Grain structure mapping

A limiting factor for the long-term reliability of power MOSFET-based devices is the electro-thermal and/or thermo-mechanical aging of the metallic parts. In this paper we assess the bonding wire and source metallization degradation of power devices, designed for applications in the automotive industry. Our approach consists in characterizing the metal microstructure before and after accelerated aging tests, by scanning electron microscopy, ion milling and microscopy, focused ion beam tomography, transmission electron microscopy and grain structure mapping. To focus on the wire-metallization bonding interface, we have set up a dedicated sample preparation that allows us to disclose the metallization under the bonding wires. This critical location is significantly different from the naked metallization, as the bonding process induces plastic deformation prior to aging. The main mechanism behind the device failure is the generation and propagation of fatigue cracks in the aluminum metallization. Away and under the wire bonds, they run perpendicularly from the surface down to the silicon substrate following the grain boundaries, due to an enhanced self-diffusion of aluminum atoms. Moreover, initial imperfections in the wire-metallization bonding (small cavities and aluminum oxide residues) are the starting point for harmful cracks that propagate along the wire-metallization interface and can eventually cause the wire lift-off. These phenomena can explain the local increase in the device resistance occurring at failure.

Introduction

The development of hybrid, mild hybrid and electric cars along with the integration of smart technologies have boosted the power electronics in automotive applications. Modern power MOSFET devices integrate several advanced capabilities in the field of security, safety, comfort etc. At the same time, during their everyday life the components have to support high current densities, because of the continuous shrinkage of new silicon technologies. The risk of failure is therefore physically increased, but has to be technically restrained to almost zero risk. In fact, power devices have to satisfy new robustness requirements to withstand the stress of repeated electro-thermal cycles for all the vehicle's life. In this context, the comprehension of the failure mechanisms under electro-thermal stress is a key point in the reliability assessment.

Because accidental failures due to intrinsic defects are now almost eliminated at the fabrication step, the most critical parts of power components are the metallic ones [1][START_REF] Glavanovics | Proc. 30th Eur. Solid-State Circuits Conf[END_REF][START_REF] Pietranico | [END_REF][4]. These parts are the source metallization (generally made out of aluminum) on top of the silicon wafer, and in particular the metallization beneath the bond wires, because it is the location of the highest concentration of current and temperature and because these parts age systematically. During the on-off operating cycles, the temperature changes result in mechanical stress due to the difference in the coefficients of thermal expansion (CTE) between the aluminum (Al) and the silicon (Si). Since aluminum has a low mechanical yield stress, the source metallization is expected to deform plastically even at moderate temperature excursions (200 °C), whereas silicon remains in its elastic domain. The result is the wellknown phenomenon of Al surface reconstruction, which leads to an increase in the drain-source resistance (R dson ) over cycles, until failure [5][6][7][8]. Whereas this reconstruction can lead to catastrophic failure has still to be determined.

In a previous study [9], we have shown that power MOSFETs aging occurs through a change in the Al grain microstructure and that plastic deformation of the top metal is due to an enhanced self-diffusion of the Al atoms along the grain boundaries [10,11]. We have also compared the initial Al metallization microstructure away from the bonding area and just under the bond wires: the latter is severely deformed, due to the bonding process. Moreover the wire-metallization interface shows initial imperfections, which probably are the starting point of crack propagation during aging.

This proves that the welded zone between the wire and the metallization is already deformed by the bonding process before aging. The thermomechanical stress, induced by the temperature cycles, weakens this interface and can lead to catastrophic events, such as bond wire lift off. Several studies focus on the bond wire failure mechanisms [START_REF] Glavanovics | Proc. 30th Eur. Solid-State Circuits Conf[END_REF]4,[12][13][14]. Bond wire heel crack failures [15] observed by experiments are found to be strongly dependent on the geometry and the thermal management of power devices. However, the fine characterization of the whole wiremetallization microstructure, which can provide useful feedback about the welding process, is poorly stated.

In the present study, we focus on the whole wire-metallization of power MOSFETs during accelerated aging tests. We characterize the microstructural evolution of the metallization just under the bonding wires, following a physical metallurgy approach: through scanning electron microscopy (SEM), focused ion beam (FIB) milling and microscopy, FIB-tomography, transmission electron microscopy (TEM) and grain structure mapping. These techniques are generally used to investigate the surface reconstruction phenomena of the naked metallization. In this study, thanks to a dedicated sample preparation that allows us to preserve the integrity of the bond interface, we could exploit them to follow the fine effects of the electro-thermal aging at this critical location.

Experimental procedure

Device under test and accelerated aging tests

The device under test (DUT) is a 12 V SmartMOS ® module from NXP Semiconductors designed for car applications. The device includes a power die, the active zone, and a control die, with diagnosis and protection functionalities (Fig. 1a). The power die consists in five MOS sectors, connected by eight wedge bond wires made out of Al. Here, we focus on the Al source metallization (3.6 µm thick), together with the Al bonding wires (381 µm of diameter) attached to it. The Al metallization contains 0.5 wt% of copper (Cu) and tungsten (W) and is coated by a silicon oxide (SiO 2 ) passivation layer (0.5 µm thick) everywhere, except in the bonding area to allow the electrical contact (Fig. 2a).

We performed specific aging tests (according to the AECQ100-12 standard [START_REF] Kelly | Automotive Electronics Council[END_REF]) at 25 °C and 70 °C, putting the component under repetitive short-circuit conditions until it failed. The supply resistance and inductance are fixed at R = 10mΩ and L = 5µH, which correspond to an output of Rshort = 50mΩ and Lshort = 5µH. During the shortcircuit, the device detects an overload and shuts down (fig 1b). Energy is dissipated at this stage, leading to an increase in the device temperature up to 250 °C-300 °C and affecting the device aging conditions.

The operating lifetime is considerably lower for the components tested at higher temperatures (5.3million cycles at 25 °C, 300k cycles at 70 °C). During the tests, only one sector undergoes repeated electrical pulses, so we could use the others as reference for the microstructural characterization of the metallic parts.

Microstructural characterization

Preliminary surface observations were performed in a Helios 600 FIB/SEM from FEI TM . Fig. 2b and2c show the comparison between the Al metallization (close to the wire bonds, where the passivation layer is absent) before and after aging, confirming a heavy surface reconstruction upon aging. The in-depth structure of the metallization away from the bonding area was observed by making cross-sectional cuts in the FIB/SEM microscope. These cross-sections are observed using ion channeling contrast, which depends on the grain orientation [START_REF] Ishitani | [END_REF]. Fig. 3 shows the comparison between a non-aged module and an aged one and reveals a strong degradation of the aged metallization, both bare and passivated, with many cavities and cracks running from the surface down to the transistor region.

However, cross-sectional observation of the Al power metallization under the bonding wire requires a preliminary sample preparation of the device: this is because it's unfeasible to cut through the entire wire thickness by ion beam. At the same time, a mechanical polishing could induce mechanical damage of the wire-metallization interface while it is the object of the investigation, and also cause the wire lift off that is usually poorly attached in the case of an aged device,. This preparation (Fig. 4) combines a mechanical grinding of the bonding wires at 45°, in order to locally reduce the wires thickness, and a final ion polishing by FIB. The cross-section ion images at the wire-metallization interface (Fig. 5) give us a first idea about the quality of the bonding process. It also allows to measure the initial deformation imposed to the Al layer by the wire bonding process and investigate the behavior of this interface during aging. High-resolution 3D images are acquired in a Helios 600 FIB/SEM, equipped with the Auto Slice and View software from FEI TM . The FIB-tomography consists in the following steps: 1) preparation of a parallelepiped of user-defined size, 2) creation of a reference mark, close to the area of interest, for the automated drift correction, 3) serial sectioning and acquisition of a SEM images stack. The sample is placed at the eucentric position, where the ion and the electron beam converge with an angle of 52° (Fig. 6). While the ion beam is used to cut and polish serial thin slices, the relative images are acquired by the electron column. The result is a stack of electron images at the wire-metal interface, which allows to follow the crack propagation inside an aged module volume. In our experiments we collected a stack of 200 SEM images with a z-spacing of 50 nm.

Finally, thin lamella (Fig. 7a, 7b) are FIB-prepared for the in-depth investigation of the source metallization and the wire-metal interface by TEM. Electron energy loss spectroscopy (EELS) are carried out in a Tecnai TM F20 microscope, in order to investigate the chemical composition of the wire-metal interface. TEM observations are also performed in a CM20FEG equipped with a Digistar TM -ACOM unit from NanoMEGAS TM . This system acquires the diffraction patterns of a TEM sample and returns the crystallographic orientation of each scanned point to build orientation maps as in Fig. 7c.

Results and discussion

Al metallization away from the bonding area

Outside of the bonding area, the initial Al metallization is characterized by a "bamboo" grain structure (Fig. 3a, 3c): the grain size is on same order of the metallization thickness and most of the grain boundaries run perpendicular to the surface down to the Si substrate.

FIB images of aged modules (Fig. 3b,3d) show that the repetitive electro-thermal cycles cause the formation of cavities and cracks running from the surface to the transistor region. They follow the grain boundaries, confirming that aging is linked to an enhanced self-diffusion of Al atoms along the boundaries [10,11]. Moreover aging seems to lower the average grain size but, as we also observed in a previous work [9], ACOM mapping (Fig. 8) shows that most of the grain domains revealed by FIB have a small misorientation (less than 10°). This means that they are subgrains formed by the eraly gathering of pre-existing dislocation inside a single initially large bamboo grain, and the apparent size reduction is only a FIB artifact.

Under the passivation, the Al diffusion process along the grain boundaries is supposed to be cut off or at least slowed down [10,[START_REF] Kraft | Dislocation plasticity in thin metal films[END_REF]. Fig. 3d shows that plastic deformation in the passivated Al metallization is less homogeneous than in the passivation-free metallization: some zones are not aged, whereas in others we observe a weak concentration of cracks and cavities.

Al metallization under the bonding wires

Under the bonding wire, the Al metallization has endured plastic deformation due to the welding process. Fig. 5 and 7c shows that the initial grain size is systematically smaller in these zones, compared to the bamboo structure (Fig. 3a,3c) of the naked metallization. One can also notice that grains in the Al wire are much larger. Plastic deformation is highly uneven. Perpendicularly to the wire-metal interface (along the applied stress of the wire bonding), this deformation ranges from 2% to 27% in average in the measured sections. To obtain these values, we have used the measurement method illustrated in Fig. 9. We measure the deformation with respect to a reference (red line) by computing the metallization areas under (zone -) or over (zone +) the reference, which corresponds to the initial position of the metallization surface. Then the deformation is computed as:

DeformedAreaRatio = DEF Area REF Area (1)
where:

DEF Area = i A i ∑ ( 2
) REF Area = cut l * ref h (3)
This measurement method takes in account the fact that the initial deformation is not uniform along the wiremetallization interface, that the interface is not straight, and that extrusion or compression of the metallization occurs during the bonding process. We can clearly see it in Fig. 10, showing two closed cross-sections of the Al metallization under the same bonding wire: in the first one (Fig. 10a) the wire-metal interface is almost straight without deformation, whereas in the second one (Fig. 10b) a severe plastic deformation significantly changes the interface profile, provoking a thinning of the metallization on the sides of the micrograph and a bulging in the middle. In the measured areas, the maximum deformation peaks to about 36%, both in depression or extrusion. Average deformation is about 10%, no matter the direction of deformation.

The surface impacted by the wire is about 0.25 mm 2 , which, taking in account the pressure applied during the ultrasonic process, leads to an average compressive stress of 35 MPa. This stress is much smaller than the yield stress of a commercially pure Al with a grain size in the range of 1-2 µm, typically in the range of 100-150 MPa [19]. To explain that significant plastic deformation occurred under this small load, we have to assume that only a fraction of the 0.25 mm 2 of the metallization surface was truly in contact with the wire. A very coarse estimate would lead to a third of the metallization being impacted by wire-bonding plastic deformation. This also explains why the surface is deformed very inhomogeneously.

Upon aging, we locally measured the deformation by taking the initial metallization thickness away from the bonding wire (in a non-aged device) as reference. The amount of plastic deformation does not change significantly, in line with previous results [9]. Fig. 12. EELS spectra from the wire (in blue), the metallization (in green) and the wire-metallization interface (in red). The red spectrum at the wire-metallization interface reveals the presence of oxygen.

Wire-metallization interface

The bonding between the metallization and the wire is characterized by initial imperfections, visible in white contrast in fig. 10 and11a. Fig. 12 shows the EELS spectra from the wire (in blue), the metallization (in green) and the wire-metallization interface (in red). The spectrum at the wire-metallization interface reveals the presence of oxygen, proving that the initial imperfections are linked to the presence of Al oxide. The most probable explanation for this presence is that this native oxide is not fully broken during the bonding process.

During aging, cracks propagate along these imperfections (Fig. 11b), causing a reduced contact between the wire and the metallization or, at worst, the wire lift-off. Cracks also run perpendicularly to the interface under the wires, as in the metallization away from the bonding area.

Crack propagation upon aging is investigated by FIB-tomography experiments. In Fig. 13, a series of six (non consecutive) SEM images from the wire-metallization interface of a 70 °C aged module are showed. In the analyzed area, the contact between the wire and the metallization is severely affected by the electro-thermal aging (image n° 1). A measure of the horizontal delamination surface using the complete stack of tomographic SEM sections leads to a delamination ratio of 78% over the tested volume. If we assume that about 33% of the wire was attached to the metallization (see paragraph 3.2), then delamination affected another 10% of the expected contact during aging. Of course, the volume tested by FIB is too small to make robust statistics, but we have gathered convergent observations showing that a significant delamination occurs between the wire and the metallization it's attached to. Going through the metallization volume, one can note some cavities appear (image n° 43) and subsequently constitute a single crack running from the bond interface to the transistor area (image n° 118). This crack splits in two new cracks (images n° 135 and 166) that propagate in the Al metallization (image n° 176). The high density of metal breaks revealed by the tomography stack can explain the increase of the R dson , which implies that subsequent electrical pulses will generate more heat, thus more deformation, and crack propagation.

The same aging mechanisms are observed in devices tested at 25 °C and 70 °C. Fig. 14 shows the comparison between two cross-section images of a 25 °C and 70 °C aged metallization under the bonding contact. In both cases, cracks propagate horizontally along the wire-metallization interface and perpendicularly to the surface down to the transistor region. However a higher test temperature can accelerate this behavior, explaining the lower life-time of the components aged at 70 °C. In these tests, the self passivation of the Al layers are observed as the thickness of the initial oxide keeps increasing upon aging, which is in agreement with the observation of oxygen content increase in the metallization of diodes undergoing electro-thermal aging [20]. However, contrary to what has been observed by these authors, this growth does not prevent the propagation of cracks.

Conclusions

The electro-thermal aging of the Al-based source metallization of power MOSFETs has been assessed by systematic observations of its microstructure.

Away from the wire bonds, the metallization is characterized by an initial "bamboo" structure (with the grain size of the same order of the metallization thickness). During aging, it undergoes a limited plastic deformation (apparent grain shrinkage, due to FIB artifacts) but several cracks are observed to run perpendicular to the surface down to the Si substrate, following the grain boundaries.

Under the wire bonds, the initial microstructure is finer than elsewhere because of the plastic deformation imposed by the bonding process. Upon aging a strong grain growth is observed. Moreover, cracks propagate not only perpendicularly to the surface, but also horizontally, along the wire-metallization interface. These cracks originate from the residual aluminum oxide, not fully broken during the bonding process. Because of this extensive horizontal crack propagation in aged devices, the remaining contact area between the wire and the metallization is reduced to a portion of the expected contact area, which appears as a probable source of device failure.

The aging of the Al metallization, both away and under the wire bonds, can be associated to a heavy generation of cracks that, cycle after cycle, are bound to propagate inside the material because of the auto-passivation of the Al around them. Two main mechanisms have been identified behind this behavior: 1) an enhanced self-diffusion of the Al atoms along the grain boundaries and 2) a poor bonding of the wires. In both cases, the result is an increase in the R dson of the devices during electro-thermal aging, leading to failure.

Fig. 1 .

 1 Fig. 1. (a) 12 V NXP's SmartMOS ® including the power and the control die. (b) Typical signals during a short cycle event in the power MOSFET by activation of overload protection.

Fig. 2 .

 2 Fig. 2. (a) SEM image of an Al bond wire welded on the Al source metallization. Away from the bonding area the metallization is coated by a SiO2 passivation layer. At higher magnification, bare metallization (close to the bond wire) of an (b) as-is module and (c) an aged one.

Fig. 3 .

 3 Fig. 3. FIB imaging of the bare metallization of an (a) as-is and (b) aged (70 °C) power module and of a SiO2 coated metallization in an (c) as-is and (d) aged (70 °C) power module. The red arrows highlight the SiO2 passivation layer. Scale bar 2.5 µm.

Fig. 4 .

 4 Fig. 4. (a) Bonding wire welded on the source metal of a power module. (b) The wire is cut and polished at 45° before being cut by FIB for interface observation. Scale bar 100 µm.

Fig. 5 .

 5 Fig. 5. FIB imaging of the wire-source metallization interface in an as-is power module. Scale bar 2.5 µm.

Fig. 6 .

 6 Fig. 6. Schematic illustration of the serial sectioning procedure by dual-beam FIB.

Fig. 7 .

 7 Fig. 7. (a) SEM image of a FIB-prepared lamella from an as-is module and (b) TEM detail of the wire-metal interface, with (c) the relative ACOM orientation maps.

Fig. 8 .

 8 Fig. 8. ACOM maps of the naked metallization of an aged module (70 °C).

Fig. 9 .

 9 Fig. 9. Schematic illustration of the measurement of the plastic deformation imposed to the Al source metallization by the bonding process.

Fig. 10 .

 10 Fig. 10. FIB images of the wire-metallization interface in two closed areas under the same bonding wire. The images refer to a non-aged module. Scale bar 2.5 µm.

Fig. 11 .

 11 Fig. 11. FIB images of the wire-metallization interface in (a) an as-processed module and in (b) an aged one (at 25 °C). Scale bar 2.5 µm.

Fig. 13 .

 13 Fig. 13. Selection of six non consecutive SEM images from the data stack of the wire-metallization interface of a 70 °C aged device (slicing distance 50 nm). Scale bar 2.5 µm.

Fig. 14 .

 14 Fig. 14. FIB images of the wire-metallization interface of a device aged at 25 °C and 70 °C. Scale 2.5 µm.


Jean-Emmanuel Pierré

Jean-Charles Passieux
email: passieux@insa-toulouse.fr

Jean-Noël Périé

J.-C Passieux

J.-N Périé

Finite Element Stereo Digital Image Correlation: framework and mechanical regularisation

Keywords: Digital Image Correlation, Stereovision, Finite Element, Mechanical Regularisation, Plate kinematics

The use of Finite Element meshes in Digital Image Correlation (FE-DIC) is now widespread in experimental mechanics. Up to now FE have been much less used in Stereo-DIC. The first goal of this paper is to explain in details how to use FE in Stereo-DIC using a formulation in the physical coordinate system. More precisely, it is shown how to perform the calibration of possibly nonlinear models, shape and displacement measurement based on a FE mesh. In addition it is shown that with such a framework it is possible to regularise the measurement with a FE model based on the same mesh. For instance, using this technique, it is shown that it is possible to measure the rotation field of a bending plate in addition to its displacement.

Introduction

Digital Image Correlation (DIC ) and Digital Volume Correlation (DVC ) [START_REF] Horn | Determining optical flow[END_REF][START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF][START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF] are now widely used in experimental mechanics. This family of methods aims at measuring 2D or 3D displacement fields from digital planar or volume images, respectively. In their initial version, the idea was to search for the best parameters of a given transformation to register small independent subsets of pixels. In the context of experimental mechanics, the use of Finite Elements in the context of DIC/DVC (FE-DIC/DVC [START_REF] Sun | Finite element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-Element" Displacement Fields Analysis from Digital Images: Application to Portevin-Le Châtelier Bands[END_REF]) has been extensively developed during the last decade. Among many advantages, to us, the main interest relies on the fact that it provides a very simple bridge between experiments and simulations, which is very convenient for validation [START_REF] Passieux | Local/global non-intrusive crack propagation simulation using a multigrid X-FEM solver[END_REF], identification purposes [START_REF] Leclerc | Integrated Digital Image Correlation for the Identification of Mechanical Properties[END_REF][START_REF] Roux | Digital image correlation and fracture: an advanced technique for estimating stress intensity factors of 2d and 3d cracks[END_REF][START_REF] Passieux | Multiscale Displacement Field Measurement Using Digital Image Correlation: Application to the Identification of Elastic Properties[END_REF] and for performing mechanically regularised measurements.

Even in situations where 2D-DIC could be used, Stereo-DIC is preferred because of both accuracy and simplicity of the setup [START_REF] Garcia | 3d deformation measurement using stereo-correlation applied to experimental mechanics[END_REF][START_REF] Sutton | The effect of out-ofplane motion on 2d and 3d digital image correlation measurements[END_REF]. Classically, the idea of Stereo-DIC is to use one stereoscopic image pair f 1 and f 2 in a reference configuration at t 0 and a second one g 1 and g 2 at t 0 + dt in a deformed configuration, see Fig. 1. 1) find ex/in-trinsic parameters from a set of images of calibration grids Shape Measurement:

2) find matching field d between reference images f 1 and f 2

3) triangulation: find X at t 0 from x 1 and x 2 = x 1 + d(x 1 )

Displacement Measurement: 4) find matching t between images f 1 and g 1 5) find matching t between images f 1 and g 2 6) triangulation: find X at t 0 + dt from x 1 = x 1 + t(x 1 ) and

x 2 = x 1 + t (x 1 )
Eventually, the 3D displacement is estimated by comparing measured shapes:

U = X -X.
The drawback of this classical Stereo-DIC formulation is that the displacement is not the unknown of a unique optimization problem as it is the case with 2D-DIC and DVC. In this situation, the bridge between measurement and simulation is less direct. In particular, applying a mechanical regularisation is not possible in this framework.

In opposition to DIC and DVC, the literature on FE-Stereo-DIC (or FE-SDIC) is much less prolific certainly because the extension is not straightforward. The first attempt was proposed in [START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3d kinematics[END_REF] where a theoretical FE mesh was used to perform the displacement measurement solution of a unique problem defined in the world coordinate system. In this paper this formulation is not fully exploited since it is not used for calibration and shape measurement. A very similar global stereo formulation based on IsoGeometric Analysis (IGA) was proposed in [START_REF] Beaubier | CADbased calibration and shape measurement with stereoDIC: Principle and application on test and industrial parts[END_REF][START_REF] Dufour | CAD-based Displacement Measurements with Stereo-DIC: Principle and First Validations[END_REF]. It is used not only for the displacement measurement, but also for calibration of a linear camera model and shape measurement.

In this paper, we develop and detail a general framework covering all the aspects of FE-Stereo-DIC. It consists in a formulation in the physical coordinate system which makes possible to calibrate a stereo rig (with possibly non-linear camera models), to measure the actual shape of the specimen and to measure the displacement field, using a FE interpolation. A technique is proposed to build the minimal number of integration points in this context. In addition, it become possible to add a regularisation based on a FE mechanical model. For instance, the method is illustrated with a measurement of the kinematic fields of a bending plate; the kinematic fields meaning in this case not only displacements but also rotations.

In Sect. 2 the three problems of calibration, shape measurement and displacement measurement will be rewritten as three optimisation problems where the unknowns are the intrinsic/extrinsic parameters, the shape and 3D displacement respectively. Then, in Sect. 3, it will be shown that, with the proposed framework, it is directly possible to add a constraint of regularity such as a mechanical regularisation term. It will be also shown that the use of a plate model in the regularisation, is likely to provide accurate rotation fields. Finally in Sect. 4, both synthetic and real test case will be presented in order to illustrate the efficiency of the method, the effects of the parameters and its robustness with respect to noise.

2 Rewriting Stereo-DIC problem in the world reference system with

Finite Elements

This section contains the three parts of the measurement of a 3D displacement in the reference system R W described by the mesh. Indeed, with an FE mesh, the theoretical shape does not necessarily correspond to the real one. Therefore, a displacement measurement needs two preliminary steps (Fig. 2). Firstly, cameras have to be located in the world reference system. This is called the extrinsic parameters calibration but intrinsic parameters can also be optimised ac-cording to the used mesh Sect. 2.1. Then, the mesh shape has to be corrected: the projection of each points of the mesh in both cameras (with their projectors P c ) have to be stereo-matching (the corresponding gray-levels must be equals). This step is the shape measurement Sect. 2.2. A very similar approach has already been published in [START_REF] Beaubier | CADbased calibration and shape measurement with stereoDIC: Principle and application on test and industrial parts[END_REF] for a CAD-based shape measurement with NURBS. Finally, the 3D displacement between two deformation steps can be measured Sect. 2.3.

Extrinsic and Intrinsic parameters calibration

In this study, the method being based on a Gauss Newton algorithm, the calibration has to be initialised with the Stereo Software Vic-3D TM . The calibration consists in finding all parameters p c of each camera's projector P c :

-6 extrinsic parameters in order to find the position of the camera's reference system in R W : 3 translations and 3 rotations.

-4 linear intrinsic parameters in order to project a point from the camera to its image : 2 focal lengths (f u , f v ) and two centers of image (c u , c v ) along both horizontal u and vertical v directions.

at least 2 non-linear intrinsic parameters : the first order of radial distortion κ and the skew angle.

More details on the non-linear models of cameras can be found in [START_REF] Garcia | 3d deformation measurement using stereo-correlation applied to experimental mechanics[END_REF][START_REF] Pierré | Unstructured finite element-based digital image correlation with enhanced management of quadrature and lens distortions[END_REF]. In this paper, the intrinsic parameters are initialised with Vic-3D TM . The first problem consists in finding the extrinsic parameters in order to locate the 3D Finite Element mesh (considered rigid) as close as possible of the real surface which is unknown.

In the following, p c denote the extrinsic parameters of the projector P c associated to the image f c (the intrinsic ones are fixed). For the sake of simplicity, we first consider a standard stereo bench, i.e. c = {1; 2}. The parameters p c are assessed by minimizing the following objective functional based on the gray-level conservation assumption:

p 1 , p 2 = arg min p 1 ,p 2 Ω f 1 (P 1 (X, p 1 )) -f 2 (P 2 (X, p 2 )) 2 dX (1)
where X stands for a 3D point located on the visible surface Ω of the specimen.

As said previously, the quadrature is done in the 3D domain and an image is not preferred over another (there is no concept of master-slave), each camera being treated symmetrically. This non-linear problem is solved with a Newton algorithm.

At iteration k, parameters p = [p 1 , p 2 ] T are searched in the form p k+1 = p k + δp.

After linearisation and differentiation, the problem reads:

M k ext δp = b k ext          M k ext,ij = Ω ∂P 1 ∂p i ∇f 1 - ∂P 2 ∂p i ∇f 2 ∂P 1 ∂p j ∇f 1 - ∂P 2 ∂p j ∇f 2 T dX b k ext,i = - Ω ∂P 1 ∂p i ∇f 1 - ∂P 2 ∂p i ∇f 2 f 1 P 1 (X, p k 1 ) -f 2 P 2 (X, p k 2 ) dX
where ∂P c ∂p i are the projectors P c gradient with respect to the extrinsic parameters (even with a non-linear model, an analytical expression is known). ∇f c is image f c gradient at P c (X, p k c ) which depends of p c . This means that the operator M k ext has to be re-assembled at each iteration.

For a Gradient algorithm, the initial guess should be close to the solution. The first guess can be obtained "by hand" by selecting some points in each image and on the mesh (it would also be possible to detect at least three particular points with a pattern recognition algorithm). Then, a small non-linear problem has to be solved. This problem has the same expression than the bundle adjustment (2) [START_REF] Triggs | Bundle adjustment-a modern synthesis[END_REF], but here, X, x 1 , x 2 = x 1 + d(x 1 ) are fixed and only the extrinsic parameters are searched for:

p 1 , p 2 = arg min p 1 ,p 2 P 1 (X, p 1 ) -x 1 2 2 + P 2 (X, p 2 ) -x 1 -d(x 1 ) 2 2 (2)
In order to illustrate these steps, on Fig. 2(a) can be seen the first guess of the position of the mesh on the real surface. Fig. 2(b) is after the extrinsic calibration (eq. ( 1)). In these two steps, the mesh X is considered fixed in its reference system and the goal is to position this system as close as possible to the unknown surface by finding the rigid translations and rotations.

In paper [START_REF] Beaubier | CADbased calibration and shape measurement with stereoDIC: Principle and application on test and industrial parts[END_REF], some intrinsic parameters of linear camera models are estimated using the same framework (using the speckle pattern of the specimen only). But the measurement of a least one distance is needed to get an absolute estimate of, for instance, the focal length. Here, the model of camera is more complex (non-linear with distortions) and the initial mesh is flat. Thus, intrinsic parameters cannot be optimised at the same time than the extrinsic ones. The idea would be to release all intrinsic parameters once the extrinsic are optimized. The functional is the same (eq. ( 1)) but with p c containing only intrinsic parameters (c u , c v , f u , f v , κ, skew)

(the extrinsic ones are now fixed). It is possible to alternate these optimizations.

In practice, and in order to be more generic, we decided to calibrate the intrinsic parameters using classic planar calibration targets. This step was done using Vic-3D TM . It was verified that these parameters were very close to the local minimum of the above graylevel based functional.

Shape measurement with an FE mesh

Once all cameras' parameters p c are known (and optimised), the shape has to be corrected. Measuring the shape means finding the position X verifying the graylevel conservation. Here, the latter is verified between the reference images of each camera. The same method as for the cameras' parameter is used but the unknowns are the 3D positions X:

X = arg min X Ω f 1 (P 1 (X, p 1 )) -f 2 (P 2 (X, p 2 )) 2 dX (3)
At iteration k, the estimation of the position is written X k+1 = X k + δX where the displacement correction δX is searched in the finite element subspace with

δX = i N i (X)q i : M k shape q = b k shape        M k shape,ij = Ω N T i (∇P 1 ∇f 1 -∇P 2 ∇f 2 ) (∇P 1 ∇f 1 -∇P 2 ∇f 2 ) T N j dX b k shape,i = - Ω N T i (∇P 1 ∇f 1 -∇P 2 ∇f 2 ) f 1 P 1 (X k , p 1 ) -f 2 P 2 (X k , p 2 ) dX
where ∇P c are the projectors' spatial gradient which have an analytical form. As for the calibration, the operator M k shape has to be re-assemble at each iteration.

After convergence, an estimation of the real shape is obtained (see Fig. 2(c)).

Remarks:

-Since this is based on a Gauss Newton, this algorithm requires an initialisation step. Usually, a multigrid initialisation with a pixel aggregation (coarse graining) is used [START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF]. But unlike what is usually done in 2D-DIC, where this step is made with a coarse mesh, here the same mesh is used with a decreasing Tikhonov regularisation in order to avoid 3D projection of fields from a mesh to another.

-As said previously, M k shape has to be re-assemble at each iteration. This has a significant cost because here, the reference images of all cameras are compared by pairs. For n cameras, there are n(n-1) 2 couple of cameras. On the other hand, in [START_REF] Dufour | Shape, displacement and mechanical properties from isogeometric multiview stereocorrelation[END_REF], a reference object f with an intrinsic texture is created. Thus, all images are compared to f which is faster when using four or more cameras because only n couples are considered.

Regularisation As always in DIC, this problem is ill-posed. With the classical approach, there are two unknowns (the components of the left-right disparity field d on the image) for only one scalar equation (gray-level conservation). Here, a 3D field is searched for, which means that for the same number of equation, there is one more unknown. Indeed, in the problem 3, the mesh can slide along the surface of the object. There can be a global sliding (see Fig. 3(b)) which imply that the dimension of the operator's kernel is at least 3 (two sliding along two surface directions and one rotation). But in the formulation, the problem is even harder because each node can moves relatively to another (see Fig. 3(c)) without changing the value of the functional. Thus, a higher degree of regularisation has to be added to the functional. For the local sliding, an isometric constraint can be imposed, like a Tikhonov regularisation term:

δX = arg min δX Ω f 1 (P 1 (X + δX, p 1 )) -f 2 (P 2 (X + δX, p 2 )) 2 dX + λ δX 2 K
where the norm • K can be associated to a mechanical operator K (truss, stiffness plate/shell or Laplace). One of the interest of using SDIC with NURBS [START_REF] Beaubier | CADbased calibration and shape measurement with stereoDIC: Principle and application on test and industrial parts[END_REF], rely on the fact that only a few control points are necessary to describe a surface which indirectly regularise the problem. Indeed, the shape interpolation can be uncoupled from the displacement interpolation which can be refined without modifying the shape. In FE, this uncoupling is not possible and the same mesh is used for both shape and displacement measurement. Thus, it must be regularised externally.

However, both NURBS and FE approaches are subject to global sliding. The previous regularisation does not prevent that problem. These shifts are particularly visible when iterating between the optimisation of the cameras' parameters and the shape measurement.

To avoid these shifts and restore the uniqueness of the problem, it has been decided to seek the correction of the nodal displacement along the normal to the surface [START_REF] Passieux | Quelques outils numériques pour la simulation et la mesure en mécanique des structures[END_REF]. The problem has then only one dof per node. In practice, the surface normal is estimated at each node (average of the neighbouring element normals) in the initial shape of the mesh.

3D displacement measurement

Once the real shape (thus, the real mesh) is known, it is possible to measure the 3D displacement U between two steps (t 0 and t 0 + dt). Unlike Classical Stereo-DIC, the three optimisation problems are turned into on single problem where U is the unique unknown. The main idea for this reformulation is to work only in the 3D mesh, chosen to be the world reference system R W . All problems are written in this 3D coordinate system and weak-form integrals of the gray-level conservation are also in R W .

For each 3D point X, the projector P c gives a point x in the image. Obviously, it is also possible to have the projection of the displaced point X = X + U(X) (at t 0 + dt) in the deformed state images g c (see Fig. 4). The point X is projected on the reference image f c and X is projected on the deformed state image g c . Their gray level should be equals and this, for each camera c. Thus, the functional reads:

U = arg min U c Ω f c (P c (X, p c )) -g c (P c (X + U(X), p c )) 2 dX (4)
which is just the gray level conservation equation for each camera. This correspond to the vertical pairing in Fig. 4. Each term of each camera being independent, it is possible to consider more than two cameras [START_REF] Dufour | Shape, displacement and mechanical properties from isogeometric multiview stereocorrelation[END_REF]. Such a formulation in the world reference system was first introduced in [START_REF] Réthoré | Robust identification of elasto-plastic constitutive law parameters from digital images using 3d kinematics[END_REF]. The benefit is to symmetrically use each camera. At iteration k, a correction of the 3D displacement field is introduced

U k+1 = U k + i N i (X)q i and the problem reads: c M c,k dic M k stereo q = c b c,k dic b k stereo x 1 Rw R1 f 1 f 2 x' x' 2 1 g 1 g 2 Rw t0+dt x2 U(X) X X'=X+U(X)
       M c,k dic,ij = Ω N T i (∇P c ∇f c ) (∇P c ∇f c ) T N j dX b c,k dic,i = Ω N T i (∇P c ∇f c ) f c (P c (X, p c )) -g c (P c (X + U k , p c )) dX
Remarks:

-Both left and right hand side are sums of independent operator on each camera.

This implies that this problem can be parallelized like what is done in domain decomposition method in 2D-DIC [START_REF] Passieux | A dual domain decomposition method for finite element digital image correlation[END_REF].

-Like shape measurement, this algorithm needs an initialisation step such as a multigrid initialisation based on a pixel aggregation (coarse graining) and a Tikhonov regularisation with a Laplacian operator.

-The correlation operator M k stereo depends on the displacement (because of the gradient of the projectors ∇P c (X + U k )) which means once again that it should be re-assemble at each iteration. In practice, this is not done because the multigrid initialisation is close to the actual displacement. Thus, no significant difference between re-assembling or not was found (except that it is far more time consuming).

To avoid a time drift of the pairing, it is also possible to introduce an additional term which minimize the difference of gray level between the deformed state images at t 0 + dt taken in pairs (this corresponds to an horizontal pairing in Fig. 4):

U = arg min U c Ω f c (P c (X), p c ) -g c (P c (X + U(X), p c )) 2 dX + c e<c Ω g c (P c (X + U(X), p c )) -g e (P e (X + U(X), p e )) 2 dX (5)
the pairing at t 0 is already verified with the calibration and the shape measurement. Of course, with this last term, a non-negligible cost is added when more than two cameras are used.

Quadrature. In all previous equations, quadrature is done in the 3D world reference system and not on the images as with classical DIC. A dedicated quadrature method has to be set up. Contrarily to simulation where a Gauss method is efficient, here, the integrand involves the grayscale gradient which has a very short variation length. Thus a Riemann integral method was therefore chosen by splitting homogeneously and homothetically each finite element. This choice was proven to be more accurate [START_REF] Pierré | Unstructured finite element-based digital image correlation with enhanced management of quadrature and lens distortions[END_REF]. In practice, a smaller number of points can be build with an equivalent accuracy. It considers an inhomogeneous partition of quadrangles and triangles, in the isoparametric coordinate system, build from homogeneous partition of the two smallest edges. On Fig. 5, we show that such a technique reduces significantly the number of integration points. According to our tests, using this method, their are rarely more than one integration point per pixel. Fig. 5 Example of a different quadrature rule in the isoparametric coordinate system in order to find as few integration points as possible Remark. In this article, this formulation in the 3D coordinate system is needed in order to add a mechanical regularisation term to the Stereo-DIC measure. But this is also interesting in 2D-DIC [START_REF] Pierré | Unstructured finite element-based digital image correlation with enhanced management of quadrature and lens distortions[END_REF] for two reasons:

-The pixel based quadrature does not integrate accurately a constant function.

The area of an element is approximated by an integer number of pixels and depends on the point of view of the camera. Even the shape of the elements are not accurate because of the linear approximation of the element edges: with distortions the projection of a segment is not necessarily linear. This has a weak impact on the standard uncertainty of the measure but can improve the systematic bias of an order of magnitude even with a linear simple model of camera.

-A more complex model of camera with distortions can be used which drastically decreases measurement uncertainties. This is usefull when identifying parameters with a FEMU method based on displacement (FEMU-U). It is also an essential method for crack propagation based on Williams series if it is not possible to use a telecentric lens.

Stereo-DIC and Plate Regularisation (R-FE-SDIC)

Once having the 3D displacement as the unique unknown, it is easier to add a regularisation term. The idea is not only to make a measurement, but also to filter the latter with a mechanical model. This is very important when trying to measure accurate rotation fields.

Furthermore, for a plate or shell model, not only displacements but also rotations at the mid-surface are considered (see Fig. 6). The 3D displacement U measured with Stereo-DIC is the displacement of the observed surface which is the top one. Here, we propose a method in order to measure directly the midplane unknowns.

Plate kinematic Regularisation

With a plate theory (such as the Mindlin-Reissner theory of plates), a simple linear projector Π can be used to determine the displacement of the upper skin depending upon the displacements V and rotations θ of the mid-surface. By extension and for a sake of simplicity, the vector of unknowns containing displacements and rotations of the mid-surface is denoted by V.

U = ΠV ⇔       U x U y U z       = Π              V x V y V z θ x θ y             
Then the Stereo-DIC problem can be written:

V = arg min V 1 -λ m α v c Ω f c (P c (X, p c )) -g c (P c (X + ΠV(X), p c )) 2 dX + λ m α m KV 2 2 (6)
where the unknown is not only the displacement of the upper skin U, but the five degrees of freedom V of the mid-surface. α v and α m are coefficients for the weighting of the terms of vision and stiffness. Thus, λ m is a dimensionless coefficient for the mechanical regularisation term set between 0 and 1 (see Sect. 4.1.3).

Finally, K is a regularisation matrix. Indeed, by adding two unknowns per node, a regularisation term is needed. But due to a formulation in the world reference system, the physical displacement is the unknown of the problem. Thus, it is easy to add a mechanical operator (such as KV = F for an elastic deformation) in order to regularise the functional. But here, for our study, it is not possible to have the external forces. Thus, the stiffness matrix K is replaced by K.

The latter regularisation matrix is defined by K = P k K with P k a diagonal matrix [START_REF] Lawson | Solving Least Squares Problems[END_REF]. The latter could be boolean, but in this study, P k is used to weight some dofs (eg to weigh the influence of rotations with respect to the membrane behaviour):

                               P ij k = 0 ∀j = i P ii k = 0
for all dofs on the non-free edges

P ii k = 0
for the membrane dofs related to V x and V y P ii k = 1 for the dofs related to V z

P ii k = 1 h 2
for the dofs related to θ x and θ y (h being the element size)

Thus, K is a matrix containing each dof except those supported by the nodes that are concerned with Dirichlet boundary conditions [START_REF] Réthoré | A fully integrated noise robust strategy for the identification of constitutive laws from digital images[END_REF].

At iteration k, the problem reads:

1 -λ m α v Π T M k stereo Π + λ m α m K T K M K δq = 1 -λ m α v Π T b k stereo - λ m α m K T K q 0
Remark. Replacing K by K means no regularisation on the degrees of freedom of the non-free edges. Thus the behaviour of the nodes of these edges is only towed by the vision term. But in this study, a plate model is used, which means that the unknowns are 3 displacements and 2 rotations at the mid-plate surface. There is an additional difficulty because each dof on the non-free edges implies a singularity in M K . On the one hand, if these dofs are not taken into account, the problem cannot be solved. On the other hand, if they are removed from M K , it means that the associated nodes are clamped. Thus, there is a difficulty that will be solved in the next section 3.2, especially for the rotations. Indeed, the displacement dofs can be considered in M k stereo but not the rotation dofs.

Mechanical and Tikhonov Regularisation

To overcome the absence of mechanical regularisation on the non-free edges (particularly for rotations of these nodes in question), it is possible to add a Tikhonov regularisation term (minimisation of a Laplacian term):

V = arg min V 1 -λ m -λ t α v c Ω f c (P c (X, p c )) -g c (P c (X + ΠV(X), p c )) 2 dX + λ m α m KV 2 2 + λ t α t ∆V 2 2 (8)
which gives:

1 -λ m -λ t α v Π T M k stereo Π + λ m α m K T K + λ t α t T T T δq = 1 -λ m -λ t α v Π T b k stereo - λ m α m K T K q 0 - λ t α t T T T q 0
where T is the Tikhonov matrix corresponding to the Laplacian term, α t is a weighting coefficient and λ t is a second dimensionless coefficient for this regularisation term set between 0 and 1 with respect to λ m (see Sect. 4.1.3). Concretely, each node excluded from K will be associated to a zero eigenvalue. Adding the Laplacian term amounts to impose that those nodes adopt a similar behaviour than their neighbour.

Remark. The normalisation coefficients α v and α m are computed with the initial guess of the displacement U 0 . For α t , the goal is to have the smallest Laplacian term, thus, it has been chosen to set it at 1. The coefficients are :

             α v = U 0 T M k stereo U 0 α m = U 0 T K T KU 0 α t = 1 4 Application

Synthetic test case

In order to quantify the measurement error with the described regularisation method, the goal here is to create some synthetic images for a plate test case. A pair of real images taken during a test is used (Fig. 8). The plate size is 26.5 × 26.5 cm 2 .

Both cameras being calibrated, the two projectors P c are known. On the initial plate, for each camera c, the 3D reference points of the plate X r = (X, Y, 0) T are also known. It therefore suffices to determine the deformed points (according to a plate model) corresponding to each reference point to create a new image. According to the small strain and small displacement assumption, each point originally located at position X r = (X, Y, 0) T will be deformed to the point

X d = (X, Y, p(X, Y ))
T with p a function giving the out-of-plane position of a point according to a chosen plate model. Thus, for each integer coordinated x (each pixel of the new image), two unknowns (X, Y ) are solution of a system of two equations that can be solved by a Newton algorithm (Fig. 7(a)):

P c ((X, Y, p(X, Y )) T ) -x = 0
Then, it is known that this deformed point was originally located at position

Mesh size influence

The idea here is to calculate the measurement error according to the mean size of elements of an FE mesh. Both FE-SDIC and R-FE-SDIC (regularised, here with a plate model) are used. Fig. 9 presents the standard deviation of the error with both method when measuring the out-of-plane displacement V z and the rotation θ y . For the latter, with an FE-SDIC method, only the 3D displacement of the upper skin U is measured and the rotations are calculated using the derivative of the out-of-plane displacement U z = V z .

It can be seen in Fig. 9(a) that when the mesh size decreases, the measurement error also decreases due to the "FE error" until a mean element size of approximately 10 mm. Then, because of a lack of pixels in the projected elements, the error increases, which corresponds to the "DIC error". In contrast, in blue, it can be noted that thanks to the mechanical regularisation, the error does not increase. Furthermore, in Fig. 9(b), it is shown that it is possible to directly measure the rotations regardless the mean element size. Furthermore, as expected, the derivative of a measured field (in red) is necessarily more noisy than a measured one (in blue).

Noise robustness

It is also possible to study the behaviour of measurement uncertainties when adding different level of noise when creating the deformed synthetic images cf. Fig. 10. For that, we added a zero-mean white Gaussian noise with a increasing standard deviation ranging from 1 to 8 graylevels. Thanks to the mechanical regularisation, the measurement is much more robust with noise.

Remarks:

-In this study, the image noise have been estimated to be less than 3 grey-levels (Allied Vision Pike FireWire 5 Megapixels).

-No filtering algorithm has been used in this work. The only regularisation is based on the use of a mechanical operator.

Regularisation coefficients

Obviously, the result from such a measurement depends on the used parameters λ m and λ t . They are the cut-off wavelength of low-pass filters [START_REF] Leclerc | Voxel-Scale Digital Volume Correlation[END_REF][START_REF] Leclerc | Digital volume correlation: what are the limits to the spatial resolution[END_REF][START_REF] Réthoré | Automatic crack tip detection and stress intensity factors estimation of curved cracks from digital images: AUTOMATIC CRACK TIP DETECTION AND SIF ESTI-MATION OF CURVED CRACKS[END_REF]. Usually, this kind of parameters are set between 0 and 1 (eq. ( 9)). Thus, the three terms of ( 8) are normalized with α v , α m and α t .

               0 ≤ 1 -λ m -λ t λ vision ≤ 1 0 ≤ λ m ≤ 1 0 ≤ λ t ≤ 1 (9)
Thus, different values of λ m can be used in order to study the measurement error. Without mechanical regularisation, a displacement can be observed but not a rotation. Thus, the rotation measurement with a low value of λ m corresponding to a low mechanical regularisation should not be less accurate than with a greater regularisation coefficient. But, as it can be seen in Figs. 11(a) and (c), it is hard to find the optimal parameter, especially when looking to the rotation. The simplest way to find it for such an approach, is to draw a L-curve [START_REF] Miller | Least squares methods for ill-posed problems with a prescribed bound[END_REF][START_REF] Lawson | Solving Least Squares Problems[END_REF] which represent both DIC and Mechanical residual. The optimal parameter will naturally be at the corner of the L-curve cf. Fig. 11(e).

On the other hand, the parameter λ t aims to take into account the degrees of freedom with a Dirichlet condition. Thus, the only condition seems to have a non-zero parameter: λ t = 0. But it can be seen in Figs.

Real test case

The same measurement is now carried out with real deformed images. As expected, the measurement of the rotation by R-FE-SDIC (mechanically regularised FE-Stereo-DIC with a plate theory) is more accurate than without mechanical regularisation, because the latter is the numerical derivative of a measured field cf. Fig. 12.

The main advantage of a Stereo method based on the world reference system is that it is possible to use a real mesh. Indeed, in a classical FE-SDIC method, the mesh has to be adapted to the measurement. Here, the goal is to adapt the measurement in order to realize a real dialogue between experiment and simulation.

The mesh chosen for a simulation can also be used for a measurement. Thus, for example, the mesh can be refined in the center for a potential default cf. Fig. 13.

In this figure, the shape measurement step is compared to the one computed by Vic-3D TM . It can be seen that the point cloud resulting from Vic-3D TM is close to the mesh except at the boundary where Vic-3D TM shape is not available.

With this kind of mesh, even the displacement field computed without regularisation would not be correct because of the mesh refinement at the center cf.

Conclusion

A full stereovision framework for measuring the 3D shape and displacements using a Finite Element mesh was developed. The formulation incorporates a gray-level based calibration step extended to non-linear camera models. It was shown that, in opposition to stereovision using NURBS functions [START_REF] Beaubier | CADbased calibration and shape measurement with stereoDIC: Principle and application on test and industrial parts[END_REF], the shape measurement is very ill-posed and requires additional regularisation. We used for that a type of isometric constraint. Thanks to this formulation in the world coordinate system, the same FE mesh can be shared between both numerical simulation and measurement. It is thus a good tool to perform quantitative comparisons between non-planar experiments and corresponding simulations that are most often based on Finite Elements [START_REF] Sztefek | Tensile stiffness distribution in impacted composite laminates determined by an inverse method[END_REF][START_REF] Sztefek | Nonlinear compressive stiffness in impacted composite laminates determined by an inverse method[END_REF]. In addition, a technique was proposed to build the minimal number of integration points in such a context.

As an application, the method was used for the measurement of the displacement and rotation fields of a plate in bending, using a mechanical regularisation based on the elastic stiffness operator associated to the FE mesh. In addition, stereo synthetic test cases are developed and analysed in order to exemplify the efficiency of the proposed method.

Fig. 1

 1 Fig. 1 Main steps of a 3D displacement field measurement of a surface with classical Stereo-Digital Image Correlation. Pairings are in green and triangulations in blue

Fig. 2

 2 Fig. 2 Steps of the new Stereo-DIC formulation with a Finite-Element mesh (FE-SDIC): (a) first guess of the theoretical CAD mesh plate (black) on the real (a priori unknown) plate (green) before any calibration or shape measurement. (b) Calibration of the extrinsic parameters by DIC (with fixed 3D points X). (c) Shape measurement (with fixed parameters pc). (d) Temporal measure of the 3D displacement U (from the black position of the mesh to the red one)

Fig. 3

 3 Fig. 3 Illustration of the ill-posed problem of shape measurement. (a) Theoretical position of the mesh after calibration and shape measurement. (b) Global or (c) local sliding during the shape measurement

FE model: V and θ observed surface h Fig. 6

 h6 Fig.6Representation of a 3D plate or shell of a thickness h. The 3D displacement U measured with an optical bench is on the observed surface whereas the plate or shell model considers displacements V and rotations θ at the mid-surface

X

  r = (X, Y, 0) T (Fig. 7(b)). It is then possible to create a synthetic deformed image by interpolating the gray level of the reference image f c at the pixel x by projecting the reference 3D point P c (X r ) (Fig. 7(c)).

Fig. 7 Fig. 8

 78 Fig. 7 Creation of a deformed synthetic image with a plate model p. (a) search for the 3D deformed point X d according to a plate model p which correspond to the pixel x. (b) projection of this deformed point to the initial shape in order to find the reference 3D point Xr. (c) interpolation of the reference image in order to associate a gray level with the pixel x of the created deformed image

Fig. 9

 9 Fig. 9 Standard deviation of the error when measuring the out-of-plane displacement Vz (a) and the rotation θy (b) for different sizes of mesh. This is computed with both methods: FE-SDIC (red) and R-FE-SDIC regularised measure with a plate model (blue)

Fig. 10

 10 Fig. 10 Standard deviation of the error when measuring the out-of-plane displacement Vz (a) and the rotation θy (b) as a function of the standard deviation of the noise. This is computed with both methods: FE-SDIC (red) and R-FE-SDIC regularised measure with a plate model (blue)

  11(b), (d) and (f) that this parameter does have an impact on the measurement. Indeed, the weight of this coefficient has an influence on the boundary. And since boundary displacement uncertainty has a tendency to spread in the mechanical regularised region, it has an effect on the solution[START_REF] Leclerc | Digital volume correlation: what are the limits to the spatial resolution[END_REF].

Fig. 14 .

 14 Fig. 14.

Fig. 11 Fig. 12 Fig. 13 Fig. 14

 11121314 Fig. 11 Schemes for studying the influence of λm (left) and λt (right). Standard deviation of the error when measuring the out-of-plane displacement Vz (a, b) and the rotation θy (c, d) for different values of λm and λt. This is computed with R-FE-SDIC regularised measure with a plate model (blue) and the value corresponding to the FE-SDIC method is plotted (red) for a comparison. L-curve (e, f) of both parameter in order to minimise both DIC and Mechanical residual

Acknowledgements This work was funded by the French "Agence Nationale de la Recherche" under the grant ANR-12-RMNP-0001 (VERTEX project).


P Delva

J Lodewyck

S Bilicki

E Bookjans

G Vallet

R Le Targat

P.-E Pottie

C Guerlin

F Meynadier

C Le Poncin-Lafitte

O Lopez

A Amy-Klein

W.-K Lee

N Quintin

C Lisdat

A Al-Masoudi

S Dörscher

C Grebing

G Grosche

A Kuhl

S Raupach

U Sterr

I R Hill

R Hobson

W Bowden

J Kronjäger

G Marra

A Rolland

F N Baynes

H S Margolis

P Gill

Test of special relativity using a fiber network of optical clocks

Phase compensated optical fiber links enable high accuracy atomic clocks separated by thousands of kilometers to be compared with unprecedented statistical resolution. By searching for a daily variation of the frequency difference between four strontium optical lattice clocks in different locations throughout Europe connected by such links, we improve upon previous tests of time dilation predicted by special relativity. We obtain a constraint on the Robertson-Mansouri-Sexl parameter |α| 1.1 × 10 -8 quantifying a violation of time dilation, thus improving by a factor of around two the best known constraint obtained with Ives-Stilwell type experiments, and by two orders of magnitude the best constraint obtained by comparing atomic clocks. This work is the first of a new generation of tests of fundamental physics using optical clocks and fiber links. As clocks improve, and as fiber links are routinely operated, we expect that the tests initiated in this paper will improve by orders of magnitude in the near future.

Special Relativity (SR), one of the cornerstones of modern physics, assumes that Lorentz Invariance (LI) is a fundamental symmetry of nature. The search for a violation of LI is motivated by two factors: (i) theoretical suggestions that LI may not be an exact symmetry at all energies and (ii) the tremendous advances in the precision of experimental tests. Indeed, a strong violation of LI at the Planck scale is likely to yield a small amount of violation at low energy, which could be measured with precise experiments [1].

Optical clocks are now the most precise measurement devices. They reach systematic uncertainties of a few 10 -18 , which can be resolved after a mere few hours of measurement with optical lattice clocks based on trapped neutral atoms [2][3][4][5][6]. Thanks to these unparalleled performances, comparing the resonance frequencies of optical clocks has led to new tests of fundamental physics, such as bounding the time variation of fundamental constants [7,8].

In this paper, we perform a test of SR using a network of distant optical lattice clocks located in France, Germany and the UK. By exploiting the difference between the velocities of each clock in the inertial geocentric frame, due to their different positions on the surface of the Earth, we are able to improve upon previous tests of time dilation. The connection between these clocks, achieved with phase-compensated optical fibers, allows * Pacome.Delva@obspm.fr for an unprecedented level of statistical resolution for the comparison of remote atomic clocks [9], making such a test possible for the first time.

LI violations are predicted by several theoretical frameworks, categorized as kinematical and dynamical frameworks (see [1] for a review). In this paper we use the Robertson-Mansouri-Sexl (RMS) kinematical framework [10][11][12][13] which contains only three parameters. It assumes the existence of a preferred frame Σ where light propagates rectilinearly and isotropically in free space with constant speed c. The ordinary Lorentz transformations from Σ to the observer frame S with relative velocity w are generalized to allow for violations of SR:

T = a -1 (t -c -1 • x) (1)
X = d -1 x -(d -1 -b -1 )( w • x) w/w 2 + wT , (2)
where a, b and d are functions of w 2 , and is a wdependent vector specifying the clock synchronization procedure in S. In the low-velocity limit:

a( w) = 1 + c -2 (α -1/2)w 2 + O(c -4 w 4 ) , (3)
where α is an arbitrary parameter quantifying the LI violation, the value of which is zero in SR.

First-order tests in w/c are based on the comparison of clocks [12]. Until recently, they gave the best constraints on the LI violating parameter α (see [14] for a review) with e.g. |α| ≤ 10 -6 obtained by comparing atomic clocks onboard GPS satellites with ground atomic clocks [15].

The three classical LI tests are the Michelson-Morley, Kennedy-Thorndike, and Ives-Stillwell experiments [10]; they are second-order tests as the LI violating signal depends on w 2 /c 2 [13]. With the advent of heavy-ion storage rings, Ives-Stillwell type experiments now give the best constraint on α [16,17]. A limit of |α| ≤ 8.4 × 10 -8 was found using 7 Li + ions prepared in a storage ring to 6.4% and 3.0% of the speed of light [16]. The experiment described in [17] uses 7 Li + ions confined at a velocity of 33.8% of the speed of light. When neglecting higher order RMS parameters, the constraint on the LI violating parameter is |α| 2.0 × 10 -8 .

In this paper, we improve upon this best previous constraint on the LI violating parameter α by a factor of around two. Our test is based on four optical lattice clocks using Sr atoms, two located at LNE-SYRTE, Observatoire de Paris, France [18,19], one at PTB, Braunschweig, Germany [20,21], and one at NPL, Teddington, UK [22]. These clocks are connected by two fiber links, one running from LNE-SYRTE to PTB operated in June 2015 [9], and one from LNE-SYRTE to NPL operated in June 2016. This paper exclusively uses the stability of the frequency comparisons between the clocks by looking for a periodic variation.

In a simplified setup, an optical clock comparison using a phase noise compensated fiber link can be described as a two-way frequency transfer between two observers A and B [23][24][25][26][27]. Observer A emits an electromagnetic signal (e.g. an IR laser) with proper frequency ν 0 , received by observer B at a proper frequency ν 1 , and reflected back to observer A, where it is received with a proper frequency ν 2 . The two-way Doppler signal (ν 2 -ν 0 ) is then locked to an ultra-stable reference -thanks to a frequency shifter inserted along the route of the fiber -such that

ν 1 ν 0 = 1 + ν 2 -ν 0 2ν 0 + ∆ , (4)
where ∆ is the "redshift signal" [14] or de-syntonization.

It contains the classical gravitational redshift between the two observer locations as well as second-order SR effects, while the first-order Doppler shift is canceled. The relative frequency difference between two clocks located at A and B is then obtained by locally beating each clock with the electromagnetic signal at each end of the link. In the low-velocity limit the de-syntonization can be written as

∆ = ∆ cl + ∆ α , (5)
where ∆ cl contains the relativistic redshift due to the static part of the gravity potential as well as temporal variations [28]. During the considered dates of clock comparaisons, time-variable gravity potential components induced by tides give approximately diurnal peak-to-peak fractional frequency variations up to 1.3 × 10 -17 between PTB and LNE-SYRTE, and up to 5 × 10 -18 between NPL and LNE-SYRTE, and are modelled at the level of two parts in 10 18 or below. The LI violating term signal is:

∆ α = αc -2 2 w • ( v A -v B ) + v 2 A -v 2 B +O(c -3 ) , (6)
where v A and v B are respectively the velocities of clocks A and B in the non-rotating geocentric celestial reference system (GCRS). They are obtained by transforming the terrestrial coordinates of the clocks, considered as constant, with the SOFA routines [START_REF]IAU SOFA Software Collection[END_REF]. w is the velocity of the Earth with respect to a preferred frame, taken as the rest frame of the cosmological microwave background (CMB). It is the sum of the Earth velocity with respect to the Solar System Barycenter (SSB), and the SSB velocity with respect to the CMB. The celestial coordinates of the SSB velocity with respect to the CMB in galactic coordinates are 263.99 • (longitude) and 48.26 • (latitude) [START_REF] Hinshaw | [END_REF], which transformed to the GCRS give 11 h 11 m 36 s (right ascension) and -6 • 54' 00" (declination) [31] with a norm of 369 km•s -1 . In June 2015 and 2016 the norm of w was w 340 km•s -1 . The first term of the LI violation in equation ( 6) varies with a period of one sidereal day as the Earth rotates around its axis. It is therefore possible to bound the LI violating parameter α by looking for daily variations in the relative frequency difference y between remote clocks, located at different longitudes (i.e. different orientation of v) and/or different latitudes (i.e. different norms of v). The second term of Eq. ( 6) is constant, and considering an upper bound of 2 × 10 -8 on the parameter α [17], is lower than 4 × 10 -20 , which is significantly below the accuracy of the clocks. Therefore we do not take it into account in our model.

We first analyze the result of the comparison between the clocks at LNE-SYRTE and NPL. Between June 10th and 15th 2016, we accumulated about 60 hours of clock comparison data between LNE-SYRTE's Sr2 and SrB lattice clocks, and NPL's Sr clock. These clocks are connected by a 812 km long cascaded optical fiber link using infra-red lasers operated at 1542 nm. The first span of 769 km connects NPL to Laboratoire de Physique des Lasers (LPL) in the north of Paris with the use of a repeater laser station at LPL [32]; the second link connects LNE-SYRTE to LPL [START_REF] Kronjaeger | 2-spans cascaded link for optical frequency transfer with 1e-18 resolution on the G ÉANT network[END_REF][START_REF] Lee | Hybrid fiber links for accurate optical frequency comparisons[END_REF]. The frequency ratio of the infra-red lasers and the Sr clock lasers at NPL and SYRTE are measured using optical frequency combs [START_REF] Johnson | [END_REF]36]. The propagation noise in the fibers is actively compensated. At LPL, a beat note is generated with light from the two stabilized links and recorded using a GPS-disciplined ultra-stable quartz oscillator and a dead-time-free frequency counter with a similar approach to the set-up described in [9]. The frequency counters at NPL, LPL, and LNE-SYRTE are synchronized to UTC(NPL), GPS time, and UTC(OP) respectively with an accuracy well below 1 ms. Figure 1 shows the relative frequency instability of the comparison.

To search for a violation of LI in the clock comparisons, we consider three different data subsets: A: Sr2 data only; B: SrB data only; C: Sr2 and SrB data com-

y NPL-SYRTE (t) = i ȳi NPL-SYRTE +2αc -2 w • [ v SYRTE (t) -v NPL (t)] , (7)
where ȳi NPL-SYRTE allows for one or two fractional frequency offsets, depending on the chosen data subset: A: i = {Sr2}; B: i = {SrB}; C: i = {Sr2, SrB}, and α is the LI violating parameter. All parameters are determined in the fitting procedure, along with correlations and uncertainties. The mean frequency offsets were removed from each comparison data subset as we are looking only for daily variations. The second line of Eq. ( 7) is the LI violation; it is very similar to a sinusoid: Q 0 sin[2π(t-t 0 )/T ] where T is one sidereal day, Q 0 = 1.60 × 10 -10 for α = 1, and t 0 = 57549.130 (MJD).

For each data subset, we used two different fitting methods: linear least-square fitting (LSQ) and affine invariant Markov Chain Monte Carlo ensemble sampler (MCMC) with 10 5 points [37]. As can be seen in Fig. 1, the Allan deviation is flat up to around 10 s averaging time. Indeed, in the short term the laser probing the narrow transition is not yet locked to the atoms, therefore the flicker floor of the free running laser is visible. This temporal correlation is taken into account in the LSQ and MCMC fits by using a non-diagonal covariance matrix. Non-diagonal elements are calculated thanks to an exponential covariance function [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], the parameters of which are fitted a priori on the empirical covariance. For the combined data set C the correlation between both data sets is taken into account in the covariance matrix.

Results for the MCMC method are given in table I for the three cases A to C. The LSQ method gives similar results, with relative differences of a few percent only. The correlations between the parameters are of the order or below 0.2. The best result is found when combining the two sets of data:

α C = (-2.83 ± 6.19) × 10 -8 (8)
Correlations and histograms of the parameters for data set C are shown on figure 2. The PTB-SYRTE comparison took place between June 4th and 24th 2015. This comparison, involving SYRTE's Sr2 clock and PTB's stationary Sr clock, is reported in [9]. We use in this paper the data of the second of the two campaigns reported in [9], representing around 150 hours of clock comparison data. An analysis of the PTB-SYRTE comparison data with a model similar to equation ( 7) (replacing the NPL clock velocity with the PTB clock velocity, and i = {Sr2}) results in a significant bias for the parameter α, five times larger than the 1σ uncertainty on α. Indeed, the power spectral density distribution of the raw data shows a significant peak at a frequency of 1 day -1 , which is around the frequency of the LI violating signal. Although this signal could be interpreted as a violation of LI, a detailed analysis shows that this effect is probably due to temperature variations in the SYRTE clock laboratory. We analyzed two independent local clock comparisons: Sr against Yb + at PTB, and Sr against Hg at SYRTE, which are not affected by a LI violation, and we used a simple model of the effect of temperature on the relative frequency difference:

y T -X (t) = α T T X (t -δt) -TX , (9)
where T X (t) is a function that interpolates the temperature at time t at some location X, TX is the mean of the temperature function T X (t -δt) evaluated for the comparison data times, α T is a temperature coefficient and δt a lag, both to be determined in the fitting procedure. A significant variation was found at 1 day -1 frequency in the local SYRTE comparison, while the comparison at PTB did not show any significant variation at this frequency.

Therefore, in addition to the LI violation, we included the effect of temperature in the SYRTE clock room, leading to the following model:

y PTB-SYRTE (t) = ȳSr2 PTB-SYRTE + y T -SYRTE (t) + 2αc -2 w • [ v SYRTE (t) -v PTB (t)] , (10)
where ȳSr2 PTB-SYRTE allows for a fractional frequency offset, y T -SYRTE is the temperature effect model given in Eq. ( 9) and α is the LI violating parameter. The relative frequency difference y PTB-SYRTE between the PTB Sr clock and the LNE-SYRTE Sr clock is corrected from the term ∆ cl . As for the model in Eq. ( 7), here the LI violating term is similar to a sinusoid with a period of one sidereal day, an amplitude Q 0 = 3.54 × 10 -10 for α = 1, and t 0 = 57177.421 (MJD). As this is a non-linear model, we use the MCMC method with 10 5 points for the fitting procedure. Note that Q 0 for the PTB-SYRTE link is more than twice the value for the NPL-SYRTE link such that it is more sensitive to a violation of LI.

The detailed result of this analysis is given in table Iline D. It shows a significant effect of the temperature on the frequency comparison of the order of 10 -16 K -1 , with a lag of around 4.8 hours. Correlations and distributions of parameters can be seen in figure 3. The lag δt SYRTE is not well constrained and its distribution is non-Gaussian. However, this does not affect the Gaussianity of the other parameters. The correlations between the parameter α and the parameters ȳSr2 PTB-SYRTE , δt SYRTE and α SYRTE T are respectively -0.48, 0.44 and 0.30. These correlations slightly degrade the uncertainty on the determination of α. The bias of the parameter |α| is below the parameter uncertainty:

α D = (-0.38 ± 1.06) × 10 -8 (11)
It is interesting to note that the bias found for the parameter α is -0.93 × 10 -8 if tides are not taken into account.

In order to check further if the choice of this temperature is significant, we repeated the same analysis with several other temperature series: (i) the SYRTE exterior (ii) the PTB exterior (iii) the PTB clock room (iv) the PTB clock laser room (v) the PTB comb room and finally (vi) a simulated sinusoid with a 1 day period. For each temperature model we fit an amplitude and a lag, as in equation (9). None of them are able to explain the residual signal, leading for the series (i-v) to biases in the parameter α ranging from 2σ up to 5.5σ, and for (vi) to a complete degeneracy of all parameters with no unique solution, which shows that the residual effect cannot be well represented by a simple sinusoidal model, as parameters are completely degenerated with parameters from the LI violating model.

The effect of temperature on the PTB-SYRTE comparison data is not yet fully understood. Further comparisons will help improve our understanding of this, thereby allowing reduction of the bias and hence the uncertainty on the determination of the parameter α. The existence of the temperature effect is however evident by the fact that it can be seen both on the distant PTB-SYRTE and the local SYRTE clock comparisons. A detailed analy-sis of the NPL-SYRTE comparisons did not show any significant systematic effects above the noise level. Consistently, a simulation of the temperature model ( 9) for the NPL-SYRTE comparisons, using the parameters determined in the PTB-SYRTE comparison, does not produce any signal above the noise level. This justifies the fact that the temperature model (9) was not used for the NPL-SYRTE comparisons.

A combination of the three data sets A, B and D has been evaluated but gives no improvement on the uncertainty of the determination of the LI violating parameter. This is due to the fact that the absolute value of the bias on α obtained in combination C (2.83 × 10 -8 ) is larger that the uncertainty obtained with data set D (1.06 × 10 -8 ).

In conclusion, by using clock comparisons between four optical clocks at NPL (UK), PTB (Germany) and SYRTE (France), linked by a leading-edge optical fiber network, we are able to put a more stringent bound on the LI violating parameter α of the RMS framework. With 1.1 × 10 -8 , α is now by around a factor of two better constrained compared to the best previous determination of this parameter, which was obtained with accelerated ions, and by two orders of magnitude with respect to the best constraint previously obtained by comparing atomic clocks. Moreover, this bound is purely limited by technical noise sources on the clock systems, which will improve in future comparisons. Projecting the comparison of distant clocks with an instability of 10 -16 / √ τ over several weeks, a reduction in uncertainty of more than one order of magnitude for α is within reach. This shows the significant potential for tests of fundamental physics with networks of optical clocks connected by optical fiber links.

FIG. 1 .

 1 FIG.1. Allan deviation of the fractional frequency difference between LNE-SYRTE's clocks Sr2 and SrB and NPL's Sr clock. After less than one day, the instability of the fractional frequency difference averages down to a few 10 -17 . This frequency instability is solely limited by the performances of the clocks, as the fiber link between SYRTE and NPL shows a fractional frequency instability of 1 × 10 -18 at 1000 s.

FIG. 2 .

 2 FIG. 2. Correlations and histograms of the parameters of model (7) in fit C of table I, corresponding to the NPL-SYRTE comparison data from Sr2 and SrB combined, using the affine invariant Markov Chain Monte Carlo ensemble sampler (MCMC) method with 10 5 points.

  81 ± 0.25 1.76 ± 0.12 -0.38 ± 1.06 TABLE I. Fitting results using the affine invariant Markov Chain Monte Carlo ensemble sampler (MCMC) method with 10 5 points. Fits A to C use the NPL-SYRTE comparison data with A: Sr2 data only; B: SrB data only; C: Sr2 and SrB data combined. Fit D uses the PTB-SYRTE comparison data.

FIG. 3 .

 3 FIG. 3. Correlations and histograms of the parameters of fit D of table I, corresponding to the fit of the PTB-SYRTE comparison data, using affine invariant Markov Chain MonteCarlo ensemble sampler (MCMC) method with 10 5 points.

ACKNOWLEDGMENTS

We thank L. De Sarlo for providing the local SYRTE Sr/Hg local comparison data, J. Achkar for providing temperature data, and P. Wolf for useful discussions. We thank S. Koke and A. Koczwara for experimental assistance on the fiber link and Harald Schnatz for longstanding support. This work would not have been possible without the support of the GIP RENATER and G ÉANT. The authors are deeply grateful to E. Camisard from GIP RENATER and G. Roberts from G ÉANT. W.-K. Lee was supported partly by the Korea Research Institute of Standards and Science under the project "Research on Time and Space Measurements", Grant No. 16011007, and also partly by the R&D Convergence Program of NST (National Research Council of Science and Technology) of Republic of Korea (Grant No. CAP-15-08-KRISS). This work was partly funded from the EC's Seventh Framework Programme (FP7 2007-2013) under Grant Agreement No. 605243 (GN3plus), and from the European Metrology Research Programme EMRP/EMPIR under SIB02 NEAT-FT, SIB55 ITOC and 15SIB05 OFTEN. This project has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme. LNE-SYRTE and LPL acknowledge funding support from the Agence Nationale de la Recherche (ANR blanc LIOM 2011-BS04-009-01, Labex First-TF ANR 10 LABX 48 01, Equipex REFIMEVE+ ANR-11-EQPX-0039), Action spécifique GRAM, Centre National d' Études Spatiales (CNES) and Conseil Régional Île-de-France (DIM Nano-K). PTB acknowledges support from the German Research Foundation DFG within RTG 1729, CRC 1128 geo-Q, and CRC 1227 DQ-mat. W.B. and S.B. would like to acknowledge the support from the Marie Curie Initial Training Network FACT.