J D Frankland
email: john.frankland@ganil.fr

D Gruyer

E Bonnet

B Borderie

R Bougault

A Chbihi

J E Ducret

D Durand

Q Fable

M Henri

J Lemarié

N Le Neindre

I Lombardo

O Lopez

L Manduci

M Pârlog

J Quicray

G Verde

E Vient

M Vigilante

Model independent reconstruction of impact parameter distributions for intermediate energy heavy ion collisions

We present a model-independent method to reconstruct the impact parameter distributions of experimental data for intermediate energy heavy ion collisions, adapted from a recently proposed approach for ultrarelativistic heavy ion collisions. The method takes into account the fluctuations which are inherent to the relationship between any experimental observable and the impact parameter in this energy range. We apply the method to the very large dataset on heavy ion collisions in the energy range 20-100 MeV/nucleon obtained with the INDRA multidetector since 1993, for two observables which are the most commonly used for the estimation of impact parameters in this energy range. The mean impact parameters deduced with this new method for "central" collisions selected using typical observable cuts are shown to be significantly larger than those found when fluctuations are neglected, and as expected the difference increases as bombarding energy decreases. In addition, we will show that this new approach may provide previously inaccessible experimental constraints for transport models, such as an estimation of the extrapolated mean value of experimental observables for b = 0 collisions. The ability to give more realistic, model-independent, estimations of the impact parameters associated to different experimental datasets should improve the pertinence of comparisons with transport model calculations which are essential to better constrain the equation of state of nuclear matter.

I. INTRODUCTION

The equation of state (EoS) of bulk nuclear matter in a wide range of densities, temperatures, and proton-neutron asymmetries is of major importance not only for nuclear physics but also astrophysics since the EoS plays a fundamental role in the understanding of core-collapse supernovae (CCSN), proto-neutron star cooling [1], and neutron star mergers as recently observed through gravitational wave data [2]. In the laboratory, precise constraints on the finite temperature EoS away from saturation density can be obtained from heavyion collisions (HIC), which can be used to explore a wide range of density, energy and asymmetry conditions, depending on the bombarding energy, neutron and proton numbers (N, Z) of projectile and target, and impact parameter, b. In particular, the freeze-out stage of central HIC in the energy range 20-100 MeV/nucleon produces transient finite systems with similar temperatures and densities as CCSN matter, albeit with smaller asymmetries [3,4].

The EoS cannot be measured directly but has to be inferred by comparing the outcome of carefully selected colli-sions with model predictions. For example, transport models [5][6][7][8][9][10][11][12][13][14] can be used to predict the dynamics of collisions at different impact parameters using different forms of the EoS, by varying the parameters of the employed effective interaction. For such a comparison to be meaningful, the calculations should ideally be run over the same impact parameter distribution P(b|S) as that of the experimental event sample S, or at least use a representative value such as the mean of this distribution, b S . This is all the more important as currently transport model calculations are not only dependent on the EoS parameters, but also on many other ingredients (related to uncertainties both of physics and of numerical implementation) which are as yet not fully under control [15,16], therefore reducing trivial bias due to mismatching of experimental and simulated impact parameters is essential to make progress.

Of course, the impact parameter for each collision cannot be measured but only inferred from the final state observables of each event. An essential feature of HIC in the 20-100 MeV/nucleon energy range are fluctuations. The fluctuations of any observable X from one collision to the next can be of the same order of magnitude as its mean value, X . Due to a combination of diminishing cross-section and increasing importance of fluctuations as b → 0, higher and higher cuts in any observable, even one strongly correlated with b, will al-ways lead to mixing between collisions over a broad range of impact paramter [17,18]. Impact parameter mixing in experimental data samples of "central" collisions is a general feature which can seriously bias comparisons with transport models if neglected [19]. For example, the authors of [20] studied how mixing can affect studies of isospin sensitive variables using the ImQMD transport model. They show that in reality the mean impact parameter selected by higher and higher cuts on the total charged particle multiplicity "saturates" at a finite value b > 0, and the difference between estimated and true centrality worsens as bombarding energies decrease below 70 MeV/nucleon (see also for example [21]). Given the current uncertainties in existing transport models, it would be an important advance to be able to quantitatively characterize the centrality of selected event samples in a model-independent way.

There are many works in the literature dealing with the characterization of experimental centrality, from the simplest and most widely used geometric approximation of [22] to an AI approach using neural networks trained by simulated data [23]. Recently, a new method for reconstructing experimental impact parameter distributions was proposed for ultrarelativistic collisions in [24,25]. It is model independent and explicitly takes into account the fluctuations in the relationship between X and b, which is adjusted in order to reproduce the experimentally measured inclusive P(X) distribution. Despite the orders of magnitude differences in beam energies and physics, the inclusive P(X) distributions for different observables used to gauge the centrality of e.g. Au+Au collisions at center of mass (c.m.) energies √ s=130 GeV or Pb+Pb collisions at √ s =2.76 TeV (see Fig. 1 of [24]) have generic properties very similar to those seen for e.g. total multiplicity or total transverse energies in HIC collisions at intermediate energies: the highest cross-section for the smallest values of X, decreasing to a wide plateau and finally a near-exponential fall-off for the largest X values. Indeed such a generic distribution is expected for any observable whose mean value decreases monotonically with b when weighted with a geometric impact parameter distribution, P(b) ∼ b.

In this paper, our aim is not to criticize or improve existing methods of selecting experimental samples of central collisions. Rather, we propose a new method to characterize any selected set of experimental data in terms of the corresponding impact parameter distribution. The method is model independent and uses experimental data as its sole input. It is based on the approach of [24,25] which we will first present along with the adaptations we have made for its use in the intermediate energy range. Before applying it to our experimental data, we will first of all validate the method with a full simulation of a typical near-Fermi energy heavy-ion reaction using the microscopic transport model AMD [26], the statistical decay code GEMINI++ [27,28] and a software filter reproducing the characteristics of the INDRA array [START_REF]KaliVeda heavy-ion analysis toolkit[END_REF][START_REF] Pouthas | [END_REF]. Then we will apply it to the very large INDRA dataset on heavy ion collisions in the energy range 20-100 MeV/nucleon. We will use it to provide a systematic estimate of the true centrality of the collisions selected by a high-transverse energy cut, as often used in previous analyses [31]. In conclusion, by providing more reliable, model-independent estimates of impact parameters associated to different experimental data, this method will improve the quality of data-model comparisons which are essential to better constrain the equation of state of nuclear matter.

II. METHOD

A. General approach

Consider an observable X whose functional dependence on the impact parameter can be written in terms of a conditional probability distribution P(X|b), which encodes both the bdependence of the mean value, X(b), and the fluctuations of X about this mean value. The inclusive distribution of X resulting from all measured collisions, having an unknown impact parameter distribution P(b), is given by

P(X) = ˆ∞ 0 P(b) P(X|b) db (1)
Introducing the quantity centrality, c b , defined as the cumulative distribution function of P(b),

c b ≡ ˆb 0 P(b ) db (2)
Eq. 1 can be rewritten as

P(X) = ˆ1 0 P(c b ) P(X|c b ) dc b = ˆ1 0 P(X|c b ) dc b (3)
as by definition P(c b ) = 1, ∀c b , and the dependency on the unknown inclusive impact parameter distribution disappears. Then Eq. 3 can in principle be used to determine P(X|c b ) by fitting the experimentally measured inclusive distribution of X, P(X). In order to make the problem tractable, the authors of [24,25] proposed to write P(X|c b ) as

P(X|c b ) = f X(c b ), θ (c b ) (4)
where f is a suitable p.d.f. with mean value X(c b ) and reduced variance determined by θ (c b ) = var(X)/X(c b ). In principle, both the mean value and the variance of the observable can depend on centrality, i.e. the impact parameter.

Once P(X|c b ) is determined by fitting the experimental P(X) distribution, it can be used to estimate the centrality distribution for any given selection of data. For example, for a cut such as x 1 ≤ X ≤ x 2 the corresponding centrality distribution is given by

P(c b |x 1 ≤ X ≤ x 2 ) = ´x2 x 1 P(c b , X) dX ´x2 x 1 P(X) dX = ´x2 x 1 P(c b |X)P(X) dX ´x2 x 1 P(X) dX
where, in the first integral, P(c b , X) is the joint probability distribution, P(c b , X) = P(c b |X)P(X) = P(X|c b )P(c b ). The second equality being nothing but Bayes' theorem, we can use it to rewrite the integrand in the numerator, remembering that P(c b ) = 1. If we also define an experimental centrality, c x , as

c x ≡ ∞ x P(X) dX (5)
i.e. the fraction of all events with X ≥ x, then we can also rewrite the denominator, giving finally for the centrality distribution corresponding to our event selection

P(c b |x 1 ≤ X ≤ x 2 ) = ´x2 x 1 P(X|c b ) dX c x 1 -c x 2 (6)
More generally, events may be selected in many different ways, not only using simple cuts, and not necessarily using the same observable X as that which is used for centrality estimation: for a generic experimental sample S the centrality distribution will be given by

P(c b |S) =
´P(X|c b ) P(X|S) P(X) dX ´P(X|S) dX (7) where P(X|S) is the sample distribution of X (i.e. a histogram of X filled from the events in the sample), and the integrals are over the full domain of X. Therefore once P(X|c b ) has been determined by fitting the experimental P(X) distribution, it can be used to estimate the centrality distribution for any given selection of data. Finally, absolute impact parameter distributions can be deduced from the calculated centrality distributions by a suitable change of variables:

P(b|S) = P(b)P(c b (b)|S) (8)
It should be noted that in this case it is necessary to assume a specific form for the inclusive impact parameter distribution P(b) and calculate the corresponding relationship between c b and b, c b (b), using Eq. 2.

B. Specific implementation

In order to apply the method to data, specific choices and approximations must be made. These concern the p.d.f. f to use in Eq. 4, the parametrization of the centrality dependence of the mean, X, and the centrality dependence of the reduced fluctuation, θ .

Fluctuation kernel, f

As observables most strongly correlated with impact parameter in heavy-ion collisions can be assumed to result from a sum of many independent microscopic processes (e.g. nucleon-nucleon collisions), a natural choice would be a gaussian or normal distribution [24]. However, as most observables related to collision violence only take positive values, X ≥ 0, the gaussian distribution has the disadvantage for small X that negative values may occur with finite probability. Therefore, as shown in [25], a better choice is the gamma distribution

f [k, θ ] = 1 Γ(k)θ k X k-1 e -X/θ (9)
with k = X/θ . The gamma distribution is only defined for X ≥ 0, being asymmetric for small X and tending towards a gaussian distribution asymptotically for large X.

Another criterion for the choice of f is the nature of the observable X: is it a continuous or a discrete variable? Although in principle the gaussian or gamma distributions can be used for both, in the case of discrete variables such as particle multiplicities, one might think that a discrete probability distribution may be more suitable. However, the most likely candidates, the Poisson distribution, the binomial distribution and the negative binomial distribution, each impose different constraints on the reduced variance θ : θ < 1 (binomial), θ = 1 (Poisson) or θ > 1 (negative binomial). As will be seen below, the experimental P(X) distributions do not contain sufficient information to constrain the dependency of both X and θ on c b , therefore rather than constrain a priori the value of θ by the choice of f , it is preferable to use the gamma fluctuation kernel for which θ can vary freely, even for observables which only take discrete values.

Centrality dependence of the mean, X

In order to reconstruct centrality distributions from experimental X distributions, the relationship between X and centrality must be monotonic. We have chosen a functional form whose parameters can vary freely while guaranteeing monotonicity and which is sufficiently general to describe the typical shapes of such relationships as predicted by various transport model calculations in this energy range (see for example [20,21,[32][33][34][35]).

In the framework of the gamma distribution of Eq. 9, it is in fact the parameter k whose centrality dependence has to be parameterized, for which we have used

k(c b ) = k max [1 -c α b ] γ + k min ( 10
)
Note that this form corresponds to a monotonically decreasing function of centrality (as expected for variables which increase with collision violence): if a monotonically increasing function of centrality is required, one can replace c b by 1c b in Eq. 10. Only 4 free parameters are required. α and γ can be directly linked to the shape of the impact parameter dependence of X, making interpretation of fit results more immediate: the value of α determines whether or not the observable's evolution with b presents a plateau for the most central collisions, i.e. when α ≥ 1 there exists a range of small impact parameters for which the derivative dk/db ≈ 0, which implies a lower limit to the observable's sensitivity to variations of bthe larger the value of α, the larger the range; the γ-parameter determines the concavity of the curve -values of γ > 1 lead to S-shaped curves with an asymptotically zero derivative at c b = 1. k max and k min determine the maximum mean value of the observable achieved in head-on collisions, i.e. for c b = b = 0:

X max = X(b = 0) = θ (k max + k min ) (11)
The 'offset' parameter k min is important because we cannot always assume that X min = θ k min is zero for the most peripheral recorded collisions. This is especially clear when considering X = N C , the total multiplicity of charged products.

All the data analyzed in this article were obtained with an online DAQ trigger corresponding to a minimum number of fired telescopes of between 3 and 5 depending on the system studied. For the lightest systems considered in our study, the maximum charged particle multiplicity can be as small as 20; in this case the role of k min is far from negligible.

Centrality dependence of reduced fluctuations, θ

As first pointed out in [24], the experimental P(X) distributions do not contain sufficient information to constrain the dependency of both X and θ on c b . The problem is underconstrained, in the sense that one cannot extract two unknown functions, X(c b ) and θ (c b ), from a single distribution, P(X). Nevertheless, the fluctuations for central collisions, θ (b = 0), can be well-constrained by the tail of the P(X) distribution, which is dominated by fluctuations around the mean value for b ∼ 0, X(b = 0). The rest of the distribution gets contributions from many different impact parameters so that fluctuations are averaged over, and P(X) only contains information about X(c b ) away from the tail. Therefore, as in [24,25], we make the approximation that the reduced fluctuation parameter θ is constant for all centralities, and equal to its value for b = 0: in other words the variance of X is proportional to its mean value for all b, with a constant of proportionality θ which is a free parameter of the fits, constrained by the tail of the P(X) distributions.

C. Test with a transport model calculation

We have used the microscopic transport model AMD (Antisymmetrised Molecular Dynamics, [8,11]) in order to simulate a typical reaction from the INDRA dataset, 58 Ni+ 58 Ni collisions at 52 MeV/nucleon. The transport model calculation was stopped at the end of the dynamical phase of the reaction, after 300 fm/c, and the resulting set of excited primary nuclei was used as input to the statistical decay code GEM-INI++ [27,28]. The detection of the final products of these collisions by the INDRA array [START_REF] Pouthas | [END_REF] was then simulated using a software replica of the apparatus [START_REF]KaliVeda heavy-ion analysis toolkit[END_REF], including a minimumbias acquisition trigger requiring simultaneous detection of 4 charged products or more (as during the experimental measurement). Our goal here is not to test the AMD model, nor to compare its predictions with data, but rather to have at our disposition a set of simulated data which correctly reproduce at least the main features of the experimental data, while linking the impact parameter of each collision with its outcome based on the pertinent microscopic physics ingredients for heavy-ion collisions in this energy range.

A first result of these calculations is shown in Figure 1(a). The reduced variance of the two observables which we will use in the following, the total charged particle multiplicity, N C , and the total transverse energy of light charged particles with Z ≤ 2, E t12 , is presented as a function of impact parameter. The variances presented here include several sources of fluctuations: fluctuations in the microscopic collision dy- /db [mb fm

σ d 0.1 ≤ X c ≤ AMD 0.0 0.1) ≤ X c ≤ P(b|0.0 0.4 ≤ X c ≤ AMD 0.3 0.4) ≤ X c ≤ P(b|0.3 0.7 ≤ X c ≤ AMD 0.6 0.7) ≤ X c ≤ P(b|0.6 0.9 ≤ X c ≤ AMD 0.8 0.9) ≤ X c ≤ P(b|0.8 (b)
Figure 2. Impact parameter distribution reconstruction for the observable E t12 for simulated data from 58 Ni+ 58 Ni collisions at 52 MeV/nucleon (see text for details). (a) Symbols: inclusive probability distribution P(E t12 ) of total transverse energy of light charged particles. Statistical uncertainties are represented by vertical bars when larger than the symbols. Dashed curve: fit using Eq. 3 and the parametrizations of Eqs. 9,10. (b) Reconstructed impact parameter distributions (curves) calculated using Eqs. 6, 8, compared to the actual distributions (symbols) for 4 different bins of experimental centrality, c X . namics, fluctuations in the secondary decay of excited primary fragments, and fluctuations due to the finite acceptance of the INDRA array and other detection effects. It can be seen that for all but the most peripheral collisions (b < 8 fm), the reduced variance of both observables can be considered approximately constant. Therefore, the approximation of II B 3, whereby the reduced fluctuation parameter θ is assumed to be independent of the impact parameter, can be seen to have some justification, at least as far as this modelization is concerned.

Figure 2(a) shows an example of the quality of fits to the inclusive data histograms which can be achieved using Eq. 3 and the parametrizations of Eqs. 9,10, here for the total transverse energy of light charged particles, E t12 . The parameter values deduced from this fit completely determine the conditional probability distribution P(E t12 |c b ), and can therefore be used to calculate centrality distributions for any arbitrary selection of data. However, to calculate impact parameter distributions, as pointed out at the end of Sec. II A, requires the knowledge of, or at least a good approximation to, the full impact parameter distribution P(b). In the present case, of course, this distribution is accessible, as shown in Figure 1(b). Initially, the AMD calculations were run with a geometric impact parameter distribution,

P(b) = 2πb
for all b ≤ 12 fm. Let us recall that the equivalent of the experimental DAQ trigger in this simulation is the detection of at least 4 charged particles in coincidence: as can be seen, for the most peripheral collisions (b > 9 fm) the probability of such events decreases to zero. Such a distribution can be well-fitted by

P(b) = 2πbP R (b) with P R (b) = 1 1 + exp b-b 0 ∆b (12)
as shown by the dashed curve in Figure 1 [36]), and suggest that in general for "minimum-bias" INDRA data, Eq. 12 with ∆b ≈ 0.4 fm gives a good approximation to the full impact parameter distribution. The formulae allowing to calculate b 0 for a given value of ∆b knowing the total measured cross-section, and the relation between centrality and impact parameter, c b (b), needed for Eq. 8, are given in Appendix A.

Examples of reconstructed impact parameter distributions calculated using the deduced form of P(E t12 |c b ) are shown in Figure 2(b). Here we have gated the simulated data according to 4 bins of experimental centrality, c X , calculated with the E t12 variable (see Eq. 5). As can be seen, the reconstructed distributions are very similar to the actual impact parameter distributions which in this case can be simply generated by applying the E t12 cuts directly to the model. Both mean impact parameter and width of the distribution for each centrality bin are very well reproduced. For the most central bin (c X ≤ 0.1, corresponding to the upper 10% of the E t12 distribution), although a small shift can be observed, leading to a slight underestimation of the mean impact parameter for these events, the reconstructed distribution is still a very good approximation to the actual one.

If this new method underestimates slightly the impact parameters for a very central cut, it is as nothing compared to the underestimation that would result from the commonlyused sharp cut-off approximation (SCA). Within the SCA, the cross-section associated with any selection of "the most central collisions" is assumed to correspond to a geometric (triangular) distribution of impact parameters between b = 0 and b = b cut , where b cut is the assumed upper limit of impact parameters retained by the cut. In this case, the fraction of all detected events retained by the cut, i.e. the upper limit of the experimental centrality c X , is simply related to the upper limit of the impact parameter by c X = (b cut /b max ) 2 . b max is the SCA upper limit of impact parameter for all detected events, which in the present case is to a very good approximation equivalent to the value of b 0 in Eq. 12. Thus the SCA predicts for a centrality cut c X ≤ 0.1 an upper limit of impact parameter b cut ≈ 3 fm, whereas in reality the impact parameters of this selection can be seen in Figure 2(b) to have a very broad distribution which extends up to 8 fm.

III. RESULTS

In the following we will present the results of applying the method presented in Sec. II to data for a wide range of different colliding systems measured with INDRA, which are summarized in Table I. The data concern the two observables which are most commonly used for centrality estimation and/or selections in this energy range, namely the total multiplicity of charged reaction products, N C , and the total transverse energy of light charged particles (LCP, isotopes of Z = 1, 2 nuclei), E t12 . N C is the impact parameter filter most commonly-used by many different groups in the intermediate energy range, while E t12 has been especially used by the INDRA collaboration as it exploits the very high, angle-and centrality-independent efficiency of the array for detection of LCP.

A. Experimental details

Impact parameter estimation and sorting of HIC in the energy range 20-100 MeV/nucleon requires powerful multidetector arrays with high granularity and 4π angular coverage. Let us briefly recall here that INDRA [START_REF] Pouthas | [END_REF]45] is one of a second generation of 4π charged particle arrays, in continued use for the study of HIC at GANIL in Caen (and briefly at GSI, Darmstadt) since 1993. Its 336 multi-layer detection modules covering 90% of the solid angle around the target, low detection and identification thresholds, and minimum-bias trigger logic based on the number of fired modules make it ideally suited for impact parameter reconstruction in this energy range. References to the original papers where details of the data-taking can be found for each measured reaction are given in Table I.

Also in Table I, as well as the mass asymmetry, projectile energy and number of recorded events, we give also the trigger multiplicity (corresponding to the minimum number of fired modules required to record an event, which may include γ-ray, electron, pion or neutron detection in the CsI scintillators) for each reaction. In the offline analysis the same condition was applied to the reconstructed events (corresponding to a minimum number of correctly identified charged products, thus excluding γ-rays etc.).

B. Results of fits to data

Examples of fits to the inclusive distributions of the observables N C and E t12 are presented in Figure 3, for the M ≥ 4 129 Xe+ nat Sn data. Fits were performed using Eqs. 3, 9 and 10 and fit parameters are given in Tables II and III for N C and E t12 , respectively. Using the measured reaction cross-sections for this data [32], the experimental and fitted P(X) distributions are presented here as differential cross-sections. To better appreciate the quality of the fits, for both low and high statistics regions of the distributions, each is presented with both linear (top row) and logarithmic (bottom row) y-axes.

Apart from the lowest N C or E t12 values, close to the trigger threshold, the shapes of the experimental distributions are globally well reproduced by each fit, especially the exponential tails for the highest multiplicities/energies. Reduced χ 2 values for each fit are reported in Tables II andIII. For E t12 this goodness-of-fit parameter is generally excellent (χ 2 ∼ 1), whereas for N C the values are less satisfactory, despite the visual impression of adequate fits: this is due to the sharp decrease of the fitted distributions at small N C , which due to the high statistics in this region dominates the overall χ 2 values. Nevertheless, it can be remarked that the deduced X min values for N C follow remarkably well the minimum multiplicity imposed by the trigger, including the increase from M ≥ 3 to M ≥ 4 for the 36 Ar+ 58 Ni data at 74 MeV/nucleon (see Table II). Fits of similar quality for both observables were obtained for all data in this study, the parameters of which are given in Tables II andIII.

C. Deduced probability distributions for impact parameter and observable

Fitting the experimental P(X) distribution determines the parameters of the conditional probability distribution P(X|c b ) or P(X|b), and hence the joint probability distribution of observable and impact parameter, P(X, b) = P(X|b)P(b). An example is shown in Figure 4(a) for the case of 129 Xe+ nat Sn collisions at 39 MeV/nucleon using the total transverse energy of light charged particles, E t12 . As shown in the figure, the joint probability distribution encodes not only the dependence of the mean value E t12 on impact parameter, given by X = θ k(c b ) (see Eqs. 9, 10), but crucially also the fluctuations around this mean value at fixed b, determined by the parameter θ . It is clear from the joint probability distribution P(E t12 , b) shown in Figure 4(a) that any selection of events using E t12 cuts will necessarily lead to broad impact parameter distributions. This is confirmed by Figure 4(b) which presents impact parameter distributions corresponding to three E t12 bins, calculated according to Eqs. 6 and 8. From the figure it can be seen that these P(b|x 1 ≤ E t12 ≤ x 2 ) are simply projections of P(E t12 , b) on to the impact parameter axis, integrated over different ranges of E t12 . The distributions of b for the three cuts are very broad and overlap to a large extent.

Finally, Figure 4(c) shows a specific projection of P(E t12 , b) on to the E t12 axis, corresponding to P(E t12 |b = 0), the deduced distribution of E t12 for head-on collisions. It is also, as expected, very broad: the variance of this distribution is determined by the product of the reduced fluctuation parameter, θ , and the mean value E t12 (b = 0) of the E t12 observable for b = 0 collisions. The values of θ for this system and all others, and for both N C and E t12 observables are given in Tables II andIII. For the fits to 129 Xe+ nat Sn data shown in Figure 3, they are much smaller for N C (θ ∼ 0.3) than for E t12 (11 ≤ θ ≤ 16). In fact, as a general result for the whole set of reactions studied here, we find that total charged multiplicity fluctuations are relatively small (θ < 1) and nearly independent of bombarding energy, whereas E t12 fluctuations are relatively large (θ 1) and increase with beam energy. Let us recall that a similar difference in order of magnitude for θ was also found with the AMD calculations shown in Figure 1(a).

Figures 5 and6 present the evolution of N C and E t12 extrapolated to head-on collisions for all studied reactions as a function of the available center of mass energy, in order to allow equivalent comparison for symmetric and asymmetric colliding systems.

The total multiplicities in Figure 5 have been normalized to the size of each system, Z tot , given by the sum of projectile and target atomic numbers. Z tot is the maximum possible charged product multiplicity for each reaction, if complete disintegration into hydrogen isotopes and (undetected) neutrons were to occur. Thus the quantity N C (b = 0)/Z tot reflects the degree to which complete disintegration would occur in head-on collisions for each reaction. For all systems this ratio increases in a similar fashion with available energy, closely following a parabolic dependence (shown by the dashed curves in Figure 5). It should be noted that the 129 Xe+ x Sn data, taken in two separate experimental campaigns, at GANIL ( 129 Xe+ nat Sn: E proj ≤50 MeV/nucleon) and at GSI ( 129 Xe+ 124 Sn: E proj ≥65 MeV/nucleon), follows a continuous evolution over the full energy range. There is also a systematic system size/mass dependence of the ratio, which is larger at a given available energy for smaller systems. Thus for E c.m. /A tot ∼20 MeV/nucleon it can be seen that the 58 Ni+ 58 Ni system reaches 50% of complete disintegration, whereas for 197 Au+ 197 Au the degree is no greater than 35%. We remark that the asymmetric system, 58 Ni+ 197 Au, does not appear to follow the same systematic trend: it could have been expected to have the same dependence on E c.m. /A tot as 129 Xe+ x Sn, which has approximately the same total mass and charge but a quasi-symmetric entrance channel.

Figure 6 shows the extrapolated E t12 for b = 0 collisions for all systems, again as a function of available energy. The normalization E t12 (b = 0)/Z tot to the total system charge here has no direct physical interpretation, but empirically we have found that it leads to a near-universal scaling of the data which collapse onto a linear dependence of the total transverse energy of LCP that would be achieved in head-on collisions as a function of the c.m. energy. Again, it can be remarked that data from campaigns in both GANIL and GSI ( 129 Xe+ 124 Sn: E proj ≥65 MeV/nucleon and 197 Au+ 197 Au) follow the same systematic trends, including this time the 58 Ni+ 197 Au system; on the other hand, here it is the 36 Ar+KCl data which do not follow the trend, in actual fact because the unscaled E t12 values for this system are almost identical to those for 36 Ar+ 58 Ni.

We know of no reason why the mean transverse energy of light charged particles E t12 for b = 0 collisions should scale in this way, over such a wide range of system sizes, asymmetries, and bombarding energies. It is our opinion that this result should be compared with the predictions of different transport models in this energy range, for which it could provide a benchmark test. E. System mass-and bombarding-energy dependence of mean impact parameters of events selected with a high-E t12 cut

As we recalled in the introduction, it is a well-known fact that when trying to select more and more central collisions using high-N C or E t12 cuts, a large amount of impact parameter mixing occurs leading ultimately to a finite limit on the true centrality of any selected sample of events. We can now quantify this saturation by providing an estimate of b for any selection of data, without relying on any specific model of collisions. We present here the mean values of reduced impact parameter for cuts corresponding to experimental centralities c X < 10% or c X < 1% performed using the total transverse energy of LCP, E t12 . It should be noted that equivalent results are found using the N C observable.

Figure 7 shows the mean value of b/b max for the impact parameter distributions we have reconstructed for all systems and bombarding energies and for the two cuts, as a function of the available center of mass energy of each reaction. The values b/b max can be seen to cluster around a near-universal energy dependency for each cut, which appears to be independent of system mass and entrance channel asymmetry. In each case, b/b max decreases with increasing center of mass energy. With the more restrictive cut of c X < 1%, the selectivity of the cuts can be seen to increase significantly, especially for the highest energies. Values of 0.35 ≥ b/b max ≥ 0.24 are obtained for the 10% centrality cut, and 0.29 ≥ b/b max ≥ 0.14 with the 1% centrality cut.

It is instructive to compare these results with the predictions of the sharp cut-off approximation (SCA) for the same cuts (see the end of Sec. II C). The SCA would expect an ) for all datasets. See Sec. II B for meaning of parameters. χ 2 is the reduced chi-square value of each fit.

IV. CONCLUSIONS

One way to improve constraints on the nuclear equation of state from comparisons between data on intermediate energy heavy-ion collisions and transport model calculations is by providing a model-independent estimation of the impact parameter distribution representative of any selected set of experimental data. To do so requires to explicitly take into account the fluctuations in the relationship between any observable X and the impact parameter b, as first shown in [24,25]. In this article we have shown how the method, first developed for ultra-relativistic collisions, can be adapted and used in the 20-100 MeV/nucleon bombarding energy range. Notably, we have proposed a new parametrization of the relationship between the mean value of an observable and the impact parameter whose parameters are simple to interpret in terms of the shape of this relationship.

We have shown, using a complete simulation of 58 Ni+ 58 Ni collisions at 52 MeV/nucleon measured by the INDRA array, calculated using the AMD transport model coupled with the GEMINI++ statistical decay code, that the method is capable of reconstructing the impact parameter distributions associated to experimental events taking account of secondary decay and finite detector acceptance effects. In this example, the impact parameter distributions for a set of experimental cuts defined over the full range of centralities were correctly reproduced, and even for the most central cut for which the method shows a slight underestimation of the actual impact parameters, the estimated impact parameter distribution is far closer to the truth than an estimation based on the commonly-used sharp cut-off approximation. These calculations also indicate that one of the major assumptions of the method, that the reduced variance of the observable is independent of centrality, is a reasonable approximation, at least for the two observables we have studied in this paper.

We then applied the approach to a very wide range of data for different collisions measured with INDRA between 25 MeV/nucleon and 100 MeV/nucleon, where in each case the data were recorded according to a "minimum-bias" trigger based on a minimum number of fired telescopes over nearly the full 4π solid angle around the target. Two commonlyused observables have been employed, the total multiplicity of charged products, N C , and the total transverse energy of light charged particles with Z ≤ 2, E t12 . Excellent fits to the inclusive distributions P(X) of each observable have been achieved for all collisions, even if for N C the failure to reproduce the distributions close to the minimum bias trigger (most peripheral collisions) biases the apparent goodness-of-fit as measured by reduced χ 2 values.

The parameters determined by the fits allow to deduce the joint probability distribution P(X, c b ) from which impact parameter distributions for any selection of data can be reconstructed, or distributions of the observable for a given range of b. The relative fluctuation of the joint probability distributions about the mean value of the observable for each b has been shown to differ according to the observable, with the total multiplicity of charged products, N C , associated with sub-Poissonian fluctuations (θ < 1), whereas the total transverse energy of light charged particles, E t12 , exhibits much larger fluctuations (θ 1), for all data studied in this article. The asymptotic values of the mean values of both observables for b = 0 collisions can be extrapolated from the fit results, and may provide new constraints for transport model calculations. We have shown, in particular, that the total transverse energy of light charged particles has mean values for head-on collisions which show a near-universal dependence on the available energy in the center of mass of the collisions. This result should be confronted with different microscopic model predictions.

Finally, we have characterized the true centrality of a commonly-used event selection employing high-E t12 cuts to retain the "most central collisions". The results are largely independent of total system mass and mass-asymmetry of the entrance channel, each system showing very similar evolution of the mean reduced impact parameter b/b max as a function of available energy for each centrality cut. The actual representative mean values of reduced impact parameters for these selections were shown to decrease with increasing bombarding energy from 0.35 to 0.24 (for a 10% centrality cut), or from 0.29 to 0.14 (for a 1% centrality cut), when the usual sharp cut-off approximation (SCA) gives mean values of 0.2 or 0.07, respectively, for these two centrality cuts, independently of the reaction bombarding energy. This overestimation of the centrality of each data sample would skew comparison with any transport model by using the wrong impact parameters as input. This is why we have tried to demonstrate in this paper that a new, model-independent method for estimating the impact parameter distributions of selected experimental events is feasible and should be used whenever possible in order to improve the constraints that can be brought on the description of nuclear dynamics and the nuclear equation of state by comparisons between experimental data and microscopic transport model calculations.

This definite integral is related to the complete Fermi-Dirac integral F j (x) = 1 Γ( j + 1) ˆ∞ 0 t j 1 + exp (tx) dt (A1) with j = 1, where Γ( j + 1) is the gamma function, Γ( j + 1) = j! for integer j. In general the value of this integral is given by a polylogarithm, Li s (z), specifically

F j (x) = -Li j+1 (-e x )
and in this particular case by the negative dilogarithm, -Li 2 (-e x ). Therefore we have for the final expression of the total cross-section which normalizes correctly the probability distribution of Eq. ( 12),

σ R = -2π(∆b) 2 Li 2 -exp b 0 ∆b (A2)
This expression can be used to find b 0 for a given total crosssection and width parameter ∆b, by numerical inversion [46].

Analytic expression for centrality

To calculate the centrality c b (b) we substitute Eq. ( 12) into Eq. ( 2 This definite integral can be calculated using the incomplete Fermi-Dirac integral F j (a, x) = 1 Γ ( j + 1) ˆ∞ a t j 1 + exp (tx) dt, a ≥ 0 with a = b/∆b, and the complete Fermi-Dirac integral F j (x) of Eq. (A1):

ˆa 0 t j 1 + exp (t -x)
dt = Γ ( j + 1) [F j (x) -F j (a, x)]

With j = 1, F 1 (x) = -Li 2 (-e x ) as above, while the incomplete FD integral can be written (by integration by parts) as

F 1 (a, x) = π 2 6 - (a 2 -x 2 )
2 + a ln 1 + e (a-x) + Li 2 -e (a-x)

The final expression for the centrality is therefore

Figure 1 .

 1 Figure 1. AMD+GEMINI++ simulations of 58 Ni+ 58 Ni collisions at 52 MeV/nucleon filtered according to the acceptance of the INDRA array (see text for details). (a) Reduced variances of the observables N C and E t12 as a function of impact parameter, b. (b) Symbols: differential cross-section distribution of all detected simulated events as a function of b. Dashed curve: fit to the distribution using Eq. 12 with b 0 = 10.1 fm and ∆b = 0.4 fm.

  (b), for parameter values b 0 = 10.1 fm and ∆b = 0.4 fm. Similar distributions have been observed previously using different models applied to different reaction systems of the INDRA dataset (see for example Figure 2(a) of

Figure 3 .

 3 Figure 3. Results of fits to the inclusive distributions of N C (left column, (a) & (c)) and E t12 (right column, (b) & (d)) for the 129 Xe+ nat Sn data. Each distribution is presented with both linear (top row, (a) & (b)) and logarithmic (bottom row, (c) & (d)) y-axes. Statistical uncertainties on the data are shown when not smaller than the symbols.

Figure 4 .

 4 Figure 4. Results obtained from the fit to the experimental E t12 distribution for 129 Xe+ nat Sn collisions at 39 MeV/nucleon: (a) reconstructed joint probability distribution P(E t12 , b) (logarithmic color scale) and the mean value E t12 (b) (black curve); (b) impact parameter distributions for three cuts in E t12 indicated in the legend; (c) distribution of E t12 for b = 0.

Figure 5 .

 5 Figure 5. Extrapolated mean values of total charged product multiplicity for b = 0 collisions deduced fits, N C (b = 0), normalized to total system charge Z tot , as a function of available center of mass energy per nucleon E c.m. /A tot for all systems. Dashed curves are parabolic fits.

Figure 6 .

 6 Figure 6. Extrapolated mean values of total transverse energy of light charged particles for b = 0 collisions deduced from fits, E t12 (b = 0), normalized to total system charge Z tot , as a function of available center of mass energy per nucleon E c.m. /A tot for all systems. The dashed gray line is a linear fit to all the data except36 Ar+KCl.

Figure 7 .

 7 Figure 7. Mean reduced impact parameter b/b max for central collisions selected with a 10% (open symbols) or 1% (full symbols) centrality cut using E t12 , as a function of available center of mass energy per nucleon, E c.m. /A tot . Dashed curves are to guide the eye.

1 +

 1 ), and making the same substitutions as above (b = t∆b, b 0 = x∆b ) we find c b (b) = 2πexp (tx) dt

( 1 +

 1 exp ((bb 0 )/∆b)) -Li 2 -e (b-b 0 )/∆b (A4)

Table I

 I

	System	Mass	E proj	Trigger	Events
		asymmetry [MeV/ multiplicity
			nucleon]
			31.54	3	3216332
	36 Ar+KCl [37]	0.00	39.97 51.66	3 3	3496188 2391311
			74.00	3	3337570
			31.54	3	8259867
			39.97	3	7234383
			51.66	3	8599855
	36 Ar+ 58 Ni [38, 39]	0.23	63.03	3	5020363
			74.00	4	7648474
			83.63	4	4657028
			95.22	4	9799670
			31.98	4	4538513
			52.00	4	4738429
	58 Ni+ 58 Ni [40, 41]	0.00	63.63 73.96	4 4	4473639 5198692
			82.00	4	5578566
			90.00	4	9144521
			31.98	4	7448285
			52.00	4	7941858
	58 Ni+ 197 Au [42]	0.55	63.63	4	4720169
			73.96	4	6685519
			82.00	4	7398023
			24.98	4	5288164
			32.00	4	3916797
	129 Xe+ nat Sn[32, 43]	0.04	38.98	4	5261377
			45.00	4	6067739
			50.13	4	5792220
			65.00	3	881642
	129 Xe+ 124 Sn(*) [31]	0.02	80.00	3	424357
			100.00	3	1328486
			40.00	3	2783629
	197 Au+ 197 Au(*) [44]	0.00	60.00 80.00	3 3	7589902 3545170
			100.00	3	10691556

. Characteristics of collisions studied in this work: mass asymmetry |A p -A t |/(A p + A t ), beam energy, DAQ trigger multiplicity and total number of recorded events. References are given to the original papers where details of the data-taking can be found. Systems marked with an asterisk were measured at GSI, all others at GANIL.

Table II .

 II Parameters of fits to total charged particle multiplicity distributions P(N C ) for all datasets. See Sec. II B for meaning of parameters. χ 2 is the reduced chi-square value of each fit. energy-independent upper limit of b cut /b max ≈ 0.3 for the 10% cut, or b cut /b max = 0.1 with a 1% centrality cut. For a geometric (triangular) impact parameter distribution, the mean value is simply related to the cut-off by b SCA = 2 /3b cut , therefore the mean values in Figure7should be compared to b SCA /b max ≈ 0.1 in the case of the 10% centrality cut, and b SCA /b max ≈ 0.07 for c X < 1%. It is clear that impact parameter estimation using the SCA for central collisions greatly underestimates representative values of b for the selected data samples. The method presented in this paper can characterize data not only in terms of far more representative mean values of b, but indeed provide an estimation of the actual impact parameter distribution which could be used as input to a transport model calculation.

	System	E proj	α	γ	θ X max X min χ 2	System	E proj	α	γ	θ	X max X min χ 2
		[MeV/nucleon]						[MeV/nucleon]			[MeV] [MeV] [MeV]
		32	0.95 1.25 0.20 14.1 3.7 13.2		32	0.35 0.76 6.1	162	3 1.2
	36 Ar+KCl	40 52	0.98 1.46 0.23 15.6 4.2 0.88 1.40 0.22 17.6 4.0	6.5 8.2	36 Ar+KCl	40 52	0.37 0.89 7.5 0.35 1.02 8.5	196 269	8 1.0 12 1.2
		74	0.89 1.60 0.21 19.8 4.3 10.2		74	0.40 1.32 11.8	389	19 3.7
		32	1.37 1.12 0.25 14.8 2.5 85.4		32	0.97 1.17 8.5	148	9 3.3
		40	1.23 1.14 0.27 16.8 2.6 46.4		40	0.83 1.17 10.0	183	11 2.5
		52	1.07 1.17 0.28 19.1 2.7 53.5		52	0.68 1.26 12.0	251	15 2.3
	36 Ar+ 58 Ni	63	0.99 1.19 0.29 20.6 2.7 28.1	36 Ar+ 58 Ni	63	0.60 1.35 13.1	328	18 1.3
		74	0.95 1.21 0.27 22.4 3.8 43.8		74	0.60 1.46 14.9	402	30 1.8
		84	0.96 1.24 0.27 23.2 3.8 25.8		84	0.60 1.52 16.1	463	32 1.7
		95	0.98 1.28 0.27 24.1 3.9 56.3		95	0.62 1.63 18.6	528	35 2.1
		32	1.24 1.19 0.28 16.3 3.7 55.4		32	0.79 1.04 9.9	186	21 1.6
		52	0.97 1.15 0.30 21.6 3.8 14.0		52	0.56 1.15 13.0	340	29 2.3
	58 Ni+ 58 Ni	63 74	0.92 1.16 0.28 24.1 3.8 0.81 1.17 0.27 25.8 3.7	5.9 4.2	58 Ni+ 58 Ni	64 74	0.55 1.30 14.9 0.52 1.40 16.7	443 541	33 1.7 37 2.8
		82	0.96 1.45 0.29 26.4 4.6	8.1		82	0.61 1.68 18.5	591	46 1.7
		90	0.93 1.40 0.28 27.4 4.2 38.2		90	0.68 1.93 20.9	629	56 1.4
		32	1.79 1.79 0.34 16.2 4.1 84.3		32	1.41 1.71 12.8	223	32 6.7
		52	1.63 1.66 0.37 23.0 4.2 49.5		52	1.08 1.45 17.2	391	36 6.1
	58 Ni+ 197 Au	64	1.48 1.64 0.37 26.4 4.3 15.1	58 Ni+ 197 Au	64	0.93 1.50 18.8	521	39 3.8
		74	1.44 1.62 0.39 28.4 4.1 45.6		74	0.92 1.63 21.6	630	44 4.8
		82	1.51 1.85 0.41 29.7 4.7 43.2		82	0.92 1.74 23.3	716	48 5.8
		25	1.26 0.95 0.32 19.6 2.8 63.5		25	0.74 0.68 11.1	241	6 2.2
		32	1.24 1.08 0.34 24.0 3.2 40.3		32	0.67 0.69 12.4	310	5 1.7
	129 Xe+ nat Sn	39	1.18 1.17 0.34 28.4 3.5 55.1	129 Xe+ nat Sn	39	0.57 0.75 13.9	408	7 2.9
		45	1.14 1.23 0.34 31.2 3.8 56.5		45	0.55 0.89 15.4	496	24 1.6
		50	1.14 1.35 0.34 34.1 4.0 36.6		50	0.57 1.06 16.0	584	34 1.3
		65	1.09 1.40 0.36 38.4 2.9	2.0		65	0.59 1.32 19.9	822	34 1.5
	129 Xe+ 124 Sn	80	1.11 1.50 0.34 42.5 3.1	1.5	129 Xe+ 124 Sn	80	0.61 1.52 23.5	1071	41 2.0
		100	1.18 1.64 0.38 45.2 3.5	2.7		100	0.62 1.65 26.8	1374	44 3.6
		40	1.23 1.27 0.42 35.5 2.1 30.1		40	1.07 1.23 24.0	521	18 24.4
	197 Au+ 197 Au	60 80	1.22 1.62 0.45 47.8 1.2 68.1 1.24 1.61 0.45 54.9 2.8 16.9	197 Au+ 197 Au	60 80	0.67 1.38 25.5 0.62 1.47 26.6	1089 1648	5 24.6 31 3.0
		100	1.26 1.64 0.49 58.5 3.1 52.0		100	0.65 1.68 34.6	2054	55 4.4
							Table III. Results of fits to total transverse LCP energy distributions
							P(E t12

ACKNOWLEDGMENTS

We would like to thank all the technical staff of GANIL for their continued support in performing the experiments. We gratefully acknowledge support from the CNRS/IN2P3 Computing Center (Lyon -France) for providing computing and data-processing resources needed for this work. We would also like to thank the ROOT [46] development team without whose software none of the analyses would be possible. We acknowledge support from Région Normandie under the Réseau d'Intérêt Normand FIDNEOS (RIN/FIDNEOS). We would also like to signal that the software necessary to perform the analyses presented in this paper will be made available to the entire community as part of the KaliVeda heavy-ion analysis toolkit [29] upon publication.

Appendix A: Inclusive impact parameter distributions for INDRA data

In order to transform deduced centrality distributions P(c b |S) into impact parameter distributions using Eq. 8 requires to calculate the centrality for each impact parameter, c b (b), (Eq. A3) and deduce the value of b 0 from the (measured) total reaction cross-section by numerical inversion of Eq. A2, assuming a typical value of ∆b ≈ 0.4 fm.

Analytic expression for total cross-section

To calculate the total reaction cross-section for a given set of parameters b 0 and ∆b, we have, by definition,


H Lebbad

B Appolaire

Y Le Bouar

A Finel

Insights into the selection mechanism of Widmanstätten growth by phase-field calculations

Keywords: phase transformation, phase field, elasticity, Widmanstätten

We investigate the growth of Widmanstätten structures in metallic alloys with phase field calculations. We show that elastic energy anisotropy associated with shear dominated transformations is alone sufficient to give rise to plate or needle like precipitates with constant lengthening rates at constant supersaturations. We show that Ivantsov analytical solutions still provide good estimates of the diffusion flux around tips even when elasticity is accounted for. A careful analysis demonstrates that, for shear dominated transformations, (i) the lengthening direction lies in the habit plane determined by one of the minima of the elastic kernel; (ii) the tip shape and size are equilibrium features, driven by elastic forces, rather than dynamically determined as usually assumed; (iii) the tip size of acicular precipitates can be rationalized by the value of the elastic kernel in the lengthening direction.

Introduction

Since they have been discovered at the beginning of the 19 th century in meteoritic Fe-Ni [START_REF] Thomson | On the malleable iron -essai sur le fer malléable trouvé en sibérie par le Prof. Pallas[END_REF][START_REF] Osmond | Microscopic Analysis of Metals[END_REF][START_REF] Mehl | Studies upon the Widmanstätten structure, i. -introduction. the aluminium-silver system and the copper-silicon system[END_REF][START_REF] Belaiew | The Widmanstätten structure in various alloys and metals[END_REF], Widmanstätten structures have been observed in many metallic alloys, such as Fe-C [START_REF] Kinsman | Thickening kinetics of proeutectoid ferrite plates in Fe-C alloys[END_REF], Al-Ag [START_REF] Laird | The growth of γ-plates in an Al-15% Ag alloy[END_REF], Ag-Cd [START_REF] Kostic | Growth kinetics of bainite plates and Widmanstätten needles in the β phase of a Ag-45at.% cd alloy[END_REF], Cu-Zn [START_REF] Purdy | Widmanstätten Precipitation from non-ideal solid solution: α in β-CuZn[END_REF][START_REF] Stephens | On Widmanstätten growth from beta brass[END_REF] and Ti alloys [START_REF] Menon | Interfacial structure of widmanstätten plates in a ti-cr alloy[END_REF][START_REF] Appolaire | Modelling of phase transformation kinetics in ti alloys -isothermal treatments[END_REF], among others. These microstructures, ensuing from diffusioncontrolled phase transformations, display acicular morphologies (plates or needles), often growing from grain boundaries of the parent phase in colonies of precipitates with parallel habit planes. They are observed at intermediate temperatures or after coolings at moderate rates. In isothermal conditions, the precipitates lengthen at constant rate.

Understanding the formation of these microstructures has long been relying on either the prediction of the lenghtening kinetics or the description of the interface structure. Along the first line, all works start with the analysis of the diffusion problem around self-similar acicular shapes growing at constant velocity, following the pioneer works of [START_REF] Ivantsov | Temperature field around a spherical, cylindrical, and needleshaped crystal, growing in a pre-cooled melt[END_REF][START_REF] Zener | Theory of growth of spherical precipitates from solid solution[END_REF][START_REF] Horvay | Dendritic and spheroidal growth[END_REF][START_REF] Trivedi | Growth kinetics of plate-like precipitates[END_REF][START_REF] Hillert | Diffusion-controlled lengthening of widmanstätten plates[END_REF]. However, the resulting models have faced two major issues. First, the diffusion problem is not closed and must be supplied with a selection criterion to obtain unique tip velocity and radius. Second, discrepancies were noted between the predictions and the few available measurements in different alloys, e.g. [START_REF] Bhadeshia | Diffusion-controlled growth of ferrite plates in plaincarbon steels[END_REF]. Starting in the 60s, Aaronson and co-workers were convinced that the above descriptions were not sufficient because they do not account for the interface structure that should play an active role in the growth process, explaining both its kinetics and the resulting morphology. Using TEM observations and bicrystallography, they identified interfacial ledges as the main vec-tor of interface propagation. Based on these works, competing models were proposed that assume that the propagation is controlled by the diffusion of solute species at the ledge risers [START_REF] Atkinson | The growth kinetics of individual ledges during solid-solid phase transformations[END_REF][START_REF] Enomoto | Computer modeling of the growth kinetics of ledged interphase boundaries --i. single step and infinite train of steps[END_REF][START_REF] Enomoto | Thermodynamics and kinetics of the formation of widmanstätten ferrite plates in ferrous alloys[END_REF]. Unfortunately, these models did not improve the agreement with the experiments, e.g. [START_REF] Hoyt | The velocity of plate precipitates growing by the ledge mechanism[END_REF]. The discussion between these two routes did not end on some real conclusion, but rather on a standstill. In fact, it can be noted that all these works have neglected the role of elastic energy associated with the lattices mismatch, although regularly mentioned as the source of discrepancies, and sometimes partially accounted for [START_REF] Dryden | Elastic interactions of transformation ledges and precipitate edges[END_REF][START_REF] Enomoto | Computer simulation of ledge migration under elastic interaction[END_REF]. Using the capabilities of phase field models coupled with elasticity and plasticity [START_REF] Wang | Shape evolution of a precipitate during strain-induced coarsening -a computer simulation[END_REF][START_REF] Finel | Phase field methods: Microstructures, mechanical properties and complexity[END_REF][START_REF] Ammar | Phase field modeling of elasto-plastic deformation induced by diffusion controlled growth of a misfitting spherical precipitate[END_REF][START_REF] Cottura | A phase field model incorporating strain gradient viscoplasticity: Application to rafting in ni-base superalloys[END_REF][START_REF] Thuinet | Phase-field modeling of precipitate evolution dynamics in elastically inhomogeneous low-symmetry systems: Application to hydride precipitation in zr[END_REF][START_REF] De Rancourt | Homogenization of viscoplastic constitutive laws within a phase field approach[END_REF], Cottura et al. [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF] have shown recently that the general features of Widmanstätten structures can be recovered by simply accounting for the anisotropic elastic energy produced by the change in crystal structure. Indeed, they have shown that there is a direction of constant lengthening rate provided that the elastic energy displays a particular dependency with respect to the interface orientation (see also [START_REF] Appolaire | Stress and strain fields associated with the formation of α" in near-β titanium alloys[END_REF] for a discussion on its relation to the invariant line), that seems to be a feature common to the different alloys where Widmanstätten structures are observed. However, although they have clearly shown that the eigenstrain carries all the necessary information for understanding the direction of growth, the selection mechanism of the tip size has remained puzzling. This is the issue that we address numerically in the present work.

After a brief explanation of the phase field model that we have used, a careful characterization of the steady plate lengthening in different conditions will be carried out so as to be critically compared to analytical solutions and available measurements. Using this analysis, a selection criterion for the tip size and velocity will be proposed.

Phase-field model

The model used in the present work is the same as in [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF]. It relies on a phase-field parameter ϕ(x, t) that distinguishes the matrix (ϕ = 0), the precipitate (ϕ = 1) and the interface (0 < ϕ < 1). The microstructure is also described by concentration c(x, t) and displacement u(x, t) fields.

The mesoscopic free energy F of the system is split into three contributions.

F = Wd 3 F = Wd 3 V dV f ch (ϕ, c) + f int (ϕ) + f el (u, ϕ) (1)
where W and d are respectively the energy density and space scales and where all quantities in the integral (and so F) are non-dimensional. The chemical contribution accounts for the free energy approximated in each phase by a parabola with respect to c, then interpolated within the interface,

f ch (c, ϕ) = K c 2 c -h(ϕ) 2 (2)
where h(ϕ) = 3ϕ 2 -2ϕ 3 and where it is assumed, without loss of generality, that the equilibrium concentrations are 0 in the matrix and 1 in the precipitate. For simplicity, we have also considered that the curvature K c of the free energies with respect to concentration are the same for both phases.

The interface contribution reads:

f int (ϕ) = ϕ 2 (1 -ϕ) 2 + λ 2 |∇ϕ| 2 (3)
with ∇ the non-dimensional gradient operator, scaled by d.

The elastic contribution reads:

f el (u, ϕ) = 1 2 ε e i j (u, ϕ) : C i jkl : ε kl e (u, ϕ) (4)
where C is the fourth order tensor of elastic constants scaled by W, supposed to be homogeneous for simplicity. In the small strain framework, the elastic strain is given by ε

e i j = ε i j -ε * i j (ϕ) = 1 2 ∇ i u j + ∇ j u i -ε * i j (ϕ).
Considering the matrix as the reference, the eigenstrain writes

ε * i j = h(ϕ) ε T i j where ε T i j
characterizes the change of crystal structure in the precipitate with respect to the matrix.

To ensure that the free energy functional decreases with time, relaxation dynamics are assumed for both c(x, t) and ϕ(x, t). The time evolutions of these fields read:

∂ϕ ∂t = -L δF δϕ (5)
∂c ∂t = 1 K c ∇ 2 δF δc (6)
where the time scale is defined as τ = d 2 /D with D the interdiffusion coefficient, and where the phase field mobility L is a positive constant.

Because the relaxation of displacements u(x, t) is much faster than diffusion, static mechanical equilibrium is considered:

∇ j • δF δε i j = 0 (7)
All equations are expressed in Fourier space to transform space differential operators into algebraic operations. Time is discretized using a semi-implicit first order Euler scheme. Considering for simplicity homogeneous elastic constants C i jkl , Eq. ( 7) can be solved very efficiently using the Fourier transform of the normalized Green function for displacements whose inverse is given by Ω -1 i j (n) = C ikl j n k n l where n = k/k is a unit vector of reciprocal space. Moreover, we can take advantage of Fourier space to express the elastic energy as:

E = 1 2 - d 3 k (2π) 3 B(n) |θ(k)| 2 (8)
where the -is the Cauchy principal value of the integral, and where θ ϕ(x, t) describes the microstructure (θ ϕ(x, t) = 0 = 0 and θ ϕ(x, t) = 1 = 1).

B(n) = ε * i j C i jkl ε * kl -n i σ * i j Ω jk (n)σ * kl n l (9)
is the elastic kernel, with σ * i j = C i jkl ε * kl . The analysis of our results will strongly rely on the anisotropy of the elastic energy density that can be visualized simply with polar plots of this kernel.

The interface energy σ int and width δ are related to the energy density W and length scale d by

σ int = Wd √ λ/(3 √ 2) and δ = 2d √ 2λ
, where λ is the non-dimensional gradient coefficient. We have selected λ = 1.92. This choice ensures that the interface width δ is large enough with respect to the discretization length d used for the numerical implementation of the gradient operator and Fourier transforms. With the present choice, the interfaces are resolved with 5 to 6 nodes. Moreover, we have selected K c = 0.77. Typical values for the interface energy and width can be obtained by selecting appropriate values for the energy density scale W and length scale d. For example, d = 0.8 nm and W = 7.65 × 10 8 J m -3 , interface energy is σ int =200 mJ m -1 and interface width is δ = 3.2 nm. Then, apart from section 4, the elastic tensor is assumed isotropic with Lamé coefficients λ 0 = 200GPa and µ 0 = 50GPa. Finally, to ensure that the overall dynamics is controlled by diffusion, the phase field mobility L must be sufficiently larger than one. We have selected L = 6.49.

Growth dynamics of 2D plates and 3D needles

2D calculations have been performed considering 5120×512 grids. Initially, circular nuclei with radius 20d are inserted at the center of the boxes, and the matrix is supersaturated with concentrations c ∞ in the range [0.1; 0.35]. According to Cottura et al. [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF], the simplest eigenstrain giving rise to Widmanstätten structures writes, in the frame of the computation box ([001] being perpendicular to the 2D section of interest):

ε * 2D = ε 0           0 0 0 0 -1 0 0 0 0           (10)
with ε 0 the amplitude of the eigenstrain, ranging from 1% to 5%. Indeed, this eigenstrain has been shown to generate an elastic energy highly anisotropic with respect to the interface orientation: the 2D polar plot of the corresponding elastic kernel B(n) displays two opposite maxima along [010] (noted e 3 ) and two opposite minima along [100] (noted e 1 ) corresponding to the habit plane normals (Fig. 1a) such as to make a figure eight. The analytical expression of B(n) is given in Appendix A.

3D calculations have been performed considering 128×128× 512 grids. Spherical nuclei with radius 20d are inserted at the center of the boxes and initial concentrations in the matrix are set so that the supersaturations range from 0.2 to 0.35. As shown by Cottura et al. [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF], the following eigenstrain in the frame of the computation box:

ε * 3D = ε 0           0 -1 0 -1 0 0 0 0 0           (11)
is typical of real systems displaying Widmanstätten structures (e.g. steels, Ti alloys or brass). The 3D polar plot corresponding elastic kernel is shown in Fig. 1b. It features two pairs of minima along direction [100], noted e 1 , and direction [010], noted e 2 , and one pair of maxima along direction [001], noted e 3 : as for the 2D case, the minima are associated with possible habit planes. In the simple case of ε * 3D given in [START_REF] Appolaire | Modelling of phase transformation kinetics in ti alloys -isothermal treatments[END_REF], the minima and maximum directions are mutually orthogonal. It is worth noting that the (e i 3 , e 3 ) sections of the 3D polar plot are similar to the 2D polar plot of ε * 2D , such that a cylindrical nucleus with axis along either e 1 or e 2 would give plates similar to the 2D calculations. However, considering quadratic calculation boxes with equivalent sides along e 1 and e 2 together with an isotropic nucleus and periodic boundary conditions will give rise to Widmanstätten needles rather than plates.

In 2D, the nuclei grow rapidly as plates with the longest direction along e 3 (Fig. 2). In 3D, due to the particular boundary conditions, the spherical nuclei grow rapidly as needles with the longest direction also along e 3 (Fig. 3). In both cases, the lengthening rate is constant (after a short transient stage) at a given supersaturation and can easily be obtained by a linear fit of length vs time.

In 2D, we have measured the radius of curvature of the tips following the protocol described in [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF]. In 3D, we have followed [START_REF] Karma | Three-dimensional dendrite-tip morphology at low undercooling[END_REF] as illustrated in Fig. 3. In both cases, the measurement uncertainty is lower than a grid spacing d. The overall 3D shape of the tip (corresponding to ϕ = 0.5) is shown together with sections perpendicular to the lengthening direction e 3 . The projections of these sections in the (e 1 , e 2 ) plane are plotted in 4. The needles display planar sides oriented normal to e 1 and (a) e 2 and an axisymmetric tip along e 3 . The area of the precipitate section (defined by ϕ > 0.5) is plotted in Fig. 4 as a function of the distance from the very tip: the slope α of the linear fit gives the radius of curvature √ α/π of the corresponding paraboloid of revolution.

The lengthening velocity V and tip radius R are used to compute the Péclet number P = RV/(2D). P is plotted in Fig. 5 (full symbols) with respect to the tip supersaturation defined as Ω = c ∞c i where c i is the interfacial concentration on the matrix side at the tip. This concentration is obtained by extrapolating the concentration profile in the matrix (where φ < 0.05) close to the diffuse interface along the lengthening axis up to φ = 0.5. We have also plotted the predictions of Ivantsov model [START_REF] Ivantsov | Temperature field around a spherical, cylindrical, and needleshaped crystal, growing in a pre-cooled melt[END_REF] for the growth of plates and needle in infinite supersaturated matrix. For a parabolic plate (2D), P is implicitly related to the supersaturation as:

Ω = Iv 2D (P) = √ πP exp(P) erfc √ P (12)
For a paraboloidal needle (3D), the implicit relation reads: where E 1 is the exponential integral. It must be stressed that theses analytical solutions assume local equilibrium at the interface without any influence of curvature or elasticity, i.e. interfacial concentrations are provided by the phase diagram.

Ω = Iv 3D (P) = P exp(P) E 1 (P) (13)
Overall, the phase field calculations (full symbols) are close to Ivantsov solutions (continuous lines). For 2D plates, all phase field outcomes (red symbols) are slightly above the Ivantsov red curve, whatever the supersaturations and eigenstrain amplitudes. For a given eigenstrain, the scaling of P with Ω is close to Ivantsov scaling (P ∝ Ω 2.5 in the range of supersaturations investigated, slightly above the scaling P ∝ Ω 2 corresponding to the lowest order expansion of Iv 2D when P → 0). For 3D needles, the phase field outcomes (black symbols) are slightly more scattered with respect to the Ivantsov black curve than in 2D, with values above Ivantsov for 4 and 5% eigenstrains, and below for 3% eigenstrain. In this last case, the scaling of P with Ω seems to be larger than Ivantsov. These differences may have several origins.

(i) First, the periodic boundary conditions are likely to influence the tip velocity. Indeed, although the 2D periodic box has been selected by successively increasing its dimensions until converging to the same results, this process has not been fully possible in 3D due to the associated computational costs (in terms of CPU and memory). Still, the ratios between the diffusion lengths D/V to the lateral sizes of the boxes are only of the order of unity for the smallest supersaturations in both 2D and 3D (see supplemental material). Consequently, this may not be the main reason for explaining the differences with Ivantsov (ii) Second, the shapes are not strictly parabolic or paraboloidal. This is pronounced for 3D with flat sides not far from the tip (Fig. 3). For 2D, we have quantitatively compared the tip shapes with parabolae, as illustrated in Fig. 6: the tips are sharper than parabolae with equivalent tip curvatures (blue circle), with stronger deviations (i.e. earlier along the long axis) from the parabolae for larger eigenstrain amplitudes. It must be emphasized that the shapes do not depend on the growth condition, i.e. supersaturation c ∞ , as illustrated in Fig. 6: this is an important outcome of our calculations that will be elaborated further in §4.

(iii) Finally, the interface is not a concentration level set, as clearly shown in Fig. 7. This observation has been overlooked in [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF] where only concentration level sets far from the tips were plotted. Although it is difficult to prove that local equilibrium holds at the interface, it can be shown that Gibbs-Thomson together with elasticity should play the major role to set the interfacial concentrations (see supplemental material).

10 -2 10 -1 10 0 Ω 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 P ε 0 = 3% ε 0 = 4% ε 0 = 5% Fe-C Cu-Zn Al-Ag
To conclude this comparison, the predicted deviations from Ivantsov solutions can be considered as moderate in view of the significant eigenstrain investigated. Additional calculations performed with more complex eigenstrains (see §4) predict a similar proximity to Ivantsov solutions, provided that the lengthening direction be associated with the maximum of elastic energy, as discussed in [START_REF] Appolaire | Modeling of phase transformations in β-metastable Ti alloys[END_REF]. This explains why simplified models based on Ivantsov solutions are very often successful when the tip radius is tuned (e.g. [START_REF] Appolaire | Modelling of phase transformation kinetics in ti alloys -isothermal treatments[END_REF][START_REF] Guo | Kinetics of ferrite to Widmanstätten austenite transformation in a high-strength low-alloy steel revisited[END_REF]).

The validity of Ivantsov solutions for describing Widmanstätten growth has often been questioned due to the disagreements observed with some experiments. Hence, we have compared our results with respect to available experiments to see whether a better agreement is found. For that purpose, we have first collected the measurements available in the literature. Then, we have selected those where all necessary parameters can really be obtained, in particular where tip radii are provided (Tab 1). The three resulting sets of values are plotted with open symbols in Fig. 5: ferrite (bcc) plates in austenite (fcc) in Fe-C [START_REF] Simonen | Lengthening kinetics of ferrite and bainite sideplates[END_REF] with pentagons; α (fcc) needles grown from β (bcc) in Cu-Zn [START_REF] Purdy | Widmanstätten Precipitation from non-ideal solid solution: α in β-CuZn[END_REF][START_REF] Stephens | On Widmanstätten growth from beta brass[END_REF] with diamonds; and γ (hcp, also noted as δ or ζ in different phase diagrams) plates growth from α (fcc) in Al-Ag [START_REF] Laird | The growth of γ-plates in an Al-15% Ag alloy[END_REF] with triangles. We have intentionally discarded some measurements previously compared to Ivantsov [START_REF] Hoyt | The velocity of plate precipitates growing by the ledge mechanism[END_REF]: (i) those concerning alloyed steels [START_REF] Guo | Kinetics of ferrite to Widmanstätten austenite transformation in a high-strength low-alloy steel revisited[END_REF][START_REF] Hall | Mechanism of formation of Mo 2 C needles in an Fe-C-Mo alloy[END_REF] because the analysis should ac-count for the diffusion of all species and possibly for the different conditions at the interface (local equilibrium vs. paraequilibrium); (ii) those related to the dissolution of Widmanstätten microstructure in Al-15%Ag alloy [START_REF] Sagoe-Crentsil | Dissolution kinetics of widmanstätten γ-Ag 2 Al precipitates[END_REF] because dissolution is not simply reverse growth. Unfortunately, we had also to discard a number of interesting investigations where the tip radius have been obtained indirectly from the measured kinetics using some modified Ivantsov solution [START_REF] Trivedi | Growth kinetics of plate-like precipitates[END_REF] and the maximum growth criterion (e.g. [START_REF] Kostic | Growth kinetics of bainite plates and Widmanstätten needles in the β phase of a Ag-45at.% cd alloy[END_REF]).

It seems that the experimental data are not in better agreement with our calculations than with Ivantsov solution, even if they concern different ranges of supersaturations. It can be noted that accounting arbitrarily for only 20% of uncertainty of Péclet number and supersaturation moderates the disagreements (Fig. 5). We believe that this uncertainty can be much higher owing to the difficulty to retrieve dynamical 3D features from 2D observations (most often post-mortem), as well as to estimate accurately supersaturations and diffusivities. In that respect, a few comments on the data in Tab. 1 are worth to be done. In [START_REF] Simonen | Lengthening kinetics of ferrite and bainite sideplates[END_REF], velocities and radii have been obtained from hundred of measurements with an accuracy of about 5 nm so that they can be considered as quite reliable (although Fig. 2 and Tab. III related to the measurements are not consistent in this article). The supersaturations have been estimated using the Fe-C phase diagram provided in [START_REF] Simonen | Lengthening kinetics of ferrite and bainite sideplates[END_REF], accounting for the Gibbs-Thomson effect. It must be stressed that this correction improves significantly the agreement with Ivantsov. Moreover, in view of the large supersaturations in the investigated Fe-C alloys, interactions between neighboring plates could be a reason why the data are below the Ivantsov curve. Finally, as already suggested in [START_REF] Simonen | Lengthening kinetics of ferrite and bainite sideplates[END_REF], carbon diffusivity is averaged over the concentration range in austenite. For that purpose, we have used a more recent estimate by [START_REF] Ågren | A revised expression for the diffusivity of carbon in binary Fe-C austenite[END_REF]. For the needles in Cu-Zn brass, the values provided in [START_REF] Purdy | Widmanstätten Precipitation from non-ideal solid solution: α in β-CuZn[END_REF] at 400 °C and in [START_REF] Stephens | On Widmanstätten growth from beta brass[END_REF] at 520 °C have been obtained from 20 measurements following a distribution surprisingly flat (Tab. I in [START_REF] Purdy | Widmanstätten Precipitation from non-ideal solid solution: α in β-CuZn[END_REF]). Moreover, the radii of the thinnest needles could not be measured and they have been attributed the resolution limit of 50 nm. It is not clear whether the supersaturations have been corrected to account for the curvature. Finally, the assessment of the interdiffusion coefficient entering P from the tracer diffusion determined by [START_REF] Kuper | Diffusion in ordered and disordered copper-zinc[END_REF] depends on the thermodynamic description of Cu-Zn β phase. We have used the thermodynamic factor proposed in [START_REF] Purdy | Widmanstätten Precipitation from non-ideal solid solution: α in β-CuZn[END_REF] although no detail is given on how it has been obtained. For Al-15%Ag, the tip velocities have been determined from in situ observations using hot stage transmission electron microscopy. Although not mentioned, it is likely that the number of measurements is not large due to the difficulty to perform such experiments. The radius of curvature could not be measured directly and it has been assumed to be equal to the measured plate thickness. In case where the tips are very sharp, such an assumption can overestimate the radii by one order of magnitude. Moreover, besides the uncertainty associated with the determination of the supersaturation and diffusivity already noted for the other alloys, it must be stressed that the plates fea-Alloy Ref. Temperature Supersaturation Growth velocity Tip radius Diffusivity Symbol (°C) (nm s -1 ) (nm) (×10 -14 m 2 s -1 ) in Fig. 5 Al ture very anisotropic shapes. Finally, it is also likely that growth in thin foils is modified with respect to growth in bulk.

Hence, although there is no doubt that the experimental data are invaluable, it must be recognized that they are excessively difficult to obtain and that their processing require many careful steps. Considering their insufficient accuracy, it seems they cannot be used to draw definite conclusion concerning their agreement with the Ivantsov predictions.

Tip shape and size

Two important outcomes of the previous section are:

• the tip shapes seem to be independent of the growth conditions (Fig. 6); • they are not parabolic or paraboloidal (Fig. 6) as usually assumed in analytical modeling of Widmanstätten plates.

Both statements are consequences of the role of elasticity anisotropy that can be examined with polar plots of the elastic kernel B(n).

To check that our results are not specific to the very simple eigenstrains in Eqs. ( 10) and ( 11), we have considered in this section the eigenstrain related to the transformation from bcc β phase to hcp α phase in Ti or Zr and their alloys. This eigenstrain can be expressed theoretically following the Burgers mechanism [START_REF] Burgers | On the process of transition of the cubic-body-centered modification into the hexagonal-close-packed modification of zirconium[END_REF] and using the lattice parameters of both phases. Shi et al. [START_REF] Shi | Predicting equilibrium shape of precipitates as function of coherency state[END_REF] have proposed the following one:

ε * Ti = 10 -2          
-8.3 -0.9 0 -0.9 12.3 0 0 0 3.5

          (14)
based on room temperature lattice parameters and accounting for interfacial misfit dislocations observed by TEM. It describes a shear dominated transformation like Eqs. ( 10) and ( 11), but it is also associated with some volume variation. It is worth noting that different shear amplitudes can be found depending on the way the lattice parameters are obtained. For example, an eigenstrain closer to Eq. ( 11) with ε 0 = 10% is obtained using lattice parameters measured at high temperature by synchrotron Xray diffraction in [START_REF] Settefrati | Etude expérimentale et modélisation par champ de phase de la formation de alpha dans les alliages de titane beta-métastable[END_REF][START_REF] Appolaire | Stress and strain fields associated with the formation of α" in near-β titanium alloys[END_REF][START_REF] Appolaire | Modeling of phase transformations in β-metastable Ti alloys[END_REF]. Hence, working with Eq. ( 14) will emphasize the generic nature of the results obtained with the simple eigenstrains.

Our aim is to demonstrate a correlation between the tip shape and the elastic kernel B(n) (Eq. ( 9)). To do so, we have modified the shape of the elastic kernel by artificially changing the values of the elastic constants. For convenience, we have used a cubic elastic tensor in the reference frame used in ( 14), and we have controlled the shape of the elastic kernel B(n) by modifying the anisotropy ratio A = 2C 44 /(C 11 -C 12 ). More precisely, the three components C 11 , C 12 and C 44 have been determined so as to give the same bulk modulus K = 1/3(C 11 + 2C 12 ) = 55 and the same apparent isotropic shear modulus G(C 11 , C 12 , C 44 ) = 13, as defined by Kröner homogenization scheme [START_REF] Kröner | Berechnung der elastischen konstanten des vielkristalls aus den konstanten des einkristalls[END_REF][START_REF] Ledbetter | Handbook of Elastic Properties of Solids, Liquids, and Gases[END_REF]. The shape change of the elastic kernel for the investigated parameter range (A ∈ [1; 5]) is presented in Fig. 8. B(n) always displays two pairs of minima. However, the overall shape of the elastic kernel is significantly modified between A=1 (Fig. 8a) and A=5 (Fig. 8b). Note that we observed numerically that the directions of the minima do not do not change when varying the value of A in the considered range.

To relate the elastic kernel to the tip shape, we have performed 2D simulations. The relevant plane for the 2D simulations, where growth is able to reach a stationary stage, is selected considering the polar plot of B(n) as follows [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF]. We first define e 1 as the direction of one pair of minima. e 2 is the vector normal to e 1 in the plane containing the two pairs of minima (as shown in Fig. 8, the pairs of minima are not normal to each others contrary to Fig. 1)). Then e 3 is defined as e 3 = e 1 × e 2 and lies in the habit plane by construction. The 2D simulation is then performed in the plane normal to e 2 . The relevant elastic kernel for the 2D simulations is the cut of the 3D kernel in the (e 1 , e 3 ) plane ( black circles in Fig. 8).

These sections of the elastic kernels in plane (e 1 , e 3 ) are plotted for A=1, 3 and 5 in Fig. 9 to highlight the change in the maxima directions when increasing A beyond a threshold value close to 4. For low A, the maxima are aligned with direction e 3 and B(n) decreases monotonously from e 3 to e 1 . For A=5, the maxima are at 45 • with respect to e 3 (and e 1 ) such that the variation of B(n) is not monotonous between e 3 and e 1 .

We now investigate how the shape of the elastic kernel impacts the tip shape and kinetics using 2D phase field simulations with different values of A between 1 and 5. As expected, for symmetry reasons, the lengthening always proceeds along the horizontal direction, i.e. along e 3 . The corresponding tip shapes during growth with c ∞ = 0.25 are plotted in Fig. 10. Whereas for A=1 the tip is still acicular as in the previous section, the difference is striking for A=5 with a blunted tip with a nearly flat side perpendicular to the lengthening direction. This comparison clearly shows that the tip shape is related to how the elastic kernel behaves far from the habit plane direction e 1 and cannot be predicted by the usual expansion of B(n) around e 1 (e.g. §8.1.2 and §8.1.3 in [START_REF] Khachaturyan | Theory of structural transformations in solids[END_REF]). Despite the lack of close analytical form to predict the shape, the present calculations provide some numerical evidence that acicular tip shapes are obtained when the lengthening direction e 3 , perpendicular to the habit plane normal e 1 by symmetry, corresponds to the maximum of the elastic kernel in plane (e 3 , e 1 ).

As mentioned at the beginning of the section, the tip size does not seem to depend on the growth condition for a given eigenstrain (e.g. Fig. 6). This is in contrast with dendritic growth during solidification [START_REF] Langer | Instabilities and pattern formation in crystal growth[END_REF][START_REF] Karma | Dendritic growth[END_REF] as in other diffusion-controlled growth problems where the selection is not ensured by the diffusion process [START_REF] Kassner | Continuum description of noiseless diffusionlimited aggregation[END_REF]. Indeed, the calculations of §3 have shown that the larger ε 0 the thinner the plates, with tip radii independent of the undercooling in the investigated range. It appears that this trend is not specific to the simple eigenstrains ε * 2D or ε * 3D , provided that there exists a plane where the elastic kernel features a figure eight with only one maximal direction. This shape is likely to be obtained for any shear dominated transformation with moderate elastic constant anisotropy, as illustrated above. As can be inferred from Fig. 10, the tip size seems to be related to the maximum value B max of the elastic kernel in plane (e 1 , e 3 ): the largest the maximum, the sharpest the tip. In order to substantiate this statement, we have performed calculations considering two different elastic kernels displaying figure eight (A = 1 and A = 2) with the same maximum, as shown in Fig. 11. To keep the same maximum, the isotropic bulk and shear moduli were rescaled. We have measured the steady state tip radii in both cases: they are equal within the precision of the measurement (their difference is about 8% of d). Hence, B max is the most important quantity to rationalize the tip radius selection, at least for figure eight B(n). Qualitatively, it is not surprising that the tip size results from the competition between the interface energy and some volumic elastic energy at the tip related to B max .

We have gathered the results for plates (2D) and needles (3D) for all investigated supersaturations and eigenstrains in Fig. 12, plotting the tip radius with respect to B max . The collapse of all the results into a single curve for a given eigenstrain whatever the growth conditions suggests that the tip size would be an equilibrium property. Thus, we have computed equilibrium shapes for the simple eigenstrains, considering various box sizes, discretizations, and volume fractions (ranging from 3 to 11% to avoid any interactions with periodic images). To rapidly obtain equilibrium shapes, we have solved only Allen-Cahn equation ( 5) under the constraint of constant phase fraction. For a given B max , all calculations have converged toward the same equilibrium plate shape with a single tip radius (with variations smaller than fractions of d) . In Fig. 12b for needles, the small discrepancies between the growth and equilibrium curves are smaller than the grid spacing d and can be attributed to the interactions with periodic images due to smaller boxes than in 2D. Using the analogy between plates and dislocations [START_REF] Nabarro | Cxxii. the synthesis of elastic dislocation fields[END_REF][START_REF] Khachaturyan | Theory of structural transformations in solids[END_REF][START_REF] Rodney | Phase field methods and dislocations[END_REF], it is easy to understand that lengthening is energetically more favorable than thickening, such that equilibrium plate thickness increases only very slowly with length. Moreover, although the equilibrium plate thickness depends slightly on volume fraction, the tip radius is smaller than the plate thickness and do not exhibit any variations larger than a few percents of d, in the range of volume fractions considered. Unfortunately, we were not able to perform similar equilibrium cal- culations for eigenstrain [START_REF] Horvay | Dendritic and spheroidal growth[END_REF] because of the significant elastic interactions between periodic images even in long boxes due to large B max . The resulting equilibrium values are plotted in Fig. 12 with black lines and white dots. The agreement between these equilibrium values and the steady growth values is excellent. We can conclude that the tip shape and size of Widmanstätten structures are equilibrium properties when only elasticity anisotropy is accounted for. It is likely to remain the main selecting phenomenon, even when other sources of anisotropy (interface energy and mobility) are present. Using realistic scales for the parameters ( §2), typical values of the tip radius decreases from about 25 nm for B max = O(10 -2 ) GPa to a few nm for B max = O(1) GPa. The highest B max are generated by eigenstrains with large shears, typically around 20% for ε * T i corresponding to the red triangles in Fig. 12a. Tip radii are usually reported to span only one order of magnitude from tens to hundreds of nm (Tab. 1 and [START_REF] Kostic | Growth kinetics of bainite plates and Widmanstätten needles in the β phase of a Ag-45at.% cd alloy[END_REF][START_REF] Trivedi | Growth kinetics of plate-like precipitates[END_REF][START_REF] Guo | Kinetics of ferrite to Widmanstätten austenite transformation in a high-strength low-alloy steel revisited[END_REF][START_REF] Hall | Mechanism of formation of Mo 2 C needles in an Fe-C-Mo alloy[END_REF][START_REF] Sagoe-Crentsil | Dissolution kinetics of widmanstätten γ-Ag 2 Al precipitates[END_REF]). Although the selection curves in Fig. 12 are not universal as they depend on eigenstrains, elastic constants and interface energies, they still give clues about the elastic energy densities that must be conveyed at the tips to reach radii with experimental order of magnitudes. It is clear that B max must be decreased by at least one order of magnitude for plates when considering eigenshear predicted by cristallographic analysis of transformations involving fcc, bcc and hcp structures. It is worth stressing that increasing interfacial energy does not drastically change this statement because we had already considered a significant value of σ int = 200 mJ m -2 . Different phenomena can contribute to decrease the elastic energy at the tip. First, anisotropic elastic constants can have a significant impact (Fig. 10). Second, concentration fields can also play a role when the eigenstrain depends on concentration, e.g. following Vegard's law. Although preliminary calculations have shown that the relaxation associated with this process is mod- erate in the case of eigenstrains with perfect invariant planes [START_REF] Lebbad | Modélisation de la croissance des structures de Widmanstätten par la méthode des champs de phase[END_REF], it could explain the slight variation of the tip radii measured by [START_REF] Simonen | Lengthening kinetics of ferrite and bainite sideplates[END_REF] at 700 °C with the carbon content (Tab. 1). Third, at high temperature, plasticity is likely to relax stresses at the tip. When described by isotropic J 2 constitutive laws, plasticity does not seem to change the selection process, keeping the same tip shape as shown in [START_REF] Cottura | Plastic relaxation during diffusion-controlled growth of Widmanstätten plates[END_REF]. Finally, misfit compensating interfacial defects such as disconnections are likely to contribute significantly to relax elastic energy so that the tip radius increases. They can be accounted for by introducing an effective eigenstrain as shown in [START_REF] Shi | Predicting equilibrium shape of precipitates as function of coherency state[END_REF]. Hence, we expect that accounting for the relaxation processes will not change the conclusions of the present work about the selection of the tip shape, but will decrease B max and increase the tip size in agreement with experimental observations. Still, below a sufficiently small elastic anisotropy, a crossover towards a dendrite-like growth with an-other selection criterion cannot be discarded. A careful investigation of the effect of these relaxation processes on the tip selection is under progress.

Conclusion

Considering a broader range of supersaturations and eigenstrains than in [START_REF] Cottura | Phase field study of acicular growth: Role of elasticity in Widmanstätten structure[END_REF], we have shown that the important features of Widmanstätten structures can be recovered by considering the anisotropic elastic energy generated by shear dominated eigenstrains: acicular morphologies (plates or needles) and constant lengthening rates at constant supersaturations. We have confirmed that the plates or needles lengthen at constant velocity along a particular direction perpendicular to two minima of the elastic kernel. We have also found that the Péclet number depends only on the supersaturation of the matrix and that this dependency is not far from Ivantsov analytical solutions for parabolic and paraboloidal shapes. One of the major outcome of the present study is that elasticity selects the tip shape and size, independently of the growth conditions. Consequently, when there is no other source of anisotropy than elasticity, the tip shape and size can be obtained by equilibrium calculations. We have shown numerically how the tip shape is related to the shape of the elastic kernel in the plane containing the habit plane normal and the lengthening direction. When the elastic kernel features in this plane a maximum in the lengthening direction, the tip is acicular. In this case, the tip size is shown to depend mainly on this maximum. Conversely, when, within this plane, the direction of the maximum of the elastic kernel deviates from the lengthening direction, the tip displays a blunted tip with a nearly flat side perpendicular to the lengthening direction.

The present work opens new perspectives for estimating the growth kinetics of Widmanstatten structures. First, our calculations have proved that Ivantsov solutions remain good estimates for the solute flow at the tips even when elasticity is playing a major role, changing the tip shape and the interface concentrations. Indeed, moderate differences with respect to these solutions have been found to be lower than common uncertainties on diffusion coefficients. This explains why the regular use of these solutions has been partly successful, e.g. [START_REF] Yin | C-curves for lengthening of Widmanstätten and bainitic ferrite[END_REF] for a recent example. Second, we can deduce from the present study that the tip size can be computed independently from elasticity equilibrium calculations rather than ensuing from an ad hoc maximum velocity criterion where the interface energy is tuned. Moreover, our finding suggests to revisit the available experimental data and to perform new experiments where the difficult measurements of the tip sizes could be avoided, provided that the eigenstrain can be determined. Still, further work is needed to substantiate the new selection criterion, in particular how relaxation processes such as plasticity can be accounted for in equilibrium calculations.

Figure 1 :

 1 Figure 1: Polar plots of the elastic kernel with eigenstrain given by (a) ε * 2D and (b) ε * 3D .

e 3 e 1 Figure 2 :Figure 3 :

 123 Figure 2: Level set ϕ = 0.5 during growth (last recorded time) for ε 0 = 4% and Ω = 0.15.

Figure 4 :

 4 Figure 4: Transverse section area S of the needle vs. the distance z from the tip, determined from the projections of ϕ = 0.5 onto (e 1 , e 2 ) plane (inset).

  solutions.

Figure 5 :

 5 Figure 5: Péclet number P vs tip supersaturation Ω: plates (2D) in red and needles (3D) in black. Ivantsov solutions are plotted with lines, phase field calculations with solid symbols, and experimental data with empty symbols (see also Tab. 1).

  c ∞ =0.20 c ∞ =0.25 c ∞ =0.30

Figure 6 :

 6 Figure 6: Tip shapes rescaled with the tip radius R (defined by the blue circle) during steady-state growth with ε 0 = 4% for different supersaturations. Dashed line: parabola with unit curvature radius.

Figure 7 :

 7 Figure 7: Concentration level sets in matrix close to interface (continuous and dashed black lines indicate level sets ϕ = 0.5 and ϕ = 0.05 respectively) during steady-state growth with ε 0 = 3% and Ω = 0.3.

Figure 8 :

 8 Figure 8: Polar plot of the elastic kernel with eigenstrain ε * Ti and isotropic elastic constants with (a) A = 1 and (b) A = 5. 2D computations are performed in plane (e 3 , e 1 ).

Figure 9 :

 9 Figure 9: Polar plot of the elastic kernel for ε * Ti and two Zener factors A = 1 (isotropic), A = 3 and A = 5.

Figure 10 :

 10 Figure 10: Tip shapes (ϕ = 0.5) during growth in the same conditions (c ∞ = 0.25) corresponding to Fig. 9 (the level sets have been shifted such that the tips are all at the same location).

Figure 11 :

 11 Figure 11: Polar plot of the elastic kernel for ε * Ti and two sets of elastic constants corresponding to A = 1 and A = 2 rescaled so as to display the same B max .

Figure 12 :

 12 Figure 12: Tip radius vs. B max for different supersaturations (colors) and different eigenstrains: (a) ε * 2D (circles) and ε * Ti (triangles) for plates (2D); (b) ε * 3D for needles (3D). The corresponding equilibrium values are plotted with white symbols.

Table 1 :

 1 Growth conditions and features inferred from in situ and post mortem measurements of Widmanstätten structures in different binary alloys.

	-15 wt.%Ag	[6]	400	0.034	6.50 to 7.35	20	1.08 [38]
			425	0.018	1.41 to 2.62	20	2.27 [38]
	Cu -44.1 at.%Zn [8]	400	0.2	160	130 to 250	1.8
	Cu -41.1 at.%Zn [8]	400	0.6	4.7 × 10 3	< 50	1.8
		[9]	520	0.5	30 × 10 3	100-200	160
	Fe -0.24 wt.%C			0.686	21 × 10 3	9.5	38.3 [39]
	Fe -0.33 wt.%C [35]	700	0.570	13 × 10 3	15.5	43.3 [39]
	Fe -0.43 wt.%C			0.440	6 × 10 3	18	49.6 [39]

Appendix A. Elastic kernel for ε *

2D

For isotropic elasticity, the Fourier transform of the Green function for displacements reads [START_REF] Khachaturyan | Theory of structural transformations in solids[END_REF]:

where µ is the shear modulus, and ν the Poisson ratio. For ε *

2D

Eq. ( 10), the eigenstress reads:

In the (001) plane we have n 3 = 0 such that B(n) simplifies to:

with n 1 = cos θ and n 2 = sin θ. Using Eq. (A.2) and Eq. (A.1), after tedious but trivial algebra B(n) becomes:

References