Alexandre Bône

Maxime Louis

Alexandre Routier

Jorge Samper

Michael Bacci

Benjamin Charlier

Olivier Colliot

Stanley Durrleman

Prediction of the progression of subcortical brain structures in Alzheimer's disease from baseline

We propose a method to predict the subject-specific longitudinal progression of brain structures extracted from baseline MRI, and evaluate its performance on Alzheimer's disease data. The disease progression is modeled as a trajectory on a group of diffeomorphisms in the context of large deformation diffeomorphic metric mapping (LDDMM). We first exhibit the limited predictive abilities of geodesic regression extrapolation on this group. Building on the recent concept of parallel curves in shape manifolds, we then introduce a second predictive protocol which personalizes previously learned trajectories to new subjects, and investigate the relative performances of two parallel shifting paradigms. This design only requires the baseline imaging data. Finally, coefficients encoding the disease dynamics are obtained from longitudinal cognitive measurements for each subject, and exploited to refine our methodology which is demonstrated to successfully predict the follow-up visits.

Introduction

The primary pathological developments of a neurodegenerative disease such as Alzheimer's are believed to spring long before the first symptoms of cognitive decline. Subtle gradual structural alterations of the brain arise and develop along the disease course, in particular in the hippocampi regions, whose volumes are classical biomarkers in clinical trials. Among other factors, those transformations ultimately result in the decline of cognitive functions, which can be assessed through standardized tests. Being able to track and predict future structural changes in the brain is therefore key to estimate the individual stage of disease progression, to select patients and provide endpoints in clinical trials.

To this end, our work settles down to predict the future shape of brain structures segmented from MRIs. We propose a methodology based on three building blocks : extrapolate from the past of a subject ; transfer the progression of a reference subject observed over a longer time period to new subjects ;

and refine this transfer with information about the relative disease dynamics extracted from cognitive evaluations. Instead of limiting ourselves to specific features such as volumes, we propose to see each observation of a patient at a given time-point as an image or a segmented surface mesh in a shape space.

In computational anatomy, shape spaces are usually defined via the action of a group of diffeomorphisms [START_REF] Beg | Computing large deformation metric mappings via geodesic flows of diffeomorphisms[END_REF][START_REF] Vercauteren | Large deformation diffeomorphism and momentum based hippocampal shape discrimination in dementia of the Alzheimer type[END_REF]17]. In this framework, one may estimate a flow of diffeomorphisms such that a shape continuously deformed by this flow best fits repeated observations of the same subject over time, thus leading to a subjectspecific spatiotemporal trajectory of shape changes [START_REF] Lorenzi | 4D registration of serial brains MR images: a robust measure of changes applied to Alzheimer's disease[END_REF][START_REF] Peyrat | Registration of 4D time-series of cardiac images with multichannel diffeomorphic Demons[END_REF]. If the flow is geodesic in the sense of a shortest path in the group of diffeomorphisms, this problem is called geodesic regression [START_REF] Fishbaugh | Geodesic regression of image and shape data for improved modeling of 4D trajectories 5[END_REF]5,[START_REF] Lorenzi | 4D registration of serial brains MR images: a robust measure of changes applied to Alzheimer's disease[END_REF][START_REF] Peyrat | Registration of 4D time-series of cardiac images with multichannel diffeomorphic Demons[END_REF] and may be thought of as the extension to Riemannian manifolds of the linear regression concept. It is tempting then to use such regression to infer the future evolution of the shape given several past observations. To the best of our knowledge, the predictive power of such a method has not yet been extensively assessed. We will demonstrate that satisfying results can only be obtained when large numbers of data points over extensive periods of time are available, and that poor ones should be expected in the more interesting use-case scenario of a couple of observations. In such situations, an appealing workaround would be to transfer previously acquired knowledge from another patient observed over a longer period of time. This idea requires the definition of a spatiotemporal matching method to transport the trajectory of shape changes into a different subject space. Several techniques have been proposed to register image time series of different subjects [START_REF] Metz | Nonrigid registration of dynamic medical imaging data using nd + t b-splines and a groupwise optimization approach[END_REF][START_REF] Wu | Estimating the 4D respiratory lung motion by spatiotemporal registration and building super-resolution image[END_REF]. They often require time series to have the same number of images, or to have correspondences between images across time series, and are therefore unfit for prognosis purposes. Parallel transport in groups of diffeomorphisms has been recently introduced to infer deformation of follow-up images from baseline matching [START_REF] Lorenzi | Geodesics, parallel transport & one-parameter subgroups for diffeomorphic image registration[END_REF][START_REF] Singh | Hierarchical geodesic models in diffeomorphisms[END_REF]. Such paradigms have been used mostly to transport spatiotemporal trajectories to the same anatomical space for hypothesis testing [START_REF] Hadj-Hamou | Longitudinal analysis of image time series with diffeomorphic deformations: A computational framework based on stationary velocity fields[END_REF][START_REF] Qiu | Parallel transport in diffeomorphisms distinguishes the time-dependent pattern of hippocampal surface deformation due to healthy aging and the dementia of the Alzheimer's type[END_REF]. Two main methodologies have emerged: either by parallel-transporting the time series along the baseline matching as in [5], or by parallel-transporting the baseline matching along the time series as in [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF]. We evaluate both in this paper.

In any case, these approaches require to match the baseline shape with one in the reference time series. Ideally, we should match observations corresponding to the same disease stage, which is unknown. We propose to complement such approaches with estimates of the patient stage and pace of progression using repeated neuropsychological assessments in the spirit of [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF]. These estimates are used to adjust the dynamics of shape changes of the reference subject to the test one, according to the dynamical differences observed in the cognitive tests.

Among the main contributions of this papers are : the first quantitative study of the predictive power of geodesic regression ; a new methodology for the prediction of shape progression from baseline ; the evaluation of its accuracy for two different parallel shifting protocols ; new evidence of the utter importance of capturing the individual dynamics in Alzheimer's disease models.

Section 2 sets the theoretical background and incrementally describes our methodology. Section 3 presents and discusses the resulting performances.

Method

Let (y j ) j=1,..,ni be a time series of segmented surface meshes for a given subject i ∈ {1, ..., N }, obtained at the ages (t j ) j=1,..,ni . We build a group of diffeomorphisms of the ambient space which act on the segmented meshes, following the procedure described in [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. Flows of diffeomorphisms of R 3 are generated by integrating time-varying vector fields of the form v(t, x)

= ncp k=1 K[x, c k (t)]β k (t)
where K is a Gaussian kernel, c(t) = [c k (t)] k=1,..,ncp and β(t) = [β k (t)] k=1,..,ncp are respectively the control points and the momenta of the deformation.

We endow the space of diffeomorphisms with a norm which measures the cost of the deformation. In the following, we only consider geodesic flows of diffeomorphisms i.e. flows of minimal norm connecting the identity to a given diffeomorphism. Such flows are uniquely parametrized by their initial control points and momenta c 0 = c(0), β 0 = β(0). Under the action of the flow of diffeomorphisms, an initial template shape T is continuously deformed and describes a trajectory in the shape space, which we will note t → γ (c 0 ,β 0 ) (T, t). Simultaneously, we endow the surface meshes with a varifold norm • which allows to measure a data attachment term between meshes without point correspondence [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF].

Geodesic regression

In the spirit of linear regression, one can perform geodesic regression in the shape space by estimating the intercept T and the slope (c 0 , β 0 ) such that γ (c 0 ,β 0 ) (T, •) minimizes the following functional :

inf c 0 ,β 0 ,T ni j=1 γ (c 0 ,β 0 ) (T, t j ) -y j 2 + R(c 0 , β 0 ) ( 1
)
where R is a regularization term which penalizes the kinetic energy of the deformation. We estimate a solution of equation ( 1) with a Nesterov gradient descent as implemented in the software Deformetrica (www.deformetrica.org), where the gradient with respect to the control points, the momenta and the template is computed with a backward integration of the data attachement term along the geodesic [START_REF] Durrleman | Sparse adaptive parameterization of variability in image ensembles[END_REF]. Once an optimum is found, we obtain a description of the progression of the brain structures which lies in the tangent space at the identity of the group of diffeomorphisms. It is natural to attempt to extrapolate from the obtained geodesic to obtain a prediction of the progression of the structures.

Two methods to transport spatiotemporal trajectories of shapes

As it will be demonstrated in section 3, geodesic regression extrapolation produces an accurate prediction only if data over a long time span is available for the subject, which is not compatible with the goal of early prognosis.

As proposed in [START_REF] Lorenzi | Geodesics, parallel transport & one-parameter subgroups for diffeomorphic image registration[END_REF][START_REF] Younes | Jacobi fields in groups of diffeomorphisms and applications[END_REF], given a reference geodesic, we use the Riemannian parallel transport to generate a new trajectory. We first perform a baseline matching between the reference subject and the new subject, which can be described as a vector in the tangent space of the group of diffeomorphisms. Two paradigms are available to obtain a parallel trajectory. [START_REF] Singh | Hierarchical geodesic models in diffeomorphisms[END_REF] advises to transport the reference regression along the matching and then shoot. In the shape space, this generates a geodesic starting at the baseline shape ; for this reason, we call this solution geodesic parallelization, and is illustrated on Figure (A1). On the other hand, [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF] advocates to transport the matching vector along the reference geodesic and then build a trajectory with this transported vector from every point of the reference geodesic, as described on Figure (B1). We will call this procedure exp-parallelization.

In such a high-dimensional setting, the computation of parallel transport classically relies on the Schild's ladder scheme [START_REF] Lorenzi | Schild's ladder for the parallel transport of deformations in time series of images[END_REF]. However, in our case the computation of the Riemannian logarithm may only be computed by solving a shape matching problem, resulting not only in an computationally expensive algorithm but also in an uncontrolled approximation of the scheme. To implement these parallel shifting methods, we use the algorithm suggested in [START_REF] Younes | Jacobi fields in groups of diffeomorphisms and applications[END_REF], which relies on an approximation of the transport to nearby points by a well-chosen Jacobi field, with a sharp control on the computational complexity. The same rate of convergence as Schild's ladder is obtained at a reduced cost.

Cognitive scores dynamics

The protocol described in the previous section has two main drawbacks. First, the choice of the matching time in the reference trajectory is arbitrary : the baseline is purely a convenience choice and ideally the matching should be performed at similar stages of the disease. Second, it does not take into account the pace of progression of the subject. In [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF], the authors propose a statistical model allowing to learn, in an unsupervised manner, dynamical parameters of the subjects from ADAS-cog test results, a standardized cognitive test designed for disease progression tracking. More specifically, they suppose that each patient follows a parallel to a mean trajectory, with a time reparametrization :

ψ(t) = α(t -t 0 -τ ) + t 0 (2)
which maps the subject time to a normalized time frame, where α > 0 and τ are scalar parameters. A high (resp. low) α hence corresponds to a fast (resp. slow) progression of the scores, when a negative (resp. positive) τ corresponds to an early decay (resp. late decay) of those scores. In the dataset introduced below, the acceleration factors (α i ) i range from 0.15 to 6.01 and the time-shifts (τ i ) i from -20.6 to 22.8, thus showing a tremendous variability in the individual dynamics of the disease, which must be taken into account. With these dynamic parameters, the shape evolution can be adjusted by reparametrizing the parallel trajectory with the same formula (2), as illustrated on Figures (A2) and (B2).

Results

Data, preprocessing, parameters and performance metric

MRIs are extracted from the ADNI database, where only MCI converters with 7 visits or more are kept, for a total of N = 74 subjects and 634 visits. Subjects are observed for a period of time ranging from 4 to 9 years (5.9 on average), with 12 visits at most. The 634 MRIs are segmented using the FreeSurfer software. The extracted brain masks are then affinely registered towards the Colin 27 Average Brain using the FSL software. The estimated transformations are finally applied to the pairs of caudates, hippocampi and putamina subcortical structures.

All diffeomorphic operations i.e. matching, geodesic regression estimation, shooting, exp-parallelization and geodesic parallelization are performed thanks to the Deformetrica software previously mentioned. A varifold distance with Gaussian kernel width of 3 mm for each structure and a deformation kernel width of 5 mm are chosen. The time discretization resolution is set to 2 months.

The chosen performance metric between two sets of meshes is the Dice coefficient, that is the sum of the volumes of the intersections of the corresponding meshes, divided by the total sum of the volumes. We only measure the volume of the intersection between corresponding structures. The Dice coefficient is comprised between 0 and 1 : it equals 1 for a perfect match, and 0 for disjoint structures.

Geodesic regression extrapolation

The acceleration factor α in equation ( 2) encodes the rate of progression of each patient. Multiplying this coefficient with the actual observation window gives a notion of the absolute observation window length, in the disease time referential. Only the 22 first subjects according to this measure have been considered for this section : they are indeed expected to feature large structural alterations, making the geodesic regression procedure more accurate. The geodesic regression predictive performance is compared to a naive one consisting in leaving the last observed brain structures in the learning dataset unchanged.

Table 1 presents the results obtained for varying learning dataset and extrapolation extents. We perform a Mann-Whitney test with the null hypothesis that the observed Dice coefficients distributions are the same to obtain the statistical significance levels. The extrapolated meshes are satisfying only in the case where all but one data points are used to perform the geodesic regression, achieving a high Dice index and outperforming the naive one, by a small margin though and failing to reach the significance level (p = 0.25). When the window of observation becomes narrower, the prediction accuracy decreases and becomes worse than the naive one. Indeed, the lack of robustness of the -although standardsegmentation pipeline imposes a high noise level, which seems to translate into a too low signal-to-noise ratio after extrapolation from only a few observations. Table 1: Averaged Dice performance measures between predictions and observations for varying extents of learning datasets and extrapolation. The [reg] tag indicates the regression-based prediction, and [naive] the naive one. Each row corresponds to an increasingly large learning dataset, patients being observed for widening periods of time. Each column corresponds to an increasingly remote predicted visit from baseline. Significance levels [.05, .01, .001, .0001] for the Mann-Whitney test. Figure 1 displays an extrapolated geodesic regression for a specific subject, with a large learning period of 72 months, and a prediction at 108 months from the baseline (Dice performance of 0.74 versus 0.65 with the naive approach).

Non reparametrized transport

Among the 22 subjects whose regression-based predictive power has been evaluated in the previous section, the two which performed best are chosen as references for the rest of this paper. Their progressions are transported onto the 73 other subjects with the two different parallel shifting methods.

In more details, for each pair of reference and target subjects, the baseline target shape is first registered to the reference baseline. The reference geodesic regression is then either geodesically or exp-parallelized. Prediction performance is finally assessed : the Dice index between the prediction and the actual observation, for the two modes of transport, are computed and compared to the Dice index between the baseline meshes and the actual observation -the only available information in the absence of a predictive paradigm.

The upper part of Table 2 presents the results. In most cases, the obtained meshes by the proposed protocol are of lesser quality than the reference ones, according to the Dice performance metric. The two methods of transport are essentially similarly predictive, although geodesic parallelization slightly outperforms the exp-parallelization for the M12 prediction.

Refining with cognitive dynamical parameters

The two reference progressions are transported through geodesic and exp-parallelization onto all remaining subjects. After time-reparametrization, the obtained parallel trajectories then deliver predictions for the brain structures. Figure 3 displays a reference geodesic and an exp-parallelized curve. The predicted progression graphically matches the datapoints, and it can be noticed that the final prediction at age 85.6 (Dice 0.73) outperforms the corresponding one on Figure 2, obtained without time-reparametrization (Dice 0.69).

Quantitative results are presented in the lower part of Table 2. At the exception of the M12 prediction, both protocols outperform the naive one. The M36, M48, M72 and M96 predictions are the most impressive ones, with p-values always lesser than 1%. This shows that the pace of cognitive score evolution is well correlated with the pace of structural brain changes, and therefore allows an enhanced prediction of follow-up shapes.

No conclusion can be drawn concerning the two parallel shifting methodologies, a single weak significance result being obtained only for the M12 prediction where the geodesic parallelization method slightly outperforms the expparallelization one with a Dice score of 0.888 versus 0.882.

Conclusion

We conducted a quantitative study of geodesic regression extrapolation, exhibiting its limited predictive abilities. We then proposed a method to transport a spatiotemporal trajectory into a different subject space with cognitive declinederived time reparametrization, and demonstrated its potential for prognosis. The results show how crucial the dynamics are in disease modeling, and how cross-modality data can be exploited to improve a learning algorithm. The two main paradigms that have emerged for the transport of parallel trajectories were shown to perform equally well in this prediction task. Nonetheless, the exp-parallelization offers a methodological advantage in that the generated trajectories do not depend on a particular choice of point on the reference geodesic, in contrast with the trajectories obtained by geodesic parallelization. It takes full advantage of the isometric property of the parallel transport, and eases the combination with time-warp functions based on the individual disease dynamics.

In future work, more complex time reparametrization could be considered as in [START_REF] Hong | Time-warped geodesic regression[END_REF]. Finally, the robustness of the proposed protocol to the choice of reference subject has not been assessed. Such a choice could be avoided by constructing an average disease model as in [START_REF] Singh | Hierarchical geodesic models in diffeomorphisms[END_REF], or by translating for shapes the method of [START_REF] Schiratti | Learning spatiotemporal trajectories from manifold-valued longitudinal data[END_REF]. We may also use this framework to estimate a joint image and cognitive model to better estimate individual dynamical parameters of disease progression.

(

  A1) Geodesic parallelization. Blue arrow: baseline matching. Red arrows: transported regression. Black dotted line : exponentiation of the transported regression. (A2) Reparametrized geodesic parallelization. Matching time and expparallel trajectory are reparametrized. (B1) Exp-parallelization. Red arrow: geodesic regression. Blue arrows: transported baseline matching. Black dotted line : exp-parallelization of the reference geodesic for the given subject. (B2) Reparametrized expparallelization. Matching time and expparallel trajectory are reparametrized.

  months) for all subjects (N = 22).

Fig. 1 :

 1 Fig. 1: Extrapolated geodesic regression for the subject s0671. Are only represented the right hippocampus, caudate and putamen brain structures in each subfigure. The three first rows present the interpolated brain structures, corresponding to ages 61.2, 64.2 and 67.2 (years). The last row presents the extrapolation result at age 70.2. On the right column are added the target brain structures (red wireframes), segmented from the original images.

Fig. 2 :

 2 Fig. 2: Exp-parallelization of the reference subject s0906 (first column) towards the subject s1080 (second column), giving predictions for ages 81.6, 82.6, 83.6, 84.6 and 85.6 (years). On the third column are added the target brain structures (red wireframes), segmented from the original images.

Fig. 3 :

 3 Fig. 3: Time-reparametrized exp-parallelization of the reference subject s0906 (first column) towards the subject s1080 (second column), giving predictions for ages 81.6, 82.6, 83.6, 84.6 and 85.6 (years). On the third column are added the target brain structures (red wireframes), segmented from the original images.

Table 2 :

 2 Averaged Dice performance measures between predictions and obervations for two modes of transport, with or without refinement by the cognitive scores. In each cell, the first line corresponds to the exp-parallelization-based prediction [exp], the middle line to the geodesic parallelization-based one [geod], and the last line to the naive approach [naive]. Each column corresponds to an increasingly remote predicted visit from baseline. Significance levels for the Mann-Whitney test [.05, .01, .001, .0001].

† Equal contributions.

Acknowledgments. This work has been partly funded by the European Research Council (ERC) under grant agreement No 678304, European Union's Horizon 2020 research and innovation program under grant agreement No 666992, and the program Investissements d'avenir ANR-10-IAIHU-06.


Guilhem Saurel

Justin Carpentier

Nicolas Mansard

Jean-Paul Laumond

A Simulation Framework for Simultaneous Design and Control of Passivity Based Walkers

In this paper, we propose a simulation framework which simultaneously computes both the design and the control of bipedal walkers. The problem of computing a design and a control is formulated as a single large-scale parametric optimal control problem on hybrid dynamics with path constraints (e.g. non sliding and non slipping contact constraints). Our framework relies on state-of-the-art numerical optimal control techniques and efficient computation of the multi-body rigid dynamics. It allows to compute both the parametrized model and the control of passive walkers on different scenarios, in only few seconds on a standard computer. The framework is illustrated by several examples which highlight the interest of the approach.

I. INTRODUCTION

Passive walkers are bipedal robots essentially powered by gravity. They exploit their natural dynamics to move forward, but in the meanwhile they are unable to exhibit quasi-static behaviors. Such mechanical systems show an excellent cost of transport (CoT). They are effective platforms, both for biomechanicians to better understand the essence of bipedal walk, and for robot designers to build efficient humanoid robots.

The purpose of this paper is to propose a generic dynamic simulation framework for optimizing both the design and the controller of such robots with respect to a given cost function. An overview of the framework is given in Fig. 1. [START_REF] Mcgeer | Passive dynamic walking[END_REF], many passive walkers have been designed and crafted. A great introduction to this field is given in Collins et al. [START_REF] Collins | Efficient bipedal robots based on passive-dynamic walkers[END_REF]. A complete methodology to build incrementally complex passive walkers is described by Wisse et al. [START_REF] Wisse | Delft pneumatic bipeds[END_REF]. It starts from a simple compass-like model and goes to Denise, a 3D dynamic walking robot with a bisecting mechanism for the torso, two arms rigidly coupled to the hip angle, knees that are unlocked during the swing phase, and ankles engineered to steer in the direction that prevents falling. Such incremental approach reveals the critical issues in the design of walker with an anthropomorphic shape. The first issue is the transition from the simplest compass model evolving in the sagittal plane to 3D system [START_REF] Collins | A three-dimensional passive-dynamic walking robot with two legs and knees[END_REF], [START_REF] Tedrake | Actuating a simple 3d passive dynamic walker[END_REF], [START_REF] Hobbelen | System overview of bipedal robots flame and tulip: Tailor-made for limit cycle walking[END_REF]. A second issue is the extension of the compass model to articulated legs [START_REF] Mcgeer | Passive dynamic walking[END_REF], [START_REF] Collins | A three-dimensional passive-dynamic walking robot with two legs and knees[END_REF], [START_REF] Hobbelen | System overview of bipedal robots flame and tulip: Tailor-made for limit cycle walking[END_REF], [START_REF] Ikemata | A physical principle of gait generation and its stabilization derived from mechanism of fixed point[END_REF], [START_REF] Westervelt | Feedback control of dynamic bipedal robot locomotion[END_REF]. The addition of ankles is addressed in [START_REF] Hobbelen | System overview of bipedal robots flame and tulip: Tailor-made for limit cycle walking[END_REF], [START_REF] Wisse | Ankle springs instead of arc-shaped feet for passive dynamic walkers[END_REF], arms in [START_REF] Collins | A three-dimensional passive-dynamic walking robot with two legs and knees[END_REF], [START_REF] Hobbelen | System overview of bipedal robots flame and tulip: Tailor-made for limit cycle walking[END_REF] and even neck and head in [START_REF] Benallegue | Contribution of actuated head and trunk to passive walkers stabilization[END_REF], [START_REF] Falotico | Head stabilization in a humanoid robot: models and implementations[END_REF]. The mass distribution {saurel,carpentier,mansard,jpl}@laas.fr for a given kinematic structure is also an important question which is addressed in [START_REF] Hass | Optimal mass distribution for passivity-based bipedal robots[END_REF], [START_REF] Remy | A matlab framework for efficient gait creation[END_REF].

A. Related Work

Starting with McGeer

Beyond those mechatronical contributions, the research also focuses on the minimal energy cosumption [START_REF] Collins | Efficient bipedal robots based on passive-dynamic walkers[END_REF], [START_REF] Bhounsule | Low-bandwidth reflex-based control for lower power walking: 65 km on a single battery charge[END_REF] and on the stability properties of passive walkers [START_REF] Ikemata | A physical principle of gait generation and its stabilization derived from mechanism of fixed point[END_REF], [START_REF] Benallegue | Contribution of actuated head and trunk to passive walkers stabilization[END_REF], [START_REF] Mombaur | Stability optimization of open-loop controlled walking robots[END_REF], [START_REF] Grizzle | Asymptotically stable walking for biped robots: Analysis via systems with impulse effects[END_REF], [START_REF] Byl | Metastable walking machines[END_REF].

By opposition to non-passive preview control approaches (e.g. Kajita [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF], capture-point [START_REF] Pratt | Capture point: A step toward humanoid push recovery[END_REF]) that require footstep planning, the control of passive walkers aims at maintaining at the best the natural steady gait that originates from the mechanical design. Passive walkers do not compute their footsteps in advance. Mechanical design and control are deeply connected and the challenge is to consider both simultaneously. The natural dynamics of passive walkers is periodic. The mechanical design tends to create limit cycles at the origin of the steady gait. The role of the controller is to maintain the system around this limit cycle in spite of environment perturbation.

In model-based optimization, the problem of finding a good controller is expressed as a constrained numerical optimization problem that benefits from solid theoretical analysis of stability [START_REF] Diehl | Stability optimization of hybrid periodic systems via a smooth criterion[END_REF]. This general numerical framework has been applied to passive bipedal walkers for the optimization of mass distribution [START_REF] Hass | Optimal mass distribution for passivity-based bipedal robots[END_REF]. It has been also instantiated to the study of passive quadrupeds [START_REF] Remy | Stability analysis of passive dynamic walking of quadrupeds[END_REF] and then to passive walkers for the creation and analysis of efficient gaits [START_REF] Remy | A matlab framework for efficient gait creation[END_REF]. In this last paper, the authors propose a dynamic simulation package whose scope is illustrated by various examples including passive walkers and running robots. Our paper falls in the same framework.

B. Problem statement

In this paper, we want to simultaneously solve the two following objectives:

1) for a given kinematic chain architecture, we want to optimize the parameters of a robot (e.g. mass distribution, body segments lengths, slope of the ground, mean forward velocity, etc.), 2) for a given robot, we want to find the best controller with respect to a cost function (e.g. the cost of transport, the minimal time, etc.). In short, for a given kinematic structure, we propose a generic approach to design all robot parameters together with the controller that optimizes a given cost function.

C. Contribution

The first originality of this approach is its ability to deal with complex architectures like human-like walkers, both in 2D and 3D. We do not restrict the motion to only lie in the sagittal plane. Moreover and contrary to similar works on this topic, our framework automatically computes the full dynamics. Therefore, there is no need of writing down the complex dynamic equations of polyarticulated systems. Hence, many passive walkers can be efficiently designed, optimized and compared.

Secondly, control can be chosen to be active or passive. We also handle periodic as well as non-periodic gaits.

Finally, we can optimize various parameters of a given walker (slope, lengths, masses, speeds, etc.) with respect to a given cost function. Fig. 1 illustrates the global architecture of our framework.

D. Outline of the paper

In Sec. II, we first introduce the dynamic contact simulator at the root of the framework. We highlight how the various problem parameters are distributed either in the mechanical model or in the controller. In Sec. III, we set up the generic optimal control formulation which allows the computation of both design and control parameters. Finally, the setup is introduced in Sec. IV and experimental results are presented in Sec. V.

II. DYNAMIC CONTACT SIMULATION

Passive walkers are intrinsically hybrid systems. They are submitted to a continuous dynamics when the stance leg is in contact with the ground, and they are also subject to impacts when the swing foot hits the ground. In addition to this hybrid dynamics, some contact constraints must be satisfied to ensure the feasibility of the entire motion.

In the following of this section, we expose the notations, the parametric model of the walkers and the contact formulation used inside the framework.

A. Notations

In this paper, we assimilate a passive walker to a free-floating base system. We denote by q ∈ SE(3) × R n its configuration vector, with SE(3) the special Euclidian group of dimension 3 encoding the placement of the robot base and n the number of degrees of freedom (DoF). The tangent velocity and acceleration of the configuration vector are denoted by q and q respectively and live in R 6+n . Finally, the torque applied at each joint is denoted by τ ∈ R n .

B. The parametric model

A passive walker is primarily a kinematic tree, i.e. a tree of joints where each joint has its own topology (e.g. revolute, free-flyer, spherical, etc.) and a particular placement regarding to its parent joint. The joints are the nodes of the tree. In addition, each joint supports a body, which is defined by its mass, its center of mass (CoM) and its inertia matrix. All the bodies together define the mass distribution of the model. The tree structure with the mass distribution correspond to the structural parameters of the system. The model of the passive walker is then parametrized by those two sets of parameters: model(tree, mass_distribution)

C. The parametric controller

A gait is characterized by its controller which is represented by a set of real parameters. For instance, a controller may be a set of splines which encode the torque trajectories or just the PID gains (stiffness and damping of the spring attached to the joint) in the case of pure passive controller.

controller(control_parameters)

D. Continuous contact dynamics

For the continuous dynamics, we make the hypothesis of punctual rigid contact with Coulomb friction cone. The dynamic equation of the polyarticulated system with constraints can be stated as:

M (q)q + b(q, q) = S τ + J c (q) f c (1)
J c (q)q + Jc (q, q) q = 0

where M (q) is the joint space inertia matrix, b(q, q) corresponds to the Coriolis, centrifugal and gravitational effects, S is a selection matrix encoding the underactutation, J c (q) is the contact Jacobian with f c the contact forces and . denotes the transpose operator.

A necessary and sufficient condition for non-sliding and non-slipping of the contact is described by the constraint that f c must remain inside the Coulomb friction cone K c . This cone reflects the fact that the normal component of the contact force is positive (the ground cannot pull), the norm of its tangential components and the normal torque are limited by the normal component.

We make the choice to solve (1)-( 2) together, and to add the conic constraint inside the optimal control problem to enforce the Coulomb contact model [START_REF] Brogliato | Nonsmooth mechanics[END_REF]. This implies that we have to impose contact phases where (1)-( 2) are enforced. Other approaches have tried to get rid of this extra hypothesis [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with coulomb friction[END_REF], [START_REF] Brogliato | Nonsmooth mechanics[END_REF] leading to difficult and still incompletely solved problems for trajectory optimisation [START_REF] Tassa | Synthesis and stabilization of complex behaviors through online trajectory optimization[END_REF], [START_REF] Posa | A direct method for trajectory optimization of rigid bodies through contact[END_REF]. In the context of this paper, selecting in advance the contact phases is not a limitation.

The joint acceleration and contact force vectors are then given by:

q = M -1 (I n -J t c Λ -1 c J c M -1 )(S τ -b) (3) -M -1 J c Λ -1 c Jc q f c = Λ -1 c J c M -1 (b -S t τ ) -Jc q (4)
with Λ c def = J c M -1 J c the so-called operational-space inertia matrix and I n the identity matrix of dimension n. The dependences on q and q have been omitted for simplicity of notation.

E. Impact dynamics

Passive walkers are also subject to impacts. Here, we make the assumption of instantaneous inelastic impacts with a restitution coefficient sets to zero, i.e. the post-impact velocity of the contact point is null. The impact dynamics then leads to a discontinuity in the joint velocity space, which is described by the two following equations:

q+ = (I n -M -1 J c Λ -1 c J c ) q- (5) λ c = Λ -1 c J c q- (6)
with q-, q+ the pre-impact and post-impact generalized velocities and λ c is the impulse resulting from the impact [START_REF] Brogliato | Nonsmooth mechanics[END_REF]. Other impact model (e.g.elastic) could also be introduced without loss of generality. This impact model is frequently used in the litterature [START_REF] Schultz | Modeling and optimal control of human-like running[END_REF], even though it is contested for its physical consistency [START_REF] Chatterjee | A new algebraic rigid-body collision law based on impulse space considerations[END_REF].

F. Efficient rigid-body dynamic computation

The optimal control solver must evaluate thousands of times the multi-body dynamics either inside the numerical integration procedure or for the evaluation of the dynamic sensitivities regarding to the model and control parameters. For that purpose, we used Pinocchio [START_REF] Carpentier | Pinocchio: fast forward and inverse dynamics for poly-articulated systems[END_REF], a whole new, open-source and efficient C++ library to model and compute the forward and the inverse dynamics of polyarticulated systems in contact. Pinocchio is written on top of the efficient Eigen C++ library [START_REF] Guennebaud | Eigen v3[END_REF] for Linear Algebra. It is based on Featherstone's algorithms [START_REF] Featherstone | Rigid body dynamics algorithms[END_REF] but they have been implemented in a way to take profit of branch prediction and cache mechanisms of modern processors [START_REF] Naveau | METAPOD Template META-programming applied to dynamics: CoP-CoM trajectories filtering[END_REF].

III. OPTIMAL CONTROL APPROACH FOR DESIGN AND

CONTROL

Optimal control is a powerful and generic mathematical tool which allows the exhibition of a particular solution among an infinite number of candidates. While geometric results only exist for a very limited class of systems, the past few years have seen the emergence of efficient and reliable numerical optimal control frameworks working on high dimensional and complex systems [START_REF] Leineweber | An efficient multiple shooting based reduced sqp strategy for large-scale dynamic process optimization. part 1: theoretical aspects[END_REF], [START_REF] Houska | ACADO Toolkit -An Open Source Framework for Automatic Control and Dynamic Optimization[END_REF].

In our framework, the simultaneous search of model and control parameters is set up as an optimal control problem with a prescribed cost function. This cost function reflects the objective of the gait and can be any real value function. Some examples of cost functions in the context of passive walkers are the cost of transport or even the minimal time. In addition, we add the possibility to set the duration of the motion as a free parameter of the problem. Other free variables, like the gait stride length for example or the slope of the ground, can be stacked to the list of parameters.

A. Notations

We denote by x def = (q, q) the state of the systems, u is the control vector and p is the vector of parameters, composed of both the model parameters and the aforementioned free variables.State and control trajectories are denoted by x and u respectively. The cost function and dynamics of the system are then written as (t, x, u, p) and dx dt = f (t, x, u, p) respectively. Here, we use a slight abuse of notations to denote with the same notation both the continuous and the impact dynamics. Finally, g(t, x, u, p) corresponds to the inequality constraints that must be satisfied along the path (equality constraint can be also considered without loss of generality). As mentioned in Sec. II-D, the function g is first and foremost composed of the Coulomb conic constraints.

B. The Optimal Control Problem formulation

The hybrid dynamics of passive walkers can be seen as a multi-phase system, each stage corresponding either to the single, double support or impact phases. Thereafter, the integer s refers to the index of the s th stage.

The generic optimal control problem for simultaneously computing model and control parameters with multi-phase dynamics can be written as:

min x,
where π in (7d) is a function which acts both on the state and control trajectories to enforce periodicity constraints. ∆t s is the duration of the phase s, and T = ∆t s is the total duration of the motion. In case of impacts, the phase duration is reduced to 0.

C. Solving the Optimal Control Problem

Two major directions exist to solve the infinite dimensional problem [START_REF] Ikemata | A physical principle of gait generation and its stabilization derived from mechanism of fixed point[END_REF]. The first direction belongs to the so-called indirect methods. It consists in exploiting the necessary conditions for optimality, namely the Pontryagin's maximum principle [START_REF] Liberzon | Calculus of variations and optimal control theory: a concise introduction[END_REF], which transforms the problem (7) into a boundary value problem working on ordinary differential equations. However, such methods are currently unable to track path constraints.

The second direction corresponds to direct methods. Direct methods first discretize the original problem into a finite dimensional nonlinear programming problem (NLP), which is then solved with standard NLP strategies. Among NLP strategies, three of them are now popular: (i) single shooting, (ii) collocation and (iii) multiple shooting. In what follows, we briefly survey these three methods. For further details, we refer the reader to the general overview on numerical optimal control methods written by Diehl et al. [START_REF] Diehl | Fast Direct Multiple Shooting Algorithms for Optimal Robot Control[END_REF].

1) Single shooting: discretizes the control and constraints according to a temporal grid. The state trajectory is recovered by integration of the discrete control trajectory along this grid. As single shooting method reduces the NLP to the search of a control trajectory, the optimization problem is of low dimension. However, the solver is hard to initialize if only an initial guess on the state trajectory is available or it may not converge at all in the context of unstable systems.

2) Collocation: disctrizes both the control and the state trajectories according to a temporal grid. In addition to the classic discritized constraints, the state trajectory is enforce to match the dynamics equation (7b) at each grid node. The problem can then be easily initialized from a given state trajectory and collocation handles well unstable dynamics. However, a very fine grid is required to make the state trajectory closer to the true dynamics of the system.

3) Multiple shooting: takes profit of both previous methods. It works on a coarser time grid, which are called multiple shooting intervals. On each interval, the control is discretized as well as the initial state value. The final state value on the interval is then obtained by integration of the system dynamics (7b). In this way, each interval is set independent from its neighbours. And the dependencies between successive intervals is shifted as equality constraints of the NLP. The NLP remains a low dimensional problem and can be easily warm-started with an initial guess on the state trajectory. Furthermore, multiple shooting is really suited for multi-phase dynamics, as each phase is set independent to the others.

Multiple shooting has been successfully applied for the modelling of human running [START_REF] Schultz | Modeling and optimal control of human-like running[END_REF]. Following several advantages listed in Sec. III-C3, we chose this strategy to solve the problem [START_REF] Ikemata | A physical principle of gait generation and its stabilization derived from mechanism of fixed point[END_REF].

D. Efficient optimal control solver

Our framework relies on the 20 years old optimization package MUSCOD-II [START_REF] Leineweber | An efficient multiple shooting based reduced sqp strategy for large-scale dynamic process optimization. part 1: theoretical aspects[END_REF], a multiple shooting solver for highly nonlinear systems submitted to path equality and inequality constraints, developed inside the Optimization and Simulation group at the University of Heidelberg. MUSCOD-II handles multi-phase systems with discontinue dynamics and periodic constraints with efficiency. While MUSCOD-II is a closed-source framework, one can depend on ACADO [START_REF] Houska | ACADO Toolkit -An Open Source Framework for Automatic Control and Dynamic Optimization[END_REF] which implements similar features.

IV. EXPERIMENTAL SETUP

The generality of our simulation framework is illustrated by various examples depicted in Tab. I. Those examples include different kinematic structures and different control schemes. Doing so, we may compare different polyarticulated topologies and, for a same topology, different control schemes, including either active or passive actuators.

Once a topology is chosen, we show that it is possible to replace the active actuation by a passive spring damper system. The motions of the first five robots are constrained to lie in the sagittal plane. Such restriction is removed for the robot with arms which is in 3D.

A. Inputs and Outputs

Fig. 1 shows the general inputs and outputs of our framework. In the experiments of Tab. I, the inputs of the system are the kinematic structure of the robot with the anthropometric parameters of the body segments (see Sec. IV-E), and the cost function, constraints and actuation type (see resp. Sections IV-B to IV-D).

In those experiments, we choose to set at a fixed value both the step duration (0.8 s) and the slope (0.05 rad) for each scenario, and let the solver optimize the step length. We limit the step length to the bounds [0.4; 1.] m.

The outputs of the system are the optimal cost of transport, together with the associated step length of the gait and the state and control trajectories of the joints during one step.

We also collect the number of iterations needed by the solver to converge and the total computation time for each experiment.

B. Cost function

We use the same cost function in each scenario. This cost function corresponds to the classic cost of transport (CoT). The CoT is a non-dimensional quantity which reflects the energy efficiency of a locomotion pattern. By definition, the CoT is the ratio between the energy consumed by the system (E) and its weight (mg) multiplied by the traveled distance (d):

CoT def = E mgd (8)
with g the gravity field value and m the mass of the system. The energy (E) consumed by the system is the sum of the potential energy mgh (with h the variation in altitude of the center of mass) with the integral of the power input

T t=0
q+ -q-M δ + τ • qdt. In this last formula, δ is the Dirac impulsion corresponding to the impact instant, • is the dot product operator and x M def = √

x M x. Finally, on a slope with angle α, the CoT is given by:

CoT = sin α + T t=0 q+ -q- M δ + τ • qdt mgd (9)
This cost function is only used as an example for the different case studies in Sec. V; in real life examples, the choice of a better suited cost function is still an open problem. Moreover, the CoT we give are not meant to be compared with examples outside of this paper.

C. Constraints

To reduce the dimension of the NLP, we compute a solution only for a half-step due to the periodicity and the symmetry between right and left segments. We then constrain the swing leg to be in contact with the ground only at the beginning and the end of the simulation. Then, the cyclicality Cases studies are based on six walkers models to measure the impact of the knees, torso, neck, arms and actuation type. For each example, the considered output are the cost of transport and the step length. The two last lines give the performance of the algorithm. In all examples the duration of a step is fixed (0.8s) as well as the slope of the ground (0.05rad).

INPUTS Model reference M A M B M C M D M E Model
and the symmetry of the motion is enforced by periodic constraints on both the initial and final configuration, velocity and torque of each joint trajectory.

D. Actuation

In a first stage, we propose an active actuation. The torque trajectory of each joint is modelized by piecewise cubic polynomials. An hardware implementation would then need an external source of power for the system, as a battery or an air canister, and an intelligent controller to deliver the desired torque.

In a second stage, we compare this active actuation with a passive one, which is simply a proportional derivative (PD) controller. In this case, we apply the torque τ j to the joint j: τ j = -K pj (q j -q 0j ) -K dj qj [START_REF] Benallegue | Contribution of actuated head and trunk to passive walkers stabilization[END_REF] where q j is the position of the joint and qj its velocity. K pj is the spring constant, q 0j the free length of the spring, and K dj the damping coefficient. Those three parameters are optimized by the solver. From those spring-damper parameters, one may craft a purely passive walker. In this case, the single source of power is the gravity field. Of course, another constraint for the numerical solver is to use the same spring and damper for each symmetrical joint.

E. Body parameters

All the robots have an anthropometric mass distribution, i.e. segment length, mass, center of mass position and inertia tensor are scaled according to a reference height and a reference mass. Those inertial parameters follow the anthropometric table proposed in [START_REF] Dumas | Adjustments to McConville et al. and Young et al. body segment inertial parameters[END_REF].

Our minimal model is composed of a pelvis, two thighs and two legs. Ankles (and knees on the model M B ) are actuated. On top of that, in the models M C , M D and M E , we add a torso and a head. This head is only actuated in the model M D . At last, we add two arms and forearms in the model M E with actuated shoulders. All actuated joints are revolute along parallel axes.

V. EXPERIMENTAL RESULTS

In the following, all the experiments are composed of 9 nodes for the single support phase and 1 single node for the impact phase.

A. Influence of the knees

The influence of the knees is highlighted by comparing models M A and M B in Tab. I.

The experimental results show that adding knees to a compass walker roughly divides by two the CoT while increases the optimal step length. The computational cost of adding knees is negligible.

B. Influence of the torso

Here, we study the impact of adding a fixed torso-head segment attached to the pelvis. This situation corresponds to the models M A and M C of Tab. I respectively, with a non-passive actuation.

Then, adding a fixed torso on top of the pelvis leads to a 40 percent decrease of the CoT and also reduces the optimal step length down to 40 cm, which is the lower bound we allowed. In the meanwhile, the computation times remain the same.

C. Influence of the neck

We consider the influence of actuating the neck by comparing the models M C and M D while keeping the actuation type.

The experiments on those models show that the increase of the CoT is lower than one percent even if we add an actuator. However, the solver needs more time to converge.

D. Comparison of actuation type

This experiment considers the model M C only and study the differences between active and passive controllers.

If we confront the results of the third and fourth columns of Tab. I, it appears that the cost is higher for the passive walker with respect to the choosen cost function, while the computational time is lower in the passive case. This can be explained by the dimensionality of the problems: in the first case, the dimensionality of the joint actuation is defined by four polynomial parameters times the number of shooting nodes while it is defined by three scalar parameters corresponding to the spring-damper model in the other case.

E. Beyond the 2D sagittal plane

Our simulation framework is general and it allows to address 3D model. In order to let the robot keep its balance, we add two actuated arms on the model M E .

In comparison to the 2D model without arms M C , the CoT is a bit higher, but remains lower than the first compass walker M A . We also notice that both the computation time and the number of iterations are higher, but the solver still converges in a similar time.

VI. DISCUSSIONS AND FUTURE WORKS

While our framework shares some common features with the Matlab package developed by Remy et al. [START_REF] Remy | A matlab framework for efficient gait creation[END_REF], its differs on various aspects. The first difference lies at the optimal control method level. In [START_REF] Remy | A matlab framework for efficient gait creation[END_REF], the authors use either collocation or shooting methods. Still, as mentioned in [START_REF] Diehl | Fast Direct Multiple Shooting Algorithms for Optimal Robot Control[END_REF] and recalls in Sec. III-C, the multiple shooting approach is the most suited in the context of multi-phase dynamical systems submitted to impacts and constraints. In addition, the choice of the C++ programming language to write both our dynamic simulation and optimal control problem is fundamental for efficient and fast computation.

Currently, we do not check the stability of the resulting motions, but only its balance at a dynamic level. As future extension, we plan to implement an approach similar to [START_REF] Mombaur | Open-loop stable solutions of periodic optimal control problems in robotics[END_REF] to converge on optimal and stable open-loop motions.

An other extension of this framework can be achieved at the level of contact modelling. Currently, punctual contacts are the standard for passive walkers. They are easy to simulate and to mechanically integrate as feet for passive walkers. However, rolling contacts may lead to more efficient gaits, as suggested by Kuo et al. [START_REF] Adamczyk | Mechanical and energetic consequences of rolling foot shape in human walking[END_REF]. We recently proposed in [START_REF] Carpentier | An analytical model of rolling contact and its application to the modeling of bipedal locomotion[END_REF] an analytical formulation of the rolling contact which is compatible with our dynamic formulation presented in Sec. II-D. The addition of this contact model inside our framework can largely improve the quality of the generated motions and enable us to deal with more complex scenarios.

Our framework is not restricted to passive walkers. On a wider scale, it can evaluate some paradigms in the design and the control of new humanoid robots [START_REF] Laumond | The Yoyo-Man[END_REF], [START_REF] Grizzle | 3d bipedal robotic walking: models, feedback control, and open problems[END_REF] and exoskeletons [START_REF] Koch | Exoopt -a framework for patient cetered design optimization of lower limb exoskelerons[END_REF]. In the near future, we will exploit this framework to evaluate and design a bipedal robot inspired from passive walkers.

Fig. 1 :

 1 Fig.1: Overview of the simulation framework. The simulator is described in Sec. II, and the solver is detailed in Sec. III.

TABLE I :

 I

The authors are with CNRS, LAAS, Univ. de Toulouse, UPS, INSA, 7 avenue du Colonel Roche, F-31400 Toulouse, France,

ACKNOWLEDGEMENT

This work is supported by the European Research Council (ERC) through the Actanthrope project (ERC grant agreement 340050), the European project KOROIBOT (FP7-ICT-2013-10/611909) and the French National Research Agency thought the Entracte project (ANR grant agreement 13-CORD-002-01).

We warmly thank Katja Mombaur and the Simulation and Optimization group at the Interdisciplinary Center for Scientific Computing (IWR) of Heidelberg University in Germany for providing the optimal control framework MUSCOD-II.


Deni Ferdian

Jacques Lacaze
email: bjacques.lacaze@ensiacet.fr

Evaluation of (Al)-Si Eutectic Reference Temperature of A3xx alloys

Keywords: Al-Si alloys, modification, thermal analysis, eutectic reference temperature

The modification level of Al-Si alloys is generally evaluated by the depression of the (Al)-Si eutectic temperature which can be recorded by thermal analysis. However, this method requires a reference temperature which should be the eutectic temperature evaluated on the relevant phase diagram. Various methods proposed to account for the effect of low level alloying elements on this reference temperature are reviewed and emphasis is put on the so-called "Mondolfo's equation" which is updated. Predictions are compared to experimental information from literature.

Introduction

Thermal analysis became a widespread technique for evaluating melt preparation of aluminum alloys in the 1980's [START_REF] Tuttle | Principle of thermal analysis for molten metal process control[END_REF]. It is in particular used for checking grain refinement and eutectic modification. Several parameters have been proposed in the literature for characterizing the thermal effects seen on the cooling curves which are associated to nucleation and early growth of (Al) phase (grain refinement) and to the (Al)-Si eutectic (eutectic modification) [START_REF] Apelian | Assessment of grain refinement and modification of Al-Si foundry alloys by thermal analysis[END_REF][START_REF] Gruzleski | The treatment of liquid aluminum-silicon alloys[END_REF]. Concerning eutectic modification, the most used method is to correlate the effectiveness of modification with increased eutectic undercooling or so called "eutectic depression". Evaluation of this undercooling requires the knowledge of a eutectic reference temperature which should be given by the relevant phase diagram and has often been evaluated experimentally as the eutectic temperature of the unmodified alloy. However, in many cases, this latter cannot be obtained in a cast shop due the possibility of prior modification treatment to the alloys or the use of returns in the charge.

One of the most used methods to evaluate the reference temperature of the (Al)-Si eutectic was proposed by Apelian et al. [START_REF] Apelian | Assessment of grain refinement and modification of Al-Si foundry alloys by thermal analysis[END_REF] who derived an equation expressing the eutectic temperature as a function of alloy's composition from the phase diagram compiled by Mondolfo [START_REF] Mondolfo | Aluminum alloys: structure and properties[END_REF]. The so-called "Mondolfo's equation" is still widely accepted as seen in the recent work by Wang and Lu [START_REF] Wang | Spheroidization of eutectic silicon in direct-electrolytic Al-Si Alloy[END_REF]. In the present work, we first reconsider the derivation of the equation and then update it according to more recent assessed phase diagram information. We then review other approaches proposed in the literature and finally compare them to available experimental data from literature.

Prediction of the solidification path of Al-Si-X alloys

Mondolfo's equation. The equation was derived to account for the effect of six elements (Cu, Fe, Mg, Mn, Ni and Zn) on the (Al)-Si eutectic with the binary (Al)-Si eutectic, e Si , set at w e,Si =12.5 wt.% Si and T e,Si =577°C. The Apelian's method considered that the effect of these elements should be additive, so that the effect of each one could be evaluated based on the related Al-Si-X ternary phase diagram. As an example, Fig. 1 shows schematically the projection of the liquidus surface of the Al-Si-Fe system (right part of the diagram) and the evolution of the (Al)-Si eutectic temperature (left part of the diagram), i.e. along the line e si -E Fe , where E Fe is the three phase invariant eutectic point between liquid, (Al), Si and β-Al 5 FeSi phase.

Assuming linearity of the temperature change (∆T e ) along the two-fold (Al)-Si eutectic line, we can write: . This is graphically represented with the arrow in the graph. As mentioned above, the method assumes further that the effect of alloying elements is additive, therefore the reference temperature for the start of the (Al)-Si eutectic reaction of an aluminum-silicon alloy is given as:

( ) ∑ ⋅ ⋅ - = X X X Si Si , e
Si , e e w a w w T T [START_REF] Apelian | Assessment of grain refinement and modification of Al-Si foundry alloys by thermal analysis[END_REF] with the sum is extended to all X alloying elements and the a X are defined as for Fe in equation [START_REF] Tuttle | Principle of thermal analysis for molten metal process control[END_REF].

Ternary Al-Si-X phase diagrams where X is Fe, Mg, Mn or Ni, are very similar, i.e with a e si -E X line that is at nearly constant Si. Data relative to these diagrams and relevant to the present work are listed in Table 1, where is the X weight content in solid (Al) in equilibrium with the ternary eutectic liquid and T E,X is the three-phase eutectic temperature of the Al-Si-X system. The partition coefficients between (Al) and the liquid, k X , have been evaluated with the ternary eutectic data and are also listed in the table. Mondolfo selected ternary phase diagrams mostly according to the extensive work by Phillips [START_REF] Phillips | Annotated equilibrium diagrams of some aluminium alloy systems[END_REF]. Recent assessment of the diagrams with Fe, Mg and Ni agree with this latter while an updated version was considered for the Al-Si-Mn [START_REF] Du | A thermodynamic description of the Al-Mn-Si system over the entire composition and temperature ranges[END_REF], Al-Si-Cu [START_REF] He | Experimental investigation and thermodynamic modeling of the Al-Cu-Si system[END_REF] and Al-Si-Zn [9,10] systems. For these three systems, original data according to Mondolfo's assessment are listed between brackets in Table 1.

While the method assumed the silicon content to be constant along the eutectic line, at 12.5 wt.%, a correction is proposed here which is an attempt to account for the difference in the silicon content of the binary and ternary eutectic invariant points. Accordingly, the a X value will be calculated in the present work as:

Si , e X , E Si E X X , E Si , e X w w w T T a ⋅ - = (3)
Using the partition coefficients in Table 1, one can calculate the solidification path of Al-Si-X ternary alloys according to lever rule and Scheil's model. This is illustrated in Fig. 2a in the case of Ni, where is seen that Scheil and lever rule solidification paths are nearly superimposed and lie very close to Apelian's method (arrow). These observations also apply to the other three elements, Mg, Fe and Mn. However, it is also clear from the graph in Fig. 2 that the method applies only to Al-Si alloys having their composition in the triangle Al-e Si -E X . Thus, this method should be restricted to alloys with a maximum content in X such that the (Al)-Si eutectic precipitates before any other secondary eutectic phase. This maximum, max X w , depends on the alloy's Si content and is roughly given as: * see text

X , E X X , E Si Si max X w w w w ⋅ = (4)
For alloys with low level in alloying elements such as A356 where the total content in elements other than Al and Si is less than 1 wt.%, the accuracy of Apelian's method is quite good in predicting the (Al)-Si eutectic reference temperature. However, the equation fails to predict it correctly for aluminum-silicon alloys with high Cu alloying content such as A319 [START_REF] Gruzleski | The treatment of liquid aluminum-silicon alloys[END_REF]. The reason for this is evidenced when considering the Al-Si-Cu phase diagram in Fig. 2b drawn according to He et al. [START_REF] He | Experimental investigation and thermodynamic modeling of the Al-Cu-Si system[END_REF]. It is seen that the (Al)-Si eutectic line moves far away to the left of the line at 12.5 wt. %Si along with increasing copper content. This leads to copper contents as estimated by the method much higher than the one when the solidification path reaches the (Al)-Si eutectic line. The maximum overestimation of the copper content is obtained for the ternary eutectic point when the method would give 2.5 (i.e. 12.5/5.0) times too high copper content. This means that the a Cu coefficient as listed in Table 1 as "original factor" should be multiplied by 0.4. Because the actual copper contents in the Al-Si alloys are much lower than the value at the ternary eutectic point, the correction can not be that dramatic in all practicality. With the experimental data analyzed below, it has been found that a multiplication factor of 0.75 gives appropriate results, i.e. a Cu was set to 1.65 after accounting for the effect of Si ternary eutectic content, see Table 1. The Al-Si-Zn phase diagram is somehow similar to the Al-Si-Cu one in that the ternary invariant eutectic is located far away from the binary (Al)-Si eutectic, and in fact is very close to the Zn-rich corner of the phase diagram. For this system and owing to the low Zn levels in A3xx alloys, it seemed wiser to use the reported isopleth section at 5.3 at.%Si [9,10] to estimate the effect of Zn on the two-fold (Al)-Si saturation line. The corresponding data is shown in Table 1 where it is seen that the final estimate of a Zn differs greatly from the one assessed with the Apelian's method.

To sum up, the following equation is proposed that incorporates the changes of silicon content along the eutectic lines as well as an updating of the ternary phase diagrams:

( ) ( )
⋅ + ⋅ + ⋅ + ⋅ + ⋅ + ⋅ ⋅ - = (5)
This equation itself shows it is restricted to silicon content higher than 1 wt.%. Therefore, the equation should be limited to silicon contents in the range 1 to 12.5 wt.% and to alloying additions that depend on the silicon content, see Eq. 4, and are at most (wt.%) 4.9 for Mg, 0.7 for Fe, 5 for Cu, 5 for Zn, 0.4 for Mn and 5 for Ni. Other equations. There are several other approaches proposed to calculate the reference eutectic temperature. Studies by Joenoes and Gruzleski [START_REF] Joenoes | Magnesium effects on the microstructure of unmodified and modified Al-Si alloys[END_REF] which focused on the magnesium effect, proposed a series of empiric calculations with a coefficient depending on the Si content in the alloys. While those equations only consider one element in the calculation and are limited to one type of chemical composition, Sthuldreier et al. [START_REF] Stuhldreier | Erfahrungen mit thermischen analyse von G-AlSi-legierungen[END_REF] considered 3 major elements, Mg, Cu and Fe based on their experimental data. A different approach was used by Djurdjevic and collaborators [START_REF] Robles Hernandez | Calculation of the liquidus temperature for hypo and hypereutectic aluminum silicon alloys[END_REF][START_REF] Djurdjevic | Thermal description of hypoeutectic Al-Si-Cu alloys using silicon equivalency[END_REF] who described the (Al) liquidus with second order polynomial for binary Al-X systems of interest and then defined a silicon equivalent (Si eq ) for each element X. Then, the enrichment of the liquid is described as with the Apelian's method and the eutectic temperature writes as follows:

( )

        ⋅ ⋅ + ⋅ - = Si 2 eq

Result and Discussion

A comparison of reference temperatures predicted with equations from the present work, Apelian et al. [START_REF] Apelian | Assessment of grain refinement and modification of Al-Si foundry alloys by thermal analysis[END_REF], Stuhldreier et al. [START_REF] Stuhldreier | Erfahrungen mit thermischen analyse von G-AlSi-legierungen[END_REF] and Djurdjevic [START_REF] Djurdjevic | Thermal description of hypoeutectic Al-Si-Cu alloys using silicon equivalency[END_REF], with thermal analysis data from literature [START_REF] Apelian | Assessment of grain refinement and modification of Al-Si foundry alloys by thermal analysis[END_REF][START_REF] Djurdjevic | Thermal description of hypoeutectic Al-Si-Cu alloys using silicon equivalency[END_REF][START_REF] Charbonnier | Microprocessor assisted thermal analysis testing of aluminum alloy structure[END_REF][START_REF] Gowri | Effect of alloying elements on the solidification characteristics and microstructure of Al-Si-Cu-Mg-Fe 380 alloy[END_REF][START_REF] Djurdjevic | Online prediction of aluminum-silicon eutectic modification level using thermal analysis[END_REF][18][19][START_REF] Farahany | Evaluation of the effect of Bi, Sb, Sr and cooling condition on eutectic phases in Al-Si-Cu alloy (ADC12) by in situ thermal analysis[END_REF][START_REF] Han | Microstructure characteristics in nonmodified and Sr modified Al-Si-Cu-Mg 319 type alloys[END_REF][START_REF] Thompson | Advanced light metals casting development: solidification of aluminium alloy A356[END_REF][START_REF] Mcdonald | Eutectic grains in unmodified and strontium-modified hypoeutectic aluminum-silicon alloys[END_REF][START_REF] Djurdjevic | Synergy between thermal analysis and simulation[END_REF][START_REF] Djurdjevic | Characterization of the solidification path of AlSi5Cu(1-4 wt.%) alloys using cooling curve analysis[END_REF] is shown in Figure 3a. In addition, comparison with calculations made with Thermo-Calc and the TCAL1 database [START_REF]Thermo-Calc Software, TCAL1 database version 1.0[END_REF] according to Scheil's model are also shown in Fig. 3b. In the literature data, two different methods are used to extract the eutectic temperature from the cooling curves which relate to recording the plateau temperature, T E,G , and the nucleation temperature, T E,Nuc respectively. The method using the nucleation temperature gives systematically much higher temperatures than that predicted by CALPHAD calculations. Though such an observation may need further investigation, it suggests that the change in the cooling rate as recorded with thermal analysis cups and associated to nucleation in the central part of the cup, may relate in fact to the heat released by the outer part of the cup where the eutectic transformation has started. In a way, the reading of thermocouple would be influenced by the heat flux released from the outer shell of the thermal analysis cup.

The graph shows that at high eutectic temperature with low alloying additions, equations by Apelian et al. [START_REF] Apelian | Assessment of grain refinement and modification of Al-Si foundry alloys by thermal analysis[END_REF] and Stuhldreier et al. [START_REF] Stuhldreier | Erfahrungen mit thermischen analyse von G-AlSi-legierungen[END_REF] present a low deviation from the bisector line. However, their predictions become scattered as the eutectic temperature is lowered due to higher alloying element. Moreover, the silicon equivalent method proposed by Djurdjevic [START_REF] Djurdjevic | Thermal description of hypoeutectic Al-Si-Cu alloys using silicon equivalency[END_REF] predicts values significantly above the experimental results in all the composition ranges investigated. [12] and Djurdjevic [START_REF] Djurdjevic | Thermal description of hypoeutectic Al-Si-Cu alloys using silicon equivalency[END_REF]. The dotted line are bisectors.

Conclusion

The so-called Mondolfo's equation derived by Apelian et al. [START_REF] Apelian | Assessment of grain refinement and modification of Al-Si foundry alloys by thermal analysis[END_REF] does agree quite well with the experimental and CALPHAD results for low Cu alloys, while giving too low values otherwise. The present approach that accounts empirically of the effect of Cu does closely follow CALPHAD-type predictions. It is finally seen that the proposed equation shows a better agreement with the experimental data than other equations in the whole range of compositions.

Fig. 1 .

 1 Fig. 1. Schematic Al corner of the ternary Al-Si-Fe system

Fig. 2 .

 2 Fig. 2. Solidification path plotted onto the liquidus projection of the ternary Al-Si-Ni (a) and Al-Si-Cu (b) phase diagrams. Arrows illustrates Apelian's method, while Scheil's and lever rule paths are hindered below them.

Fig. 3 .

 3 Fig. 3. Graph showing the correlation between experimental data (a), TCAL1 database (b) and calculated eutectic temperatures according to present work, Apelian et al. [2], Stuhldreier et al.[START_REF] Stuhldreier | Erfahrungen mit thermischen analyse von G-AlSi-legierungen[END_REF] and Djurdjevic[START_REF] Djurdjevic | Thermal description of hypoeutectic Al-Si-Cu alloys using silicon equivalency[END_REF]. The dotted line are bisectors.

Table 1 .

 1 Ternary eutectic points data of selected Al-Si-X ternary systems. See the text for definitions; compositions are in wt.% and temperature in Celsius. Values assessed by Mondolfo are between brackets when they have been updated.

	Element	w	, Si E	X	w	, X E	S	w	, X E	(	) Al	k X	T E,X	Original factor a x	w w	Si , X , E e Si	Present work factor a x
	Mg	12.95	0.85	4.96	0.171	555	4.43	1.036	4.59
	Fe	12	0.05		0.7	0.071	576	1.43	0.96	1.37
	Cu*	6 (5)	4.5 25.6 (27) 0.175 521 (525)	1.93	0.75	1.65
	Zn*	4.2 (3)				40 (80)		535 (440)	1.7	0.336	0.35
	Mn	12. 7 (12) 0.03	0.4 (1) 0.075 576 (574)	3	1.016	2.54
	Ni	11	0.05		5		0.01	557	4	0.88	3.52

Acknowledgments

The authors wish to acknowledge the financial support received from the European Space Agency (contract #4200014347) through the MICAST program and DGHE -Indonesia Education Ministry.


Wissam Sabbah

Pierre Bondue

Oriol Avino-Salvado

Cyril Buttay

Hélène Frémont

Alexandrine Guédon-Gracia

Hervé Morel

High Temperature Ageing of Microelectronics Assemblies with SAC Solder Joints

In some applications, electronic systems are expected to operate at high ambient temperature (e.g. 150 °C). In this paper, we investigate the failure mechanism and microstructure evolution of solder-free (SAC) solder joints at a maximum temperature of 175 °C. It is found that no new failure mechanisms are triggered, and that ageing tests for solder can be accelerated at 175 °C. In particular, the growth rate of the interfacial intermetallic compound (IMC) is found to be consistent with that observed at lower temperatures.

Review of high-temperature microelectronics assemblies

Electronic systems are increasingly used in harsh environment. For example, in automotive applications, some systems are expected to operate in close proximity to the engine, at ambient temperatures exceeding 120 °C [START_REF] Kassakian | The Future of Electronics in Automobiles[END_REF]. Other examples include the control of power electronic systems [START_REF] Buttay | State of the art of high temperature power electronics[END_REF]. Therefore, there is a need for microelectronics assemblies that are able to operate at 150 °C or more.

Some applications (such as down-hole data logging or aerospace) involve very high ambient temperatures (in excess to 200 or even 250 °C). They require specific assembly technologies, such as ceramic substrates, high temperature solders, etc. In particular, component attaches rely either on gold-based solder alloys, or on silver sintering [START_REF] Manikam | Die Attach Materials for High Temperature Applications: A Review, Components, Packaging and Manufacturing Technology[END_REF][START_REF] Khazaka | Survey of High-Temperature Reliability of Power Electronics Packaging Components, Power Electronics[END_REF].

This paper focuses on temperatures ranges up to 175 °C, where it is desirable and conceivable to use more standard technologies such as PCBs and lead-free solders. These technologies are typically used at lower temperatures, and most papers only address temperatures up to 125 °C [START_REF] Hai | Effects on the reliability of lead-free solder joints under harsh environment[END_REF]. In some rare cases, the authors go up to 150 °C [START_REF] Berthou | Microstructure evolution observation for SAC solder joint: Comparison between thermal cycling and thermal storage[END_REF].

Therefore, there is a need to investigate the behaviour of lead-free solders at 175 °C. SnAg 3 Cu 0.5 (SAC305) solder is considered in this article. In particular, we will study the growth of the Intermetallic Compounds (IMCs) [START_REF]Lead-Free Electronic Solders[END_REF], which have an effect on the thermo-mechanical strength of the solder joints [START_REF] Pin | Numerical study of thermomechanical fatigue influence of intermetallic compounds in a lead free solder joint[END_REF].

In this paper, we present such investigations, using two test conditions: temperature cycling, and isothermal ageing. They are described in the next section. The results are presented in section 3 and discussed in section 4.

Test Protocol

The results presented in this paper were produced by two different teams. Schematically, one of the teams (Laboratoire * Corresponding author Ampère, Lyon) performed high temperature storage tests, while the other (IMS, Bordeaux) was in charge of the temperature cycling tests. Special care was taken at the analysis stage to produce consistent and comparable data.

Sample Preparation

The results are based on Printed Circuit Board (PCB) material, on which ceramic capacitors are assembled with SAC 305 solder joints. The board material is polyimide (Arlon 85N), with a very high glass transition temperature (> 250 °C), and 6 layers of copper with a thickness of 70 µm. Pad finish is ENIG (Ni: 6 µm, Au: 80 nm). The boards were assembled by the company Elvia PCB (France).

The capacitors are Vishay VJ series (1210 package), with Ni terminals plated with 100 % Sn. The solder attach is a lead-free alloy (SAC305, HFA8.9 Series from Indium Corp.). Soldering was performed by the company Tronico (France). No cleaning of flux residues was performed afterwards. Ceramic capacitors were chosen for this study because they are fragile and large components, so they tend to be sensitive to thermomechanical stresses.

Before testing, the boards were cut into smaller coupons, using a diamond cut-off wheel (3000 RPM, 1 mm/s feed speed). No degradation of the PCB was observed (no burring of copper layers nor any trace of delamination). After cutting, the coupons were cleaned in water and ethanol, and dried to remove any lubricant residue. Before starting the actual tests, the samples were dried at least 6 h at 120 °C [START_REF]IPC, Ipc-1601a -printed board handling and storage guidelines[END_REF][START_REF] Charrier | Moisture pick-up and minimum drying conditions of epoxy and polyimide-glass printed circuit board[END_REF].

Characterization

The mechanical strength of the solder joints was measured through shear testing of the capacitors (5 capacitors were used for each data point). At Ampère, this was performed using a XYZTEC Condor Sigma tester, equiped with a 200 kgf load cell, a tool height of 350 µm and a tool speed of 200 µm/s. At IMS, a Dage 3000 tester, with a 100 kgf load cell, a tool height of 200 µm and a tool speed of 100 µm/s was used. The shear stress was calculated by dividing the rupture shear force by the maximum surface area of the solder joint (4.35 mm 2 per pad).

For each data point, two capacitors (i.e. 4 solder joints) were cross-sectioned and prepared for observation of the solder micro-structure. Once encapsulated and cut using a low-speed diamond saw, the samples were prepared as follows:

• at Ampère: Grit 1200 SiC paper (60 s), 10 µm diamond disk (120 s), followed by 3 and 1 µm diamond suspension on hard clothes (180 s each), and a coloidal silica finish on soft cloth (60 s). Final cleaning in water and ethanol;

• at IMS: Grits 1200 (150 s), 2400 (180 s) and 4000 (240 s) SiC paper, followed by 3 (240 s), 1 (300 s) and 1/4 µm (180 s) diamond suspensions on hard clothes, and a colloidal silica (45 s) on soft cloth. This is specified in [START_REF] Berthou | Microstructure evolution observation for SAC solder joint: Comparison between thermal cycling and thermal storage[END_REF].

Observation of the micro-structure was performed mainly using SEMs (Tescan Vega 3 at Ampère, JEOL JSM-6100 at IMS, both equiped with EDS sensors for elemental analysis). One particular aspect was the measurement of the Intermetallic Compounds (IMC) layer thickness at the solder/capacitor and solder/PCB interfaces. For this, an interface length of 30-50 µm was observed. At Ampère, we performed a series of 20 to 30 manual thickness measurements over this length. At IMS, the area of the IMC was identified through image processing (brightness threshold), and this area was divided by the length of the interface. Initial comparisons produced similar results for both methods.

Finally, electrical tests (capacitance measurements) were performed regularly at IMS to detect any damage to the capacitors. At Ampère, the assemblies were photographed before and after tests, using a Leica M205C binocular and a camera which settings were recalled to ensure consistent color rendering.

Reliability tests conditions

As a preliminary test, un-populated PCBs coupons were submitted to "step stressing": these coupons were placed in an oven, at a temperature of 150 °C. This temperature was increased by 20 °C every 24 h, up to 250 °C. At each step, one coupon was removed from the oven, and characterized. In particular, we measured its change in weight (Ohaus Adventurer Pro 210 scale, with a 1 mg resolution).

Fig. 1 shows that no noticeable weight loss was observed at 210 °C and below. A small weight loss (0.7 %) was observed after an additional 24 h step at 230 °C, and a slightly larger weight loss (1.5%) after 24 h at 250 °C. Overall, this weight loss remains very limited compared to that observed on FR4 boards in the same conditions (8% at 230 °C). This is confirmed by a second series of tests (not shown here), where PCBs using the same Arlon 85N material were stored in air for 1000 h at 190 °C without any noticeable weight loss. Such results are consistent with tests performed at higher temperature [START_REF] Shaddock | Reliability of high temperature laminates[END_REF]. Therefore, one can consider that the tests presented hereafter (around 3000 h at 175 °C) should not result in degradation of the PCB material.

Based on these preliminary results, it was decided to conduct the reliability tests using the following parameters:

• Cycling tests: using a Rapid Temperature Change (RTC) chamber (Climats Excal 1423-TE), with temperature ramp-up and ramp-down slopes of 10 °C/min, a minimum temperature of -55 °C and a maximum temperature of 175 °C (both maintained for 15 min).

• Storage tests: in a forced convection oven (Nabertherm N30/45HA), at a constant temperature of 175 °C.

Solder Joint

Results

The shear-test measurements are presented in Fig 2 . Initially, the solder attach has a very good strength, as in most cases the rupture occurs partially in the capacitors. For the storage tests, the mechanical strength slowly decreases from ≈ 15 MPa at t=0 down to ≈ 9 MPa at t=3000 h.

For the cycling tests, however, the change is much more dramatic, with a strong reduction (80%) in shear strength after 2250 cycles (equivalent to 562 h at 175 °C).

Cross-sections of some solder joints are presented in Fig. 3 (optical microscopy). In Fig. 3a, the solder joint was not submitted to any cycling test (initial conditions). This cross section is representative of the solder attaches, with some voiding. Indeed, a noticeable number of voids was found in most solder joints.

After 750 cycles (Fig. 3b), a crack is visible in the solder layer under the capacitor. Another observation is the change in the shape of the solder joint. As cycling progresses, the solder boundary roughens and becomes thinner on top of the meniscus. SEM observations (Fig. 4) show a change in the interfacial microstructure of the solder joint: in initial conditions (Fig. 4a), the IMCs form a lamellar structure. After 310 h storage (Fig. 4b), the interfacial IMCs become continuous. In cycling tests (Figs. 4c and4d for substrate and capacitor side, respectively) the interfacial IMCs also become continuous, but some voiding appears at the IMC/bulk solder interface. No such crack is observed in storage (310 h and 1500 cycles are equivalent in terms of time spent at 175 °C). This indicates that the cracks are generated by the thermo-mechanical stresses. This also means that the continuous IMC layers have little influence on the mechanical strength, as presented in Fig. 2.

For longer storage periods, however, cracks start to appear between the capacitor and the substrate (Fig. 5). These cracks develop where the solder is the thinnest, and occur once the IMCs have grown to occupy the entire solder volume. However, this cracking has little influence on the shear strength, as visible in Fig. 2 (no drop in shear strength for samples stored for 984 h). This can be explained because the solder area under the capacitor terminals only represent a small part of the total joint. Another visible effect of storage ageing is the cracks observed on some capacitors (Fig. 5b). However, these are not observed on all capacitors. Only two capacitors were cross-sectioned per time step, so no statistics could be made regarding the effect of storage time on capacitor cracking.

As storage time increases, the IMCs continue to grow, and voids develop all along the IMC/solder interface (Fig. 6), as was observed for the cycled samples in Fig. 4d. This is correlated with the continuous decrease in the shear strength, as plotted in Fig. 2.

The thickness of the IMC layers, measured at the solder/substrate interface and at the solder/capacitor interface is shown in Fig. 7, as a function of storage time. The difference in the growth rate is related to the finish of the parts: the capacitors have a Sn finish, while the PCB has a NiAu finish.

Examples of some rupture surfaces observed for capacitors stored for 2062 h at 175 °C are given in Fig. 9. Several regions can be identified in Fig. 9c: region A corresponds to the solder located below the capacitor electrode. In this area, the shearing of the capacitor caused no smearing, indicating that the solder was completely separated from the capacitor electrode, as can be observed in Figs. 5 and6. In region B, which corresponds to the "vertical" solder/capacitor interface in Figs. 5 and6, some smearing can be observed. This indicates that the interface was not completely cracked. Also, many voids can be observed. Region C corresponds to the solder area located below the capacitor ceramic body (where no attach occured before or after ageing). Finally, in the right side, some leftovers of the capacitor electrode (Ni) can be seen. The solder layer underneath was larger (because the capacitor has rounded edges), and did not crack completely as the IMC grew, resulting in a better solder/metallization interface.

Discussion

In thermal cycling, the capacitors showed a very fast reduction in the shear strength of their attach (6-fold reduction in 500 None of the the polyimide PCBs was found to be damaged by the long-term exposure to 175 °C. This is consistent with the specifications of the material, and with the preliminary tests presented in Fig. 1 with less than 1% weight loss after 24 h at 230 °C. All damages were found to be in relation with the solder layer: cracks propagated either at the solder interfaces, or in the capacitors, close to the solder joint.

The evolution of the solder interfacial microstructure is as follows: interfacial IMCs grow at different speeds at the capacitor and PCB interfaces; cracks propagate at the solder/IMC interface below the capacitor, where the solder joint is thinner and depletes quicker (Fig. 5); Finally the cracks develop along the complete Sn/IMC interface (Fig. 6). The IMC growth rate is plotted in Fig. 7. Previous tests [START_REF] Berthou | Microstructure evolution observation for SAC solder joint: Comparison between thermal cycling and thermal storage[END_REF] were performed at lower storage temperatures (80, 125 and 150 °C). After 1000 h, the IMC thicknesses were found to be 2.25, 2.5 and 3.6 µm respectively on a PCB with NiAu finish. This is consistent with the with SAC305 balls mounted on small PCB pads. Also, in our case, the cracks which appeared in the capacitors attach limited the growth of the IMC. With a different geometrical configuration, where the solder joint would be allowed to change in thickness without cracking, the mean IMC growth rate might be found to be faster. The cracks occurred both in cycling and storage conditions, but they appeared sooner for cycling, and are related to two different mechanisms. For cycling, the rupture occurred along the brittle IMCs regardless of their position (Figs. 4c andd). In storage, they started under the capacitor electrodes, and were caused by the depletion of the solder joint due to the interfacial IMC growth.

Conclusion

In this article, we assessed the ageing of SAC305 solder joint submitted to cycling or storage conditions with a maximum temperature of 175 °C. It was found that the joints degrade rapidly, especially in cycling conditions. Except when short term operation is allowed (for example for testing purposes), the reliability of these joints is probably insufficient for most applications.

However, the ageing tests performed at 175 °C did not show any new failure mechanism compared to previous results obtained at lower temperatures. Therefore, it is possible to run ageing tests at 175 °C on SAC soldered assemblies to perform a faster reliability assessment.

The growth rate of the interfacial IMCs is consistent with the rates measured at lower temperatures, albeit on a different structure (solder bumps of a BGA package). A more accurate investigation would require testing of the same structure at different temperatures. This would allow to identify the activation energy of the ageing mechanism.

Figure 1 :

 1 Figure1: Weight change of unpopulated PCB test coupons during "step stress": each successive temperature step is maintained for 24 h, and one test coupon is removed from the oven at each step, for characterization (including weight measurement).

Figure 2 :

 2 Figure 2: Shear strength of the capacitors as a function of the time spent at 175 °C, for test vehicles submitted to storage or cycling conditions. The error bar corresponds to the standard deviation of the measurements.

Figure 3 :

 3 Figure 3: Cross section of solder joints, before testing (a) and after 750 cycles between -55 °C and 175 °C (b). A few voids are visible initially in Fig. (a). After cycling, a crack develops, first under the capacitor terminals, and eventually across the full solder joint.

Figure 4 :

 4 Figure4: Enlargement of the solder area between the capacitor terminal and the PCB substrate. From top to bottom: cross section before testing (0 h), cross section after 310 h storage at 175 °C, substrate/solder interface after 1500 cycles (corresponding to 375 h at 175 °C), and capacitor (left)/solder interface, also after 1500 cycles.

Figure 5 :

 5 Figure 5: Microsections after 984 h storage at 175 °C. The growth of the IMC generated stresses that resulted in crack, either in the solder joints, in the capacitor metal finish, or in the capacitor body itself.

Figure 6 :

 6 Figure 6: Microsections of the solder joints after 2562 h storage at 175 °C, showing the propagation of the voiding along the IMC/solder interface.

Figure 8 :

 8 Figure 8: Enlargement of the Substrate/Solder interface after 3000 h storage at 175 °C. Many cracks appear at the interfacial IMC area. Ag forms precipitates (clear areas) in the bulk of the solder joint.

Figure 9 :

 9 Figure 9: Rupture surface on one side of the capacitor (sample stored for 2062 h at 175 °C): Observation area; overview of the breaking facies after shear test of the capacitor (travel direction of the shear tool was from right to left), zoom on the right area of the rupture surface.

Acknowledgement

This work was funded as part of the GENOME-PREMICES initiative, in the scope of the "Programme d'Investissements d'Avenir" (Investment for the future) of CORAC (French Research Council for Civil Aeronautics).