Matthieu Degeiter

Yann Le Bouar

Benoît Appolaire

Mikael Perrut

Alphonse Finel

Instabilities in the periodic arrangement of elastically interacting precipitates in nickel-base superalloys

Keywords: single-crystal superalloy, microstructure, elasticity, stability analyses, phase field, pattern formation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

In multiphase alloys, the coherent coexistence of misfitting phases generates internal elastic fields. Given their long-range and usually anisotropic character, these fields alter the kinetics of diffusion-controlled phase transformations, and also control the shapes and mutual arrangement of the precipitates. During microstructure evolution, the formation of long-range spatially-correlated patterns is therefore observed [START_REF] Khachaturyan | Theory of Structural Transformations in Solids[END_REF][START_REF] Bouar | Origin of chessboard-like structures in decomposing alloys[END_REF].

The microstructure of single-crystal nickel-based superalloys, obtained by precipitation of the L1 2 -ordered γ phase in the FCC γ matrix, is formed under the influence of elasticity. Besides modifying the interface thermochemical equilibrium [START_REF] Larché | Thermochemical equilibrium of multiphase solids under stress[END_REF][START_REF] Gurtin | The continuum mechanics of coherent two-phase elastic solids with mass transport[END_REF][START_REF] Schneider | Phase field elasticity model based on mechanical jump conditions[END_REF], elasticity causes the shape of a growing γ precipitate to gradually change from spheroid to cuboid [START_REF] Thompson | The equilibrium shape of a misfitting precipitate[END_REF]. In addition, interacting γ precipitates tend to align along the elastically soft directions of the matrix during coarsening [START_REF] Su | The dynamics of precipitate evolution in elastically stressed solids-II. Particle alignment[END_REF][START_REF] Thornton | Large-scale simulations of Ostwald ripening in elastically stressed solids: I. Development of microstructure[END_REF][START_REF] Thornton | Large-scale simulations of Ostwald ripening in elastically stressed solids: II. Coarsening kinetics and particle size distribution[END_REF]. For high precipitate volume fraction, this reconfiguration process leads to the formation of peculiar arrangements of nearly periodically aligned cuboids [START_REF] Chellman | The coarsening of γ precipitates at large volume fraction[END_REF][START_REF] Caron | Improvement of creep strength in a Nickelbase single crystal superalloy by heat treatment[END_REF][START_REF] Caron | High γ solvus new generation nickel-based superalloys for single-crystal turbine blade applications[END_REF].

However, the resulting microstructures systematically display defects in the precipitate alignments, as highlighted in Fig. 1. At high volume fraction (Fig. 1-a), two types of defects are observed on both horizontal and vertical alignments: branches, when an alignment separates into two modulations, and macro-dislocations, when a modulation terminates. Interestingly, what might initially appear as a nearly periodic arrangement of roughly aligned cuboids is thus highlighted as an intricate net-work of interweaving defects in the precipitate alignments. Other kinds of patterns are observed at lower volume fractions (Fig. 1-b). They may be identified as adjacent precipitate alignments along [100] and [010] directions, forming right-angled structures called chevrons or herringbone patterns [START_REF] Ardell | On the modulated structure of aged Ni-Al alloys[END_REF][START_REF] Loomis | The influence of molybdenum on the γ phase in experimental nickel-base superalloys[END_REF][START_REF] Biss | The effet of molybdenum on γ coarsening and on elevated-temperature hardness in some experimental nickel-base superalloys[END_REF][START_REF] Grosdidier | Precipitation and dissolution processes in γ/γ single-crystal nickel-based superalloys[END_REF]. All these alignments defects and patterns are observed both experimentally and numerically (see Ref. [START_REF] Cottura | Role of elastic inhomogeneity in the development of cuboidal microstructures in Ni-based superalloys[END_REF][START_REF] Hu | A phase-field method for evolving microstructures with strong elastic inhomogeneity[END_REF][START_REF] Zhu | Morphological evolution during phase separation and coarsening with strong inhomogeneous elasticity[END_REF] for forks and dislocation patterns at high volume fraction, and Ref. [START_REF] Thornton | Large-scale simulations of Ostwald ripening in elastically stressed solids: I. Development of microstructure[END_REF][START_REF] Vaithyanathan | Coarsening of ordered intermetallic precipitates with coherency stress[END_REF][START_REF] Wang | Field kinetic model and computer simulation of precipitation of L1 2 ordered intermetallics from F.C.C. solid solution[END_REF] for herringbone patterns at lower volume fractions).

Though it is well-established that the anisotropy of the elastic moduli is responsible for the formation of these specific arrangements, and despite a significant number of studies devoted to the γ/γ microstructure, the origin and evolution of the pattern defects remain to be clarified.

In this contribution, we address the question of the origin of the pattern defects in the γ/γ microstructure, and show that their formation is driven by specific shapedependent instabilities of the periodic arrangement of interacting precipitates. This paper is organized as follows. The general formalism used to perform the static stability analysis of a periodic arrangement of misfitting inclusions is presented in Section 2. This formalism is then applied in Section 3 to successively analyze the importance of the precipitate shape, the volume fraction and the elastic anisotropy. A close attention is paid to the instability modes of a cubic arrangement of cubic precipitates. In Section 4, an elastic phase field model is used to analyze the dynamics resulting from of the two main instability modes predicted by the Preprint submitted to Acta Materialia [START_REF] Diologent | Comportement en fluage et en traction de superalliages monocristallins à base de nickel[END_REF]; (b) Chevron patterns in a binary model nickel-based superalloy [START_REF] Li | Regular γ/γ phase interface instability in a binary model nickel-based single-crystal alloy[END_REF]. static stability analysis. The specific patterns generated by these modes are presented and compared to the microstructures observed in Ni-base superalloys. We finally investigate the consequences of an inhomogeneity on the C shear modulus, which has been shown to significantly influence the microstructure evolution in these alloys [START_REF] Cottura | Role of elastic inhomogeneity in the development of cuboidal microstructures in Ni-based superalloys[END_REF].

Stability analysis of periodic arrangements

The subsequent analytical development is based on Khachaturyan's theory of elasticity [START_REF] Khachaturyan | Theory of Structural Transformations in Solids[END_REF], as well as on the work of Khachaturyan and Airapetyan [START_REF] Khachaturyan | Spatially Periodic Distributions of New Phase Inclusions Caused by Elastic Distorsions[END_REF]. Considering a coherent two-phase mixture in stress-free boundary condition, and assuming small strains and elastic homogeneity, the strain energy of the microstructure reads

E el = 1 2 V dVλ ijkl (ε ij (r) -ε 0 ij (r))(ε kl (r) -ε 0 kl (r)) (1)
where λ ijkl is elastic tensor, ε ij (r) is the total strain, and

ε 0 ij (r) = ε 0 ij θ(r)
is the stress-free strain tensor associated to the transformation eigenstrain ε 0 ij . θ(r) is the shape function equal to unity if r lies inside a precipitate, and zero otherwise. To simplify the notation, we assume that the volume V is orthorhombic with periodic boundary conditions.

At mechanical equilibrium, the strain energy takes the form [START_REF] Khachaturyan | Theory of Structural Transformations in Solids[END_REF]:

E el = V 2 k =0 B(n)|θ(k)| 2 (2)
where n = k/k is a unit vector of reciprocal space. The Fourier transform of θ(r) is defined by

θ(k) = 1 V V dV θ(r)e -ikr . (3)
The elastic kernel B(n) is:

B(n) = ε 0 ij λ ijkl ε 0 kl -n i σ 0 ij Ω jk (n)σ 0 kl n l ( 4
)
where σ 0 ij = λ ijkl ε 0 kl is the transformation stress tensor, and Ω jk (n) is the normalized Green tensor defined by

Ω -1 ij (n) = λ iklj n k n l .
Let a distribution of N identical precipitates be periodically arranged along the x, y and z directions of a simple cubic macrolattice, with a the macrolattice parameter. The aim of the present calculation is to analyze the stability of this arrangement with respect to small disturbances of the precipitate positions, that we will call displacements (not to be confused with material displacements). The microstructure is completely described by the shape function:

θ(r) = R θ 0 (r -R) (5)
where θ 0 is the shape function of a single precipitate sitting at the origin of the macrolattice and where the sum runs over all the precipitate positions R. In reciprocal space, θ becomes:

θ(k) = θ 0 (k) R e -ikR . (6)
The strain energy of the microstructure becomes

E el = V 2 k B(n)|θ 0 (k)| 2 R,R e -ik(R-R ) (7)
where we have set B(0)=0.

The double sum on the right-hand side of Eq. 7 can be divided into a sum over R = R corresponding to the selfenergy of the precipitates, and a sum over R = R which carries the dependence of the total strain energy to the spatial arrangement of precipitates. As a conclusion, the stability of the arrangement is related to the configuration energy defined by

E conf el = V 2 k B(n)|θ 0 (k)| 2 R,R R =R e -ik(R-R ) (8)
We now consider that the position R of each precipitate is close to some position R 0 in a perfect cubic macrolattice and we define the displacements vectors u(R 0 ) = R -R 0 .

Assuming small displacements, the configuration energy change between the perfect and the distorted arrangements is, to the second order,

∆E conf el = V 2 k B(n)|θ 0 (k)| 2 R0,R 0 R0 =R 0 - 1 2 k i k j u i (R 0 ) -u i (R 0 ) u j (R 0 ) -u j (R 0 ) e -ik(R0-R 0 ) . ( 9
)
As usual with periodic structures, reciprocal lattice vectors are decomposed as k = H + τ where H is a vector of the reciprocal lattice and τ is a vector of the first Brillouin zone of the macrolattice.

Using this decomposition we have

R0 e -ikR0 = R0 e iτ R0 = Nδ τ =0 . (10)
In addition, since u(R 0 ) is defined at the lattice sites, we define its forward and backward Fourier transforms as:

v i (τ ) = 1 N R0 u i (R 0 )e -iτ R0 (11)
u i (R 0 ) = τ v i (τ )e iτ R0 . (12)
Introducing Eqs. ( 12) and [START_REF] Chellman | The coarsening of γ precipitates at large volume fraction[END_REF] in Eq.( 9) leads to the following closed form for the configuration energy change:

∆E conf el = N 2 τ κ ij (τ )v i (τ )v * j (τ ) (13)
where the curvature tensor for a unit cell of the periodic arrangement is defined by

κ ij (τ ) = V N H B H + τ H + τ (H i + τ i )(H j + τ j )|θ 0 (H + τ )| 2 -B H H H i H j |θ 0 (H)| 2 . ( 14
)
It follows from ( 13) that the configuration energy change is a quadratic form of independent Fourier components of the displacements. Therefore, the configuration of identical precipitates on a macrolattice is stable with respect to the precipitate displacements if the tensor κ ij (τ ) is positive for each Fourier mode τ . As usual, the positiveness of the tensor is analyzed by computing its eigenvalues. If all eigenvalues are positive, the configuration is stable. If any eigenvalue is negative, the configuration is unstable when the arrangement is perturbed by a displacement of mode τ along the corresponding eigenvector.

For the cubic macrolattice considered here, the first Brillouin zone is the set of vectors τ such that -π/a < τ i ≤ π/a. For the numerical computation of Eq. ( 14), the summation over reciprocal lattice sites H has to be truncated. In our calculations, the summation is carried over

H i = 2πα i /a, where α i is an integer such that -n RL ≤ α i ≤ n RL .
The number of replicated Brillouin zones n RL is chosen depending on the precipitate shape, to ensure sufficient precision. Given the Brillouin zone symmetry, the eigenvalues of κ ij (τ ) are computed for the wave vectors τ along the high symmetry directions. These directions are defined between the symmetry points Γ(0, 0, 0), X(π/a, 0, 0), M (π/a, π/a, 0) and R(π/a, π/a, π/a). Since |τ | = 2π/Λ, where Λ is the wavelength, Γ corresponds to an infinitewavelength perturbation (i.e. a translation), which induces no change in configuration energy (κ ij (Γ) = 0). The perturbations at the other symmetry points have wave vector components equal to 0 or π/a and span over two precipitates.

Stability of a cubic array of precipitates

In this section, stability analyses of the periodic distribution of elastically interacting precipitates are carried out. The importance of the precipitate shape, of the volume fraction and of the elastic anisotropy are successively analyzed.

Shape-dependent instabilities in the arrangement

We first consider the case of spherical inclusions of radius R for which the shape function in Fourier space is

θ 0 (k) = 4 3 πR 3 V 3 kR cos(kR) -sin(kR) (kR) 3 . ( 15
)
We choose the macrolattice parameter a = 350 nm, and the volume fraction f v = 37 %. The transformation eigenstrain is ε 0 = 0.48 % and the elastic constants are C 11 = 250 GPa, C 12 = 160 GPa and C 44 = 118.5 GPa.

All the above parameters have been selected to allow a direct comparison with a previous work devoted to the stability of cubic arrangements of spherical precipitates [START_REF] Khachaturyan | Spatially Periodic Distributions of New Phase Inclusions Caused by Elastic Distorsions[END_REF]. The computation of the curvature tensor at each Brillouin zone wave vector was carried out with n RL = 500, to ensure sufficient precision. The three eigenvalues κ 0 α of κ ij (τ ) are plotted in Fig. 2a against the Brillouin zone wavevector τ of the cubic macrolattice, along the high-symmetry directions. The longitudinal κ 0 L (i.e. eigenvector parallel to τ ) and transverse transverse κ 0 T1 and κ 0 T2 (i.e. eigenvector normal to τ ) modes are indicated. Finally, Fig. 2-b displays, for each branch, a sketch of the deformation of a cubic array of spheres for wave vectors with ζ = 2π/(12a). The branches of the curvature spectrum in Fig. 2-a are positive for all Brillouin zone wave vectors along the high-symmetry directions, implying that all the perturbations considered lead to an increase in the configuration energy. We have also verified this property for all vectors of the first Brillouin zone (not shown). So the periodic distribution of spherical inclusions is stable with respect to displacements from their ideal positions. This conclusion is similar to that of Khachaturyan and Airapetyan [START_REF] Khachaturyan | Spatially Periodic Distributions of New Phase Inclusions Caused by Elastic Distorsions[END_REF] 1 .

Γ τ 1 X M τ 2 Γ τ 3 R -2 0 2 4 6 8 κ 0 α (τ ) (J.m -2 ) [ζ00] [ζζ0] [ζζζ] L T 1,2 L T 1 T 2 L T 1,2 (a) (b) κ 0 L (τ 1) κ 0 T1 (τ 1) κ 0 T2 (τ 1) κ 0 L (τ 2) κ 0 T1 (τ 2) κ 0 T2 (τ 2) κ 0 L (τ 3) κ 0 T1 (τ 3) κ 0 T2 (τ 3)
We now consider the case of cubic precipitates. The shape function in Fourier space becomes

θ 0 (k) = c 3 V sin(k x c/2) k x c/2 sin(k y c/2) k y c/2 sin(k z c/2) k z c/2 , ( 16
)
where c is the edge length of the cube. We have considered edge lengths which are multiple of a/(2n RL + 1) and reciprocal lattice vectors τ with components which are multiple of 2π/(a n p ), where n p is an integer. This choice corresponds to a simulation box of length a n p discretized in n p (2n RL + 1) grid points, with periodic conditions. Because the edge lengths are a multiple of the grid spacing, this choice ensures a rapid convergence such that the number of replicated Brillouin zones can be reduced to n RL = 50. All other parameters were kept the same as for the periodic distribution of spheres, and the corresponding stability landscape is presented in Fig. 3-a.

As clearly visible on the spectrum in Fig. 3-a, three branches of the stability landscape become negative when considering cubes rather than spheres. Along the [100] direction, the longitudinal branch L is negative for all wavelengths. For the [ζζ0] wave vectors, the transverse branch T 2 is positive for short wavelengths (i.e. close to M ), changes its sign for ζ around 2π/(3.4a) and remains negative for greater wavelengths. Finally, along the [111] direction, the transverse branches are slightly negative for large wavelengths. These calculations demonstrate that the periodic arrangement of cubic precipitates on a cubic macrolattice is unstable for the volume fraction f v = 37%.

The instability, which displays the most negative value of κ 0 α , sits along the longitudinal [ζ00] branch, illustrated at the top of Fig. 3-b. Two features can be brought about to explain why. First, it has the property to maintain the alignments of precipitates along the cubic directions, which are the elastic soft directions. This instability mode also has the characteristics of avoiding the alignments of precipitates along the 110 directions. These directions are elastically unfavorable, as can be deduced from the configurational force maps [START_REF] Su | The dynamics of precipitate evolution in elastically stressed solids-II. Particle alignment[END_REF][START_REF] Degeiter | Etude numérique de la dynamique des défauts d'alignement des précipités γ dans les superalliages monocristallins à base de nickel[END_REF]. Note that the instability of the γ/γ microstructure with respect to this mode has already been reported in two-dimensional simulations [START_REF] Degeiter | A new analysis of the microstructure of Ni-based single-crystal superalloys : relevant topological parameters for efficient microstructural modeling[END_REF]. In pattern formation, the long-wavelength longitudinal instability of lamellar patterns is called Eckhaus instability [START_REF] Cross | Pattern formation outside of equilibrium[END_REF]. In this context, this instability has been shown to induce the formation of topological defects such as branches and macro-dislocations [START_REF] Sakaguchi | Defect Creation by the Eckhaus Instability[END_REF][START_REF] Rasenat | Eckhaus instability and defect nucleation in two-dimensional anisotropic systems[END_REF][START_REF] Newell | Instabilities of Dislocations in Fluid Patterns[END_REF][START_REF] Hu | Convection under rotation for Prandtl numbers near 1: Linear stability, wave-number selection, and pattern dynamics[END_REF][START_REF] Peña | Transverse instabilities in chemical Turing patterns of stripes[END_REF]. During the formation and growth of cuboidal microstructures, the longitudinal instability may also be at the origin of the formation of topological defects, in particular through the pinching of precipitate rows, in the same way as the pinching observed in Rayleigh-Bénard convection patterns [START_REF] Sakaguchi | Defect Creation by the Eckhaus Instability[END_REF].

The second instability is the long wavelength [ζζ0] transverse branch κ 0 T2 , illustrated at the bottom of Fig. 3-b. The direction of the displacements of the precipitate positions,

Γ τ 1 X M τ 2 Γ R -2 -1 0 1 2 3 4 5 κ 0 α (τ ) (J.m -2 ) [ζ00] [ζζ0] [ζζζ] L T 2 T 1,2 (a)
(b) given by the eigenvector of the curvature tensor, is [1 10]. This instability preserves the (001) planes of the macrolattice. Inside these planes, the macrolattice is alternatively rotated clockwise and counter-clockwise. The rotation is supplemented by an elongation of the macrolattice (at constant volume) either in the [100] or in the [010] direction.

κ 0 L (τ 1) κ 0 T2 (τ 2)
The resulting configuration at the bottom of Fig. 3-b is similar to a herringbone pattern and could initiate the formation of chevron patterns. This point will be further addressed in Sec. 4. This instability mode slightly disrupts the favorable alignments between the first neighbor precipitates in the elastically soft directions. On the other hand, this mode has the advantage of displacing the second neighbor precipitates, initially aligned along the elastically unfavorable 110 directions, in a transverse direction, which increases their distance to the unfavorable alignment.

The third instability is due to the transverse branches along the [ζζζ] direction. Because this instability is significantly weaker than the two others, it will not be further studied in this work.

To sum up, the calculations with 37% of precipitates have highlighted two important points:

• When considering cubic anisotropy and elastic soft directions along the cubic directions, the perfect arrangement of cubic precipitates is not necessarily stable.

• The approximation of spherical precipitates, used in a previous work, is not accurate enough to analyze the stability of the precipitate arrangement in Ni-base superalloys.

Rôle of the volume fraction on the stability

The effect of the volume fraction of γ precipitates on the stability landscape is now analyzed. For that purpose, curvature spectra were computed for increasing value of f v , ranging from 5% to 95% within steps of 5%. The three eigenvalues were organized according to the longitudinal L and the transverses T 1 and T 2 branches. This classification, in principle not valid between X and M , can easily be extended in this region using the continuity of the branches. Eigenvalues for the longitudinal L and transverse T 2 modes corresponding to the main instabilities are plotted with respect to f v in Fig. 4-a and Fig. 4-b, respectively.

For volume fractions greater than 0.2, the stability landscape is not very sensitive to the volume fraction: indeed, instability regions similar to that obtained for f v = 0.37 (Fig. 3-a) are observed. In particular, the longitudinal [ζ00] mode is always unstable as indicated by the negative eigenvalues for all wavelengths between Γ and X whatever the volume fraction (Fig. 4-a). On the contrary, the wavelengths range of negative eigenvalues for the transverse [ζζ0] T 2 mode (between M and Γ) shrinks by about a factor of 2 when f v increases from 0.3 to 0.9. Thus, according to this stability analysis, chevron and herringbone lengths are expected to increase with the volume fraction.

Rôle of elastic tensor anisotropy

The rôle of the elastic tensor anisotropy on the stability of the periodic arrangement is now investigated. We have used the elastic constants C 11 = 250 GPa, C 12 = 160 GPa, and we have computed C 44 so that the anisotropy parameter ξ = (C 11 -C 12 -2C 44 )/C 44 takes the targeted values of: 0, -0.5, -1, -1.5. The transformation eigenstrain was kept at ε 0 = 0.48 %, and the volume fraction at f v = 37 %.

The elastic anisotropy appears clearly in the polar plots of the elastic kernel B(n) shown in Fig. 5-a. Minima are located along the 100 directions and maxima along the 111 directions. 110 directions are also high energy directions. Because the elastic interaction energy of two coherent precipitates aligned along a direction n is, at long distance, proportional to B(n), these polar plots illustrates the favorable alignments along the cubic directions as well as the unfavorable alignments along the 110 and 111 directions.

The stability spectra are presented in Fig. 5-b. At ξ = 0, the curvature spectrum is flat and equal to 0, which means that elastic energy is insensitive to perturbations. This is consistent with the Bitter-Crum theorem [START_REF] Bitter | On Impurities in Metals[END_REF][START_REF] Cahn | A simple model for coherent equilibrium[END_REF]: the strain energy of an elastically homogeneous and isotropic coherent mixture depends neither on the shapes nor on the spatial distribution of the precipitates, but only on their volume fraction.

For ξ < 0, the spectra all display the same three in-(a) stability ranges observed in Fig. 3-a, with the same hierarchy from the most destabilizing mode to the weakest:

ξ = 0 ξ = -0.5 ξ = -1 ξ = -1.5 Γ X M Γ R -2 0 2 4 6 κ 0 α (τ ) (J.m -2 ) [ζ00] [ζζ0] [ζζζ] (b) L T 1,2 L T 1 T 2 L T 1,2 ξ = 0 ξ = -0.5 ξ = -1. ξ = -1.5
[ζ00] longitudinal, [ζζ0] transverse T 2 and [ζζζ] transverse. In Fig. 5-b, it also appears that a decrease of ξ leads to an extension of instability regions as well as a decrease of the related eigenvalue. More precisely, whatever the negative value of the anisotropy, the longitudinal mode along [ζ00] is unstable for all wavelength, but a decrease of the anisotropy progressively reduces the instability range of the other mode to larger wavelengths. Increasing the anisotropy |ξ| increases the absolute value of the negative eigenvalues as well as the range of the unstable wavelengths. Hence, the longitudinal [ζ00] modes are still all unstable.

We therefore come to the conclusion that, whereas the anisotropy of the elastic tensor is at the origin of the alignments of precipitates along the cubic directions, it is also responsible for the instability of the simple cubic macrolattice.

Instability-induced microstructure evolutions

Phase field model

The stability analysis of the previous section relies on the energy variation when the reference configuration is disturbed by displacing precipitates. This calculation has revealed the instability of the reference microstructure. However, as it is purely energetic and does not include any kinetic processes, this calculation gives no information on the growth rate of the destabilizing modes. In fact, this approach features three limitations: (i) it gives no information on the growth rate of the destabilizing modes; (ii) it is limited to infinitesimal precipitate displacements; (iii) it relies on prescribed shapes for the precipitates, while elastic interactions are known to induce morphological evolutions [START_REF] Su | The dynamics of precipitate evolution in elastically stressed solids-II. Particle alignment[END_REF]. To overcome these limitations, we have undertaken a kinetic study of the microstructure evolution using a phase field model.

Phase field models have been extensively used to analyze the microstructure evolution in nickel-base superalloys [START_REF] Vaithyanathan | Coarsening of ordered intermetallic precipitates with coherency stress[END_REF][START_REF] Wang | Field kinetic model and computer simulation of precipitation of L1 2 ordered intermetallics from F.C.C. solid solution[END_REF][START_REF] Boisse | Phase-field simulation of coarsening of γ precipitates in an ordered γ matrix[END_REF][START_REF] Boussinot | Phase-field modeling of bimodal microstructures in nickel-based superalloys[END_REF][START_REF] Gaubert | Coupling phase field and viscoplasticity to study rafting in Ni-based superalloys[END_REF][START_REF] Zhou | Large-Scale Three-Dimensional Phase Field Simulation of γ rafting and creep deformation[END_REF][START_REF] Mushongera | Effect of Re on directional γ -coarsening in commercial single crystal Ni-base superalloys: A phase field study[END_REF][START_REF] Cottura | Coupling the Phase Field Method for diffusive transformations with dislocation density-based crystal plasticity: Application to Ni-based superalloys[END_REF][START_REF] Goerler | Topological phase inversion after longterm thermal exposure of nickel-base superalloys: Experiment and phase-field simulation[END_REF]. We have used the classical phase field model detailed in Refs. [START_REF] Degeiter | A new analysis of the microstructure of Ni-based single-crystal superalloys : relevant topological parameters for efficient microstructural modeling[END_REF][START_REF] Boussinot | Phase-field modeling of bimodal microstructures in nickel-based superalloys[END_REF]. In this model, the two phases alloy is described by a single concentration field and a homogeneous free energy density approximated by a standard double-well potential. This model is not able to describe the coalescence of ordered precipitates but is sufficient to analyze the microstructure evolutions considered in the present study. Also we mention that the model is constructed in such a way that, under stress, interfaces perpendicular to the elastically soft directions adopt a predefined profile (here (1+tanh(2x/w))/2, where w is the interface width) with a predefined surface tension [START_REF] Boussinot | Phase-field modeling of bimodal microstructures in nickel-based superalloys[END_REF]. The parameters of the phase field model are taken from [START_REF] Degeiter | A new analysis of the microstructure of Ni-based single-crystal superalloys : relevant topological parameters for efficient microstructural modeling[END_REF] and correspond to a simplified Ni-base superalloy at 950 • C. In this section, the elastic constants are C 11 = 197 GPa, C 12 = 144 GPa and C 44 = 90 GPa, and the precipitate volume fraction is f v = 37% as in §3.1. The eigenstrain is ε 0 = 0.059 and the interface energy, assumed isotropic, is σ = 10 mJ/m 2 . Times are given in reduced units using the characteristic time t 0 = 8d 2 /D where d is the discretization step and D is the diffusion coefficient. The interface width is w = 3.2d, corresponding to about 7 points in the interface defined as the interval between 1% and 99% of the concentration jump.

Elastically homogeneous case

Two kind of calculations have been performed to analyze the two main instabilities identified previously. In both cases, the initial configuration is set as an arrangement of cubic precipitates, whose positions have been shifted with respect to the periodic reference according to a destabilizing mode. The amplitude and wavelength of the modes were chosen so that the displacements are multiple of the grid spacing d. To smooth the initial interfaces, the concentration field is convolved with a gaussian function.

We first consider a microstructure with 8 precipitates inside a 3.7×0.5×0.5 µm 3 box discretized with d = 3.6nm, and we impose periodic boundary conditions along the three cubic directions. This microstructure is initially perturbed with a [ζ00] longitudinal mode. If a denotes the initial periodicity of the macrolattice, the wavelength of the initial mode is Λ = 8a, which corresponds to ζ = 2π/8a, and the initial amplitude is 28.8 nm (8 grid spacings). The evolution is presented in Fig. 6 at three different times, the last one preceding coagulation of neighboring precipitates.

At t * =0 in Fig. 6-a, the precipitates are cubes with diffuse interfaces, and the imposed slight variation in their longitudinal spacings is difficult to see: the spacings between the precipitates on the left-and right-hand sides of the simulation box are smaller and larger than in reference arrangement, respectively.

As shown in Fig. 6-b andc, and in agreement with the stability analysis, the dynamics does not bring the microstructure back to the reference alignment of equidistant cubic precipitates. More precisely, the center of mass of the precipitates does not significantly move and rather we observe a morphological evolution of the precipitates:

(i) The precipitates that have been made closer to each other by the perturbation along [100] elongate progressively in this direction to form rectangular cuboids. This evolution is similar to the one observed in 2D simulations by Su and Voorhees [START_REF] Su | The dynamics of precipitate evolution in elastically stressed solids-II. Particle alignment[END_REF] during the interaction of two precipitates. Note that the coexistence of all three kinds of shapes has been reported in the 3D characterization of a Ni-base superalloy [START_REF] Lund | A quantitative assessment of the three-dimensional microstructure of a γ-γ alloy[END_REF].

Just before the first coalescence event (t * = 18000), the shapes of precipitates that are about to merge exhibit con-cave sides, as a result of the configurational forces that are attractive near the corners and repulsive near the center of the faces [START_REF] Lund | The effects of elastic stress on coarsening in the Ni-Al system[END_REF].

It is worth stressing that the simulation in Fig. 6 is a simplified case with a 1D [100] longitudinal instability with periodic images along the [010] and [001] directions. As mentioned earlier, this kind of instability could be associated with the formation of branches and macro-dislocations because it is similar to the pinching mechanism observed in the Eckhaus instability in Rayleigh-Bénard convection cells pattern [START_REF] Sakaguchi | Defect Creation by the Eckhaus Instability[END_REF]. Although the full description of such a process would require much larger simulation boxes, we believe that the existence of the [100] longitudinal instability could be responsible for the formation of the defects highlighted in Fig. 1-a.

The second mode investigated is the T 2 transverse along the [ζζ0] direction, imposed on a reference arrangement of 4×4×1 precipitates, periodically distributed in a simulation box of volume 1.8×1.8×0.5 µm 3 discretized with d=3.6nm. The initial perturbation mode spans over 4 precipitates along [100] and [010] (ζ = 2π/4a) with an amplitude of 36 nm (10 grid spacings).

The evolution of this perturbed configuration during coarsening is presented in Fig. 7.

As expected from the linear stability analysis, the precipitate centers do not go back to the sites of the reference cubic macrolattice. In fact, the phase field calculation reveals that the instability induces a significant shape evolution of the precipitates while their center of mass remains at their initial perturbed position during the simulation. Due to the small wavelength considered in this simulation, only two types of morphologies are observed: plate-like and irregular. The formation of platelet-like precipitates creates large matrix channels alternatively in the vertical and horizontal directions (see 2D sections in Fig. 7-b andc). During the microstructure evolution, we also observe that the precipitates that are close to each other align along cubic directions, which are the elastic soft directions (e.g. [START_REF] Su | The dynamics of precipitate evolution in elastically stressed solids-II. Particle alignment[END_REF]). At the end of the simulation, a herringbone pattern is formed, which is similar to chevron patterns observed in Fig. 1-b. Note that chevron patterns have also been reported in two-dimensional phase field simulations [START_REF] Vaithyanathan | Coarsening of ordered intermetallic precipitates with coherency stress[END_REF][START_REF] Wang | Field kinetic model and computer simulation of precipitation of L1 2 ordered intermetallics from F.C.C. solid solution[END_REF].

As a conclusion, beyond the stability analysis, the phase field calculations have shown that the development of the instabilities does not proceed by the movement of the precipitates but by a change of their morphologies. The instabilities lead to the coexistence of precipitates with different shapes and to the formation of specific microstructures that can be observed in electron micrographs.

Inhomogeneity on C elastic constant

Linear stability analysis

The influence of an elastic inhomogeneity ∆C = (C γ -C γ )/C γ is now investigated by comparing calculations When considering an elastic inhomogeneity, the curvature tensor cannot be computed with Eq. ( 14). We have therefore used the FFT mechanical solver of the phase field code to compute the energy of configurations perturbed by sinusoidal waves of increasing amplitudes. The energy increases with respect to the amplitudes are fitted with quadratic functions so as to retrieve the components of the curvature tensor. In the following, we have only analyzed the two main instabilities identified in the elastically homogeneous case. Note finally, that the FFT solver that we have used is best suited to diffuse interfaces. Therefore, as for the initial configuration of the phase field calculations (see Sec.4), the interfaces have been smoothed by applying a small gaussian filter. Considering sharp interfaces would require more sophisticated solvers [START_REF] Ruffini | Three-dimensional phase-field model of dislocations for a heterogeneous face-centered cubic crystals[END_REF][START_REF] Valdenaire | Crystal plasticity: Transport equations and dislocation density[END_REF] that are beyond the scope of the present work.

Γ X M Γ R -2 -1 0 1 2 3 4 5 κ 0 α (τ ) (J.m -2 ) [ζ00] [ζζ0] [ζζζ]
(a) The continuous lines are obtained from Eq. ( 14). Blue circles and red squares correspond to curvatures computed using the FFT solver with ∆C = 0% and ∆C = 50%, respectively.

∆C = 0 % ∆C = 50 % Γ X M Γ R -2 0
The spectra are presented in Fig 8-a and b for the two volume fractions f v = 37 % and f v = 60 % and compared to the homogeneous spectrum obtained with Eq. [START_REF] Loomis | The influence of molybdenum on the γ phase in experimental nickel-base superalloys[END_REF].

First, we have checked that smoothing the interfaces does not change drastically the spectra for the homogeneous elastic constants. This is the case as demonstrated by the blue dots obtained with diffuse interfaces very close to the black lines. The differences are only visible near the edges of the Brillouin zone.

Then, the influence of the elastic inhomogeneity ∆C can be deduced by comparing the spectra for ∆C = 0 (blue circles) and for ∆C = 50% (red squares). For the [ζ00] modes, our calculations show that increasing ∆C induces a decrease of all the branches. In particular, the longitudinal mode becomes more unstable. For the [ζ ζ 0] modes, the consequences of the elastic inhomogeneity ∆C are more complex and depend both on the volume fraction and on the considered branch. When focusing on the destabilizing transverse T 2 mode, we observe that an increase of ∆C leads to a very slight increase of the destabilizing wavelength range for f v = 37 % whereas a decrease is observed for f v = 60 %. This suggests that increasing ∆C would delay the formation of chevron and herringbone patterns in high volume fraction alloys.

Microstructure evolution

In this subsection, we carry out phase-field calculations to further investigate the two main instabilities of the periodic arrangement of cubic precipitates with an elastic inhomogeneity on the C shear modulus. This allows us to go beyond the static stability analysis performed in the previous section and, in particular, such calculations have the advantage to handle shape evolutions, a key ingredient as already shown in §(Sec. 4). In this section, we use a high volume fraction (f v = 60 %) corresponding to usual monocrystalline nickel-based superalloys.

First, we handle the longitudinal mode along the [100] direction, with a wavelength Λ=16a and an amplitude of 36 nm (Fig. 9). This mode is applied on a microstructure consisting of 16×1×1 precipitates periodically arranged in a simulation box of volume 7.4×0.5×0.5 µm 3 discretized with d=3.6nm.

The movement of the precipitates is analyzed using a new method proposed in Ref. [START_REF] Degeiter | A new analysis of the microstructure of Ni-based single-crystal superalloys : relevant topological parameters for efficient microstructural modeling[END_REF]. It relies on the phase ψ x (r) (not to be confused with phase field φ(r)) which is a continuous equivalent of the discrete phase shift ψ x (R 0 ) = k x u x (R 0 ) induced by the movement of one precipitate along [100]. This quantity is plotted at different time steps (Fig. 9-bottom) to follow the precipitates movement along [100].

To highlight the consequences of elastic inhomogeneity, the microstructure evolution with ∆C = 50% (Fig. 9-b) is compared to the microstructure obtained with homogeneous elasticity (Fig. 9-a).

In the homogeneous case (Fig. 9-a), the evolution is similar to the one observed in Fig. 6 for a lower volume fraction: the precipitate positions remain the same while their shapes are evolving up to a point where a coagulation event occurs (t * ≈ 1.5 10 4 ). In the inhomogeneous case (Fig. 9b), the microstructure evolution is much slower. At the time when a coagulation event occurs in the homogeneous simulation, no significant evolution of the precipitate positions or shapes can be measured in the inhomogeneous simulation. For much longer times, the precipitates move and their shape evolves, but the amplitude of the initial mode neither grows nor decays during the simulation, and appears to be kinetically frozen. The movement of the precipitates results from the development of other longitudinal modes with shorter wavelength, as it is clearly visible on the phase profiles presented in the bottom of Fig. 9-b. Note that a similar behavior was reported in 2D simulations [START_REF] Degeiter | A new analysis of the microstructure of Ni-based single-crystal superalloys : relevant topological parameters for efficient microstructural modeling[END_REF].

The behavior of the elastically inhomogeneous system can be rationalized using the following results from the literature. First, when the precipitates are harder than the matrix, as it is the case here for the C component, the elastic inhomogeneity has been shown to favor compact precipitate shapes [START_REF] Barnett | The strain energy of a coherent spherical precipitate[END_REF][START_REF] Lee | A study on coherency strain and particle morphology via discrete atom method[END_REF][START_REF] Kato | Elastic strain energy of sphere, plate and needle inclusions[END_REF][START_REF] Li | Two-and three-dimensional equilibrium morphology of a misfitting particle and the Gibbs-Thomson effect[END_REF]. Hence, shape changes observed in the homogeneous simulations are less favorable in the inhomogeneous case. Second, elastic inhomogeneity has been reported to slow down the evolution kinetics. This can explain why the evolution is slower when ∆C = 50% (Fig. 9-b) than when ∆C = 0 (Fig. 9-a).

We now analyze the evolution of a configuration initially perturbed with the transverse T 2 mode along the [ζζ0] direction. We consider a periodic microstructure composed of 4 × 4 × 1 precipitates in a 1.8 × 1.8 × 0.5 µm 3 box discretized with d = 3.6nm. The initial amplitude is 36 nm and the wave vector component ζ = 2π/4a is chosen close to the limit of the instability region predicted by the linear analysis (Fig. 8). More precisely, the wavelength is 6% and 16% lower than the stability limit in the homogeneous and inhomogeneous cases, respectively. Therefore, within the framework of the linear stability analysis considering constant precipitate shapes, the configuration is stable in both the homogeneous and inhomogenous cases.

The final configuration of the phase field simulations, presented in Fig. 10-a, demonstrates that, in the simulation with homogeneous elasticity, the initial perturbation induces precipitate shape changes and leads to the formation of a herringbone pattern, as in Fig. 7. It means that the energy gain associated with the morphological evolution exceeds the slight energy increase that is predicted by the linear analysis performed at constant precipitate shape.

The microstructure evolution is very different with inhomogeneous elasticity (Fig. 10-b). Indeed, the initial perturbation progressively vanishes and a perfect cubic macrolattice is recovered. Two reasons can be given to explain the difference in behavior between homogeneous and inhomogeneous calculations. First the relevant eigenvalue of the curvature tensor is slightly higher for the inhomogeneous case (Fig. 8-b). as mentioned above, when the precipitates are harder than the matrix, the elastic inhomogeneity favors compact precipitate shapes.

As a conclusion, we have shown with two examples that the inhomogeneity of the shear modulus C may qualitatively change the microstructure evolution resulting from some instability. Therefore defects related to these instabilities (chevron patterns, macro-dislocations, ...) are expected to be sensitive to this inhomogeneity coefficient.

Conclusions

In this contribution, the origin of the defects in the precipitate alignments of the γ/γ microstructure in nickelbase superalloys has been addressed, by means of semianalytical stability analysis and phase field calculations.

In the context of microstructure formation in nickelbased superalloys, our calculations and phase field simulations have revealed that the perfectly periodic arrangement of cuboidal precipitates is unstable whatever the volume fraction. Specifically, the stability analysis together with the phase field simulations have shown that configurational instabilities lead to formation of defects closely related to experimentally observed branches and herringbone patterns. Also, our phase field simulations have shown that precipitate shape evolutions are correlated to these defects. Thus, during an isothermal annealing, even if the microstructure approximates well on average a cubic arrangement of precipitates, it must necessarily contain defects. To characterize microstructure evolution in detail, it is therefore important to study the dynamics of these defects.

Preliminary observations on 2D phase field calculations of the γ/γ microstructure in Ref. [START_REF] Cottura | Role of elastic inhomogeneity in the development of cuboidal microstructures in Ni-based superalloys[END_REF] have revealed that, during an isothermal heat treatment, branches and macro-dislocations migrate in the microstructure according to mechanisms close to dislocation climbing (see also Ref. [START_REF] Vaithyanathan | Coarsening of ordered intermetallic precipitates with coherency stress[END_REF][START_REF] Bouar | Coarsening of the chessboard-like structures: a TEM study and computer simulations[END_REF]) and gliding mechanisms. In a future work, we intend to extend this analysis by examining the dynamics of these arrangement defects using the recently proposed S-PFM methods, which allows the simulation of large-scale microstructures [START_REF] Finel | A Sharp-Interface Phase Field Method[END_REF].

Figure 1 :

 1 Figure 1: Defects in the cuboidal microstructure: (a) Branch and macro-dislocation patterns in the AM1 superalloy [17]; (b) Chevron patterns in a binary model nickel-based superalloy [18].

Figure 2 :

 2 Figure 2: (a) Eigenvalue spectrum of the curvature tensor for spherical precipitates with fv = 37 %. Longitudinal (L) and transverse (T) branches are indicated, and the dotted lines highlight wave vectors with ζ = 2π/(12a). (b) Sketches of the corresponding deformed arrangements, in which spherical precipitates are represented as blue spheres.

Figure 3 :

 3 Figure 3: (a) Eigenvalue spectrum of the curvature tensor for cubic precipitates with fv = 37%. The transverse branch labelled T 2 between M and Γ corresponds to the eigenvector along [1 10] and the dotted lines highlight wave vectors with ζ = 2π/(12a). (b) Two dimensional sketches of the corresponding deformed arrangements in a (001) plane, for which the displacement vectors lie in the (001) plane.

Figure 4 :

 4 Figure 4: Stability landscapes of the (a) longitudinal L and (b) transverse T 2 (bottom) branches of the periodic arrangement of cubic precipitates, for volume fractions varying from 5% to 95%. The curvatures are scaled using a bi-symmetric log transformation [33] sgn κ 0 α (τ , fv) log 10 1 + |κ 0 α (τ , fv)/c| , with c = 1/250, to capture the wide range of variations.

Figure 5 :

 5 Figure 5: (a) Polar plots of the normalized elastic kernels B(n). (b) Eigenvalue spectra of the curvature tensor for different values of the elastic anisotropy ξ. Precipitates have a cubic shape and the volume fraction is fv = 37%. The spectrum for ξ = 0 is superposed with the axis of abscissas. Also, to help the reading, labels are added to the longitudinal and transverse modes for ξ = -1.5.

Figure 6 :

 6 Figure 6: Evolution of a microstructure initially perturbed along a κ 0 L (ζ, 0, 0) destabilizing mode with ζ = 2π/8a. The volume fraction is fv = 37%, and the snapshots correspond to reduced times (a) t * = 0, (b) t * = 9000, and (c) t * = 18000.

  (ii) The precipitates that have been made farther apart from each other along [100] become thinner along this direction and larger along [010] and [001] to give platelet precipitates. With their periodic replicas, they form rafts perpendicular to the [100] elastic soft direction. (iii) The other precipitates adopt irregular shapes.

Figure 7 :

 7 Figure 7: Evolution of a periodic microstructure initially perturbed along a κ 0 T 2 (ζ, ζ, 0) destabilizing mode with ζ = 2π/4a. The volume fraction is fv = 37%, and the phase field simulation results correspond to the reduced times (a) t * = 0, (b) t * = 12000, and (c) t * = 24000.

Figure 8 :

 8 Figure 8: Eigenvalue spectra of the curvature tensor for periodic arrangements of cubic precipitates at (a) fv = 37% and (b) fv = 60%.The continuous lines are obtained from Eq. (14). Blue circles and red squares correspond to curvatures computed using the FFT solver with ∆C = 0% and ∆C = 50%, respectively.

Figure 9 :

 9 Figure 9: Phase field results for the longitudinal perturbation along [100] with (a) homogeneous elasticity and (b) inhomogeneous elasticity. Top: final microstructures at (a) t * f = 15000, and at (b) t * f = 350000. Middle: Cross sections of the microstructure evolution. Bottom: horizontal profiles of the topological phase ψx(r) at t * = 0 (red) and at different time steps (black).

Figure 10 :

 10 Figure 10: Snapshots (left: 3D image, right: 2D cross section) of cubic macro-lattices perturbed by a transverse T 2 mode along [ζζ0] (ζ = 2π/4a) with (a) homogeneous elasticity and (b) inhomogeneous elasticity.

Note that the branches between Γ and X presented in Fig.3aof Ref.[START_REF] Khachaturyan | Spatially Periodic Distributions of New Phase Inclusions Caused by Elastic Distorsions[END_REF] are questionable, especially for the behavior close to the limits of the first Brillouin zone and for the relative position of the longitudinal and transverse branches.

The existence of numerous instability modes at low volume fraction is consistent with the observation, in low volume fraction γ/γ alloys[START_REF] Ardell | On the modulated structure of aged Ni-Al alloys[END_REF][START_REF] Vaithyanathan | Coarsening of ordered intermetallic precipitates with coherency stress[END_REF], of isolated precipitate alignments that cannot be described as a perturbation of a perfect cubic macrolattice of precipitates.

Acknowledgements

The authors thank Didier Locq and Pierre Caron (ON-ERA) for fruitful discussions and for providing electron micrographs of single-crystal superalloys.