Silvère Akamatsu

Sabine Bottin-Rousseau

Gabriel Faivre

Efim A Brener

Scaling theory of two-phase dendritic growth in undercooled ternary melts

Keywords: numbers: 47.54.+r, 61.72.Mm, 81.10.Aj, 81.30.Fb

Two-phase dendrites are needle-like crystals with a eutectic internal structure growing during solidification of ternary alloys. We present a scaling theory of these objects based on Ivantsov's theory of dendritic growth and Jackson and Hunt's theory of eutectic growth. The additional introduction of the relationship ρ ∼ λ (ρ: dendrite tip radius; λ: eutectic interphase spacing) suggested by recent experimental results [Akamatsu et al., Phys. Rev. Lett., 104, 056101 (2010)] leads to a complete solution of the selection problem and to the scaling rule ρ ∼ v -1/2 (v: dendrite tip growth rate).

The term of two-phase dendrite is used in materials sciences to designate needle-shaped crystals with a fine two-phase internal structure like those observed during univariant solidification of three-component alloys [1,2]. Few detailed studies have so far been devoted to these objects. The first theoretical questions that have to be dealt with are whether two-phase dendrites can grow in a stationary way and, if so, whether anisotropy effects are crucial in the process. An answer to these questions has recently been provided by the finding of stationary spiral two-phase dendrites (for short, sp dendrite) during directional solidification of a transparent nonfaceted ternary eutectic [3,4]. The properties of the sp dendrites are illustrated in Fig. 1. Like one-phase dendrites, they can be subdivided into a tip region characterized by a smooth outer shape, and a tail region where this shape is disrupted by morphological instabilities. In the tip region, the outer shape is a paraboloid. The two-phase substructure is generated by a spiral eutectic pattern located at the tip and rotating at a constant rotation speed. Except close to the tip, an interphase spacing λ, which is approximately equal to the helix pitch, can be defined. More details are given below. For now, we simply stress the fact that the spiral mechanism makes possible a perfectly steady two-phase growth along the curved solid-liquid interface, and that this mechanism in no way involves anisotropy effects. This has recently been confirmed by phase-field numerical simulations of sp dendrites performed in a model ternary system in which the level of interfacial anisotropy could be varied [5].

In this article, we present a scaling theory of the steady state of two-phase dendrites in free growth in a ternary system without surface tension anisotropy. Our main objective is to cast light on the mechanism of selection and the range of existence of this growth structure in terms of morphological and thermodynamical parameters. The paraboloid outer shape and the regular eutectic substructure of sp dendrites suggest to view the diffusion field * Electronic address: akamatsu@insp.jussieu.fr outside two-phase dendrites as being composed of two parts: a long-range field, which only depends on the average composition of the liquid along the envelope of the dendrite, and obeys a theory similar to that of one-phase dendrites; a short-range field, which is driven by the differences in concentration in the liquid in front of the two eutectic phases, and is relevant to a theory similar to that of lamellar eutectic growth. We begin with a brief summary of these theories. Consider the free growth of a one-phase dendrite in a AX binary alloy, where X is the solute. As a first step, the theory asumes that surface-tension effects are negligible. This implies that the solid-liquid interface is an isoconcentration surface, and that the solute mole fractions (for short, concentrations) of the liquid and the solid are those of the liquidus and solidus at the temperature T of the system. Let x and x s be these concentrations, respectively, and x ∞ be the solute concentration far from the dendrite. Let v be the growth speed of the dendrite tip, ρ the tip radius, D x the diffusion coefficient of X in the liquid, and

∆ x = x -x ∞ x -x s (1)
the degree of supersaturation of the liquid. Ivantsov demonstrated that all the paraboloids such that

vρ D x = F Iv (∆ x ) , (2)
where

F Iv (∆ x ) ∼ ∆ x /ln(∆ -1
x ) when ∆ x << 1 (which is the case we consider here), are steady-state solutions of the equations of growth without surface-tension effects [6]. In the next step an additional equation of the form ρ 2 v = constant is provided by introducing slightly anisotropic surface-tension or kinetic effects. As a result, v and ρ are uniquely determined and linked by a ρ ∼ v -1/2 scaling relationship at given ∆ x . This has been established for both 2D [7][8][9][START_REF] Pomeau | Amar in Solids far from equilibrium[END_REF] and 3D dendrite growth [START_REF] Ben Amar | [END_REF]12]. Similar selection and scaling rules for v and ρ have been found for various other needle-like crystal growth shape, such as the doublon (also called double finger) [13,14] and the two-phase finger [15,16]. Unlike the dendrite, these last growth shapes are not crucially dependent on interfacial anisotropy and are not selected in orientation. On the other hand, like the dendrite, their shape does not depart much from an Ivantsov paraboloid, except at the tip, where a specific local structure ensures their stabilization and speed selection.

Consider now a AB binary alloy, where B is the solute, having a eutectic plateau at some temperature T E and a global concentration c ∞ falling well inside this plateau. The growing solid contains two different crystal phases, namely, a A-rich α phase and a B-rich β phase. At small solidification rates the system admits steady states consisting of a planar αβ growth front with a uniform spacing λ and a small average undercooling with respect to T E . We summarize the Zener-Hillert-Jackson-Hunt (JH) theory of these states [17]. This theory was formulated for directional solidification in the G/V → 0 limit and thus also holds in free growth. It is valid under the condition that λV /D << 1, where D is the diffusion coefficient of B in the liquid, which was verified during the experimental observation of sp dendrites. The average concentration c α of the liquid in front of the α-liquid interfaces is larger than the average concentration c β in front of the β-liquid ones at the temperature T of the growth front. This forms the basis of the exchange of solute between solid phases during growth. To a good approximation, c α is on the α liquidus of the alloy and therefore is in equilibrium with a point c α s of the α solidus. The same applies to c β mutatis mutandis. JH showed that the mass conservation equation at the interface leads to

P λV D = ∆ sol , (3)
where

∆ sol = c α -c β c β s -c α s ( 4
)
and P is a small (< 0.1) numerical factor occurring from the summation over the periodic structure. Note, incidentally, that the fraction η of β phase in the solid is related to c ∞ by the global mass conservation equation c ∞ = ηc β s + (1 -η)c α s . Furthermore, capillary effects (namely, the interface equilibrium conditions at the α-βliquid junctions) generate a mean curvature of the solidliquid interface. This gives an additional "capillary" contribution to the average undercooling of the growth front, which reads

d 0 λ = ∆ cap , (5)
where d 0 is a material-dependent capillary length. The spacing value

λ JH = √ P -1 d 0 DV -1 at which ∆ sol = ∆ cap = ∆ E 2 (6)
is an important characteristic length for eutectic growth patterns. To be sure, there is no "strong selection" (in the sense that λ tends towards λ JH , or a value close to λ JH , over time at constant V ) of the spacing in binary eutectic growth, contrary to what was long believed. However, it is also true that any local structure containing λ values deviating from λ JH by a factor of more than about 2 is short-lived so that, broadly speaking, λ is almost always close to λ JH , a fact which is sometimes referred to as "weak selection" [18][19][20].

Consider finally the directional solidification of a αβ solid in a ABX ternary alloy, as studied in Ref. [3]. The rejection of X into the liquid during growth generated a concentration gradient that caused a large-scale (compared to λ) Mullins-Sekerka-like instability of the αβ growth front at values of V larger than a threshold value V c [21,22]. Spiral two-phase dendrites were observed at V >> V c and appeared in the form of isolated objects. Their tip region was thus growing under nearly free-growth (virtually zero G/V and infinite primary spacing) conditions [23][24][25]. The theory presented herein should therefore be applicable to them. Their growth direction, like those of doublons and two-phase fingers, was not fixed, but history dependent, meaning that different sp dendrites had different values of their tip growth rate v at fixed V . The explored v-range was 0.1 -1 µms -1 . The tip radius ρ of the sp dendrites was measured by fitting a parabola to the contour of the images. The extension of the tip region was between 5 and 10 ρ. The spacing λ was determined as the product of v and the period of rotation of the spiral pattern. The measured λ values turned out to be within 10% of the λ JH ones for the AB binary alloy (λ 2 JH V ≈ 10.2 µm 3 s -1 ) over the whole experimental v range. The ρ values showed a similar variation, with ρ ≈ 0.75λ JH in average. Most importantly, two sp dendrites growing simultaneously, side by side, had the same values of λ, ρ and tip temperature. There is thus a clear indication of a strong selection of both the tip radius and the eutectic spacing of sp dendrites according to a λ ∼ ρ ∼ v -1/2 law.

We now come to the theory of two-phase dendrite growth. A letter c will designate a concentration of B, a letter x a concentration of X and a pair (c, x) a composition of the alloy. Concentrations in the solid phases will be tagged by a subscript s and average concentrations along the growth front by a bar. As a first step, we consider the long-range concentration field of a sp dendrite. For ease of exposition, we assume that there is no cross-diffusion in the liquid and that the diffusion coefficients of B and X, called D and D x , respectively, have the same order of magnitude, as is in fact often the case in the experiments. Then an Ivantsov equation holds for each component separately: we have

vρ D = F Iv (∆ c ) , vρ D x = F Iv (∆ x ) , (7)
where

∆ c ∼ c -c ∞ c -cs ; ∆ x ∼ x -x ∞ x -xs . (8)
The average concentration x is defined by the equation

x = ηx β + (1 -η)x α , (9)
where η is unknown. Similar equations hold for c, cs and xs . By elimination of vρ, Eqs. (7) lead to a relationship between ∆ c and ∆ x , and thus between the average concentrations. Knowing that D x ≈ D, we rewrite Eqs. (7) in the form

∆ c ∼ ∆ x , (10)
and

vρ D ∼ ∆ x ln ∆ -1 x . ( 11
)
Regarding the short-range part of the concentration field, we approximate the spiral substructure near the dendrite tip by some effective lamellar structure. The resulting uncertainty in the definition of λ and v is within the limits of the scaling theory proposed here. Then, according to Eqs. ( 3) and ( 5), and assuming that λ/λ JH is kept constant at a value close to 1 in agreement with the experimental measurements:

P λv D ∼ ∆ E , (12)
and

d 0 λ ∼ ∆ E , (13)
where ∆ E is given by Eqs. ( 4) and ( 6). There are apparently twelve unknown concentration variables in Eqs. ( 10) to ( 13), but their actual number is only two due to the four relationships defining the average concentrations and the six relationships giving, for each phase, three equilibrium concentrations as a function of the fourth one. Thus, at this stage there are six unknowns, say, x α , x β , η, v, ρ and λ, at fixed (c ∞ , x ∞ ), and four relationships between them. We now look for additional relationships. For the sake of clarity, we begin by considering the particular case of a system having A↔B symmetry (Fig. 3a). In other words, we assume that the isothermal section of the phase diagram is symmetrical with respect to the x = 1 -2c axis and that (c ∞ , x ∞ ) is on this axis. Incidentally, these conditions were approximately fulfilled in the experimental system in which sp dendrites were observed. Under these conditions, the sp dendrites must satisfy the symmetry requirements η = 1/2 and x α = x β . The number of unknowns is now four, say, x, v, ρ, λ. On the other hand, the A↔B symmetry also imposes that x = 1 -2c and xs = 1 -2c s , so that Eq. ( 10) becomes an identity and is no longer part of the equations of the problem. One additional relationship arising from the interaction between the long-range and the short-range dynamics must be found. The experiments clearly suggest that there is a linear, or almost linear, relationship between λ and ρ. We therefore write:

λ ∼ ρ . (14)
Given that λ ∼ λ JH , Eq. ( 14) also implies that λ ∼ ρ ∼ v -1/2 . Presumably, Eq. ( 14) arises from the fact that sp dendrites cannot grow with ρ << λ for geometrical reasons while, on the other hand, sp dendrites with ρ >> λ are unstable with respect to a decrease in ρ that would make them grow faster. So, the only ρ vs λ-range in which sp dendrites can stabilize is the one given by Eq. ( 14). From Eqs. ( 11), ( 12) and ( 14) one finally obtains

∆ E ∼ P ∆ x ln ∆ -1 x . ( 15
)
Both the quantities P and 1/ln(∆ -1 x ) are substantially smaller than unity. Thus ∆ E /∆ x << 1 and x u -x << 1. We can therefore approximate ∆ x as

∆ x ≈ x u -x ∞ x u -x su . ( 16
)
Plugging this relationship into (15), it can be seen that ∆ x and ∆ E , i.e. the driving forces for growth, now depend only on control and material parameters. Expending the equilibrium equations to the first order in x -x u near the univariant points, one obtains from Eq. ( 8)

x u -x ∼ c β su -c α su 2|∂c α /∂x u | -1 ∆ E , (17)
which yields x as a function of the control and material parameters. Then one gets from Eqs. (12) ( 13) and ( 14)

v ∼ D P d 0 ∆ 2 E ( 18
)
and

λ ∼ ρ ∼ d 0 ∆ E , (19)
which completely solves the selection problem. As a last remark, we put Eq. ( 17) into the form

x u -x ∼ K (x u -x ∞ ) . (20)
It can be seen from Eqs. ( 15), ( 16) and ( 17) that K is the product of P/ln ∆ -1

x and some material parameters that generally are on the order of unity. Thus, in general, K << 1 and x is much closer to x u than x ∞ .

In the case of systems without any particular symmetry, there are two more independent unknowns (say, η and x α -x β ) and one more relationship than in symmetrical systems. One additional relationship is thus required for determining all the unknowns uniquely. However, this relationship is about the details of the short-range field and is of little consequence as regards the domain of existence of two-phase dendrites. As an exploratory suggestion, we therefore set x α -x β to zero. Keeping the same symmetrical phase diagram as above, we introduce a small departure δc ∞ of the alloy concentration from the axis of symmetry and calculate the departures of η, c and cs from their values at δc ∞ = 0. Expending Eqs. ( 8) and [START_REF] Pomeau | Amar in Solids far from equilibrium[END_REF] to first order and using the fact that ∆

E << ∆ x , one gets δc ∼ -δc ∞ ∆ E /∆ x , δc s ∼ δc ∞ /∆ x and δη ∼ δc ∞ (c β su -c α su )∆ x . (21)
It is interesting to note that δη, δc s and δc ∞ have the same sign and that |δc s | is substantially larger than |δc ∞ |. When, for instance, δc ∞ is positive, the dendrite is, on average, more enriched in B than the liquid. The concentration range of existence of sp dendrites, as determined by Eq. ( 21) and the condition |δη| < 0.5, is delimited by the conjugation lines (c σ , x σ )-(c σ s , x σ s ). Thus, as x is close to x u , the range of existence of sp dendrites roughly coincides with the three phase domain (Fig. 3b). The operating point of the dendrite tip remains close to the univariant groove of the phase diagram as the alloy composition goes through this domain.

In conclusion, we have obtained a complete solution of the selection problem for two-phase dendrites on the condition of postulating a relationship linking the tip radius of the dendrite to the interphase spacing of the eutectic microstructure. This relationship reads ρ ∼ λ and is based on the geometry of the source of the eutectic patterns located at the tip of the dendrite. This very natural length scale selection by the underlying eutectic structure replaces the much more delicate surface tension anisotropy mechanism that takes place in one-phase dendrites (the so-called microscopic solvability criterion). No particular condition appears to be required for the alloy to display two-phase dendrites beyond the fact that its composition is close to a univariant groove of the phase diagram and that it solidifies in a nonfaceted way. This suggests that sp dendrites should be of common occurrence in univariantly solidified ternary eutectics.

FIG. 1 :

 1 FIG. 1: Spiral two-phase dendrites. a and b: In situ optical micrographs taken during directional solidfication of a transparent ternary alloy (for experimental details, see [3]). a) bulk sample. Horizontal dimension: 105 µm. b) Thin samples. Horizontal dimension: 65 µm. c: Schematic 3D representation of the solid-liquid interface. One of the solid phase is grey, the other white. (d) Longitudinal section of the internal microstructure (dotted line: parabolic envelope).

FIG. 2 :

 2 FIG.2: Isothermal section of the phase diagram of the ABX alloy in Gibbs' triangular coordinates. Grey areas: two-phase domains. Triangle: three-phase domain.

5 FIG. 3 :

 53 FIG. 3: Isothermal section of a phase diagram having A↔B symmetry. a) Symmetrical system. b) Symmetry-broken system. Closed circles: (c∞, x∞). Closed squares: (c, x). Closed triangles: (cs, xs). Open squares: (c α , x α ) and (c β , x β ). Open triangles: (c α s , x α s ) and (c β s , x β s ).

Acknowledgments

We thank M. Plapp and G. Boussinot for fruitful discussions. EAB benefited from an invited-professor position at UPMC.


Supriyo Ghosh

Abhik Choudhury

Mathis Plapp

Sabine Bottin-Rousseau

Gabriel Faivre

Silvère Akamatsu

Interphase anisotropy effects on lamellar eutectics: A numerical study

In directional solidification of binary eutectics, it is often observed that two-phase lamellar growth patterns grow tilted with respect to the direction z of the imposed temperature gradient. This crystallographic effect depends on the orientation of the two crystal phases α and β with respect to z. Recently, an approximate theory was formulated that predicts the lamellar tilt angle as a function of the anisotropy of the free energy of the solid(α)-solid(β) interphase boundary. We use two different numerical methods -phase-field (PF) and dynamic boundary-integral (BI) -to simulate the growth of steady periodic patterns in two dimensions as a function of the angle θ R between z and a reference crystallographic axis for a fixed relative orientation of α and β crystals, that is, for a given anisotropy function (Wulff plot) of the interphase boundary. For Wulff plots without unstable interphase-boundary orientations, the two simulation methods are in excellent agreement with each other, and confirm the general validity of the previously proposed theory.

In addition, a crystallographic "locking" of the lamellae onto a facet plane is well reproduced in the simulations. When unstable orientations are present in the Wulff plot, it is expected that two distinct values of the tilt angle can appear for the same crystal orientation over a finite θ R range. This bistable behavior, which has been observed experimentally, is well reproduced by BI simulations, but not by the PF model. Possible reasons for this discrepancy are discussed.

I. INTRODUCTION

The solidification of nonfaceted binary alloys of near-eutectic composition produces selforganized composite materials, the microstructure of which results from the coupled growth of two distinct crystal phases from the melt. In directional solidification -crystal growth at constant velocity V in a fixed thermal gradient G -the two most frequent morphologies are parallel lamellae (platelets) of the two phases and fibers (rods) of one phase dispersed in a continuous matrix of the other phase. They are a trace left behind in the solid by the dynamics of self-organized two-phase patterns at the solidification front. The diffusion-controlled growth theory of these patterns is well established [1,2]. More complex microstructures have also been observed in experiments [3,4], and further explored with the help of numerical simulations using boundary-integral (BI) [5,6] and phase-field (PF) methods [7][8][9][10][11][12], with good agreement between experimental and numerical results.

There are, however, major experimental observations that cannot be accounted for by the existing theories. Most importantly, it has been known for a long time that the solidification dynamics of eutectic patterns may strongly depend on the crystal orientation of the two solid phases [13]. Such crystallographic effects during eutectic growth have been neglected so far, both in theories and models, in spite of their practical importance [14,15]. We will focus here on lamellar eutectics. In previous works [16,17], a distinction has been made between two types of eutectic grains. A eutectic grain is defined as a region of substantially uniform crystal orientation of the two solid phases α and β, and thus with a constant orientation relationship between the two crystals. On the one hand, in floating (eutectic) grains, the dynamics of the lamellar solidification front patterns is well described by the standard theory of regular eutectics. In particular, spatial inhomogeneities of the lamellar spacing are smoothed out with time by a "spacing-diffusion" process [8,9] generically present in out-of-equilibrium pattern forming systems [START_REF] Manneville | Dissipative Structures and Weak Turbulence[END_REF]. On the other hand, in locked grains, eutectic lamellae grow in a direction that is essentially aligned to a certain crystallographic plane, and are inclined (or tilted) with respect to the main growth axis z [17,[START_REF] Hecht | [END_REF]. The strength of this crystallographic locking effect varies between different eutectic grains. In strongly locked grains, the spacingdiffusion process is absent [16]. Examples of floating and locked eutectic patterns observed in thin-sample directional solidification experiments are shown in Fig. 1. The corresponding schematic views of the interfaces are depicted in Fig. 2. In the locked case, tilted lamellar microstructures are left behind in the solid by a eutectic growth front pattern that drifts laterally at a constant velocity V d (the pattern is in a steady state in the traveling reference frame). The magnitude of the drift velocity can be large, i.e. comparable to the pulling velocity, which leads to the freezing of strongly tilted lamellae, as shown in Fig. 1b. In a steady-state condition, the (lamellar) tilt angle θ t is defined by tan θ t = V d /V . On the basis of in situ directional solidification observations using thin samples of metallic and transparent organic eutectic alloys, a conjecture was formulated recently that permits to relate the value of θ t to the anisotropy of the free energy of the interphase boundaries (interfacial anisotropy) [20,21].

The main underlying hypotheses are that (i) only the solid-solid interfaces are anisotropic (i.e., in a nonfaceted alloy, the anisotropy of the solid-liquid interfaces has a negligible effect on the lamellar growth dynamics), and (ii) the solid-liquid interface keeps virtually the same shape -with mirror symmetry about the mid-plane of a lamella -as for standard (nontilted) lamellae. The tilted pattern shown in Fig. 1b satisfies these conditions. In essence, under this symmetric-pattern (SP) approximation, the conjectured theory states that the Cahn-Hoffman surface tension vector σ (to be defined below) is aligned with z (Fig. 2b). (SP) approximation. α, β: solid phases. L: liquid. z: growth direction parallel to the thermal gradient. x: direction of the isotherms. θ sp : SP-approximation lamellar tilt angle. The lateral drift velocity is given by tan θ sp = V d /V , with V the pulling velocity. Other symbols: see text.

The goal of the present work is to test the SP approximation by numerical simulations, in which the anisotropy of the interphase boundaries can be freely chosen. We use two different numerical models, in two dimensions. The first one is a sharp interface code, which uses the dynamic boundary-integral (BI) formalism previously developed by Karma and Sarkissian [6]. This method combines an evaluation of the solute diffusion field by the boundary-integral method with an explicit front-tracking algorithm for the interface evolution. For simplicity, the solid-liquid interfaces are assumed to be isotropic. The anisotropy of interphase boundaries is incorporated by changing the local equilibrium condition at the trijunctions from the Young to the Young-Herring law. The second method is a multi-phase-field model, based on recent grand-canonical formulations of alloy solidification [22,23]. The interphase anisotropy is directly introduced into the underlying free-energy functional.

We characterize a given eutectic grain by the anisotropy function, or Wulff plot (i.e.

the polar plot of the surface free energy) of the interphase boundary. It is important to note that the details of this function depend solely on the relative orientation of the α and β crystals. In two dimensions, changing the in-plane orientation of a eutectic grain with respect to the growth axis z while keeping the relative αβ orientation fixed is equivalent to rotate the Wulff plot globally by a single angle θ R . Our results can be classified into two main categories, according to unstable orientations being absent or present in the Wulff shape, respectively. For all the anisotropy functions that we have tested and that do not exhibit unstable orientations, the two numerical methods give almost identical results for the tilt angle θ t as a function of the rotation angle θ R . Importantly, the value of θ t follows quite closely the variation of θ sp vs θ R , where θ sp is the tilt angle predicted by the SPapproximation. In addition, a strong locking of the lamellae onto a certain direction is well reproduced in both BI and PF simulations by using an anisotropy function with a peaked minimum. Moreover, the tilt angle is found to be largely independent of the lamellar spacing and the pulling velocity, in good agreement with the SP approximation. For anisotropies that are large enough to create orientations that exhibit a Herring instability [24], the two numerical methods give different results. In the BI simulations, a phenomenon of bistability is observed, that is, for a finite range of eutectic-grain orientations, there are two stable lamellar patterns with different tilt angles. Those two branches of steady-state solutions are essentially the same as the ones predicted by the SP approximation. Manifestations of a bistable behavior have indeed been observed in thin-sample solidification experiments [20].

In contrast, no such bistability is observed in the phase-field simulations.

The remainder of the article is structured as follows. In Sec. II, we recall some facts about eutectic grains and anisotropic interfaces, present the equations of the eutectic growth problem in the sharp-interface formulation and review the SP-conjecture for the prediction of the tilt angle. In Sec. III, we briefly outline our simulation methods and specify how we have incorporated the interface anisotropy in the existing models. In Sec. IV, we describe our results for various choices of the anisotropy function. In Sec. V, we will discuss separately (i) a way to estimate the accuracy of measurements of the interphase boundary Wulff plot using experimental observations with the so-called rotating directional-solidification method [21], and (ii) possible reasons for the absence of hysteretic behavior in the high-anisotropy case in the PF simulations. Conclusions and perspectives are presented in Sec. VI.

II. BACKGROUND A. Sharp-interface problem

Consider the directional solidification of a binary eutectic alloy in an externally imposed temperature gradient G along the z axis with a fixed pulling velocity V . We make the approximation that the temperature field is independent of the solid-liquid front shape (frozentemperature approximation), which is a good approximation for thin-sample solidification in which heat conduction takes mainly place in the sample walls. In a two-dimensional system (an appropriate model for thin-sample directional solidification) in which the temperature gradient is directed along the z axis, the temperature field is hence given by

T ( x, t) = T E + G(z -V t), (1)
where we have chosen the origin of the z axis at the eutectic temperature T E at time t = 0, and x is a position vector in the (x, z) plane.

A dimensionless concentration field u is introduced,

u( x, t) = C( x, t) -C E C β -C α , (2)
where C( x, t) is the space-and time-dependent composition of the alloy, and C E , C α , and C β are the equilibrium compositions of the three phases (liquid, solid α and solid β) that are in coexistence at T = T E . In the one-sided model of solidification, this field obeys the diffusion equation in the liquid with the diffusivity D,

∂ t u = D ∇ 2 u, (3)
whereas no diffusion takes place in the solid. At the solid-liquid interfaces, the conservation of solute implies

V n u ν = -Dn • ∇u, (4)
where V n is the normal growth velocity of the interface,

u ν = (C ν -C E )/(C β -C α ) for ν =
α, β, and n is the unit normal vector to the interface pointing into the liquid. In writing down this expression, we have made the simplifying assumption that the concentration differences between the phases do not depend on temperature (parallel liquidus and solidus lines).

For isotropic interfaces, the classic problem of eutectic growth is completed by the local equilibrium condition,

u int =    -(ζ(x) -V t)/l α T -d α 0 κ, α-L interface (ζ(x) -V t)/l β T + d β 0 κ, β-L interface (5)
Here, ζ(x) indicates the position of the solid-liquid interface in the z direction,

l ν T = G |m ν |(C β -C α ) (6)
is the thermal length of phase ν, with m ν being the liquidus slope at the eutectic point,

d ν 0 = γ νL T E L ν |m ν |(C β -C α ) (7)
is the capillary length of phase ν, with L ν the latent heat of melting per unit volume and γ νL the surface free energy, and κ is the interface curvature, counted positive for a convex solid.

At the trijunction point, local equilibrium implies the balance of surface tensions (Young's law). For isotropic interfaces,

γ αL tαL + γ βL tβL + γ αβ tαβ = 0, (8)
where tµν are the unit vectors tangent to the µ-ν interface at the trijunction point, and pointing away from the trijunction.

B. Anisotropic interphase boundaries

A lamellar-eutectic solid consists of eutectic grains made of crystals of the solid phases α and β with uniform orientations. A (eutectic) grain hence constitutes a heterophase bicrystal.

The relative orientation between the lattices of the two phases (constant within a grain) determines the interphase boundary energy and its anisotropy, which may therefore vary between grains. For a description of interphase configurations in directional solidification, we need to specify the orientation of the bicrystal and that of the interphase boundary with respect to the temperature gradient and the sample plane of the directional solidification setup. In two dimensions, two angles with respect to the temperature gradient axis are sufficient for a complete specification of these orientations.

Consider first a bicrystal in a fixed orientation with respect to the temperature gradient.

Let n be the unit normal vector of the interphase boundary, and θ the angle between n and the x axis (we have n x = cos θ and n z = sin θ). Furthermore, let the anisotropic interphase boundary energy be given by

γ αβ (θ) = γαβ a c (θ), (9)
where γαβ is a constant and a c (θ) is a dimensionless function.

It may be useful for the following to recall a few standard definitions and well-known facts. The vector tangent to the interface is given by t = -dn(θ)/dθ. With this definition, the angle between the interface direction and the z axis is also equal to θ (as depicted in Fig. 2b). The Cahn-Hoffman ξ and σ vectors [25] are defined by

ξ = γ αβ n -γ αβ t (10)
and

σ = γ αβ t + γ αβ n, (11)
where γ αβ = dγ αβ (θ)/dθ. With their help, the equilibrium shape of a β inclusion inside an α matrix and the anisotropic equilibrium condition at trijunction points can be obtained in a simple way. The Wulff plot is defined by r(θ) = γ αβ (θ)n. The minimum-energy shape (Wulff shape) is traced by plotting the vector ξ(θ). For low anisotropies that satisfy γ αβ (θ)+γ αβ (θ) > 0 for all orientations, the Wulff shape is smooth. If the interface stiffness γ αβ (θ) + γ αβ (θ) becomes negative for a given value of θ, a flat interface of this orientation is unstable with respect to the formation of a hill-and-valley structure (Herring instability [24]). In this case, the plot of ξ as a function of θ has self-intersections, and the Wulff shape is given by the inner convex part only. The other parts, often called "ears", consist of three segments delimited by turning points (see for example Ref. [26] for details and illustrations). The interface stiffness is negative only on the middle segment, but all orientations located on the "ears" are missing from the physically observable convex equilibrium shape.

When the bicrystal is rotated with respect to its reference configuration by an angle θ R , the interphase energy becomes

γ αβ (θ) = γαβ a c (θ -θ R ). (12)
We choose the reference configuration (θ R = 0) such that an interphase orientation of minimal energy is aligned with the growth direction. Note that positive and negative θ R thus correspond to rotations to the left and to the right with respect to a minimum-energy direction. Below, we will essentially test two types of anisotropy functions. The first one is of the standard form used for the modeling of crystals with m-fold symmetry, namely,

a m (θ) = 1 -m cos mθ , ( 13
)
where m is the m-fold anisotropy coefficient. It should be noted that for interphase boundaries, at which two centrosymmetric crystals meet, a 2-fold anisotropy is always expected (since two opposite orientation vectors describe the same surface). The second type of anisotropy function that we will study is motivated by the observation of strong locking of growth directions onto certain crystallographic planes. This usually occurs when the two solid phases exhibit an epitaxial orientation relationship [START_REF] Hecht | [END_REF]. In this case, the locking planes correspond to sharp cusp-like minima in the γ αβ (θ) function. In order to avoid the additional difficulties related to the regularization of the cusp (see for example [27]), we prefer to use a deep and narrow, but smooth minimum, which we model by a Gaussian in the anisotropy function, that is,

a g (θ) = 1 -g exp -(θ/w g ) 2 , ( 14
)
where g is the amplitude, and w g is the width of the Gaussian. For both types of anisotropy, a finite range of forbidden orientations appears for large enough values of the anisotropy coefficient. In our simulations, we also used linear combinations of the functions a m (with m = 2, 4) and a g ,

a c (θ) = 1 -g exp -(θ/w g ) 2 -2 cos 2θ -4 cos 4θ (15)
with various values of the parameters m , g and w g . This form can reproduce well typical anisotropies that have been obtained from experiments [21] and molecular dynamics simulations [28].

Let us comment on how the sharp-interface problem introduced previously needs to be modified in order to take the interphase boundary anisotropy into account. Since we suppose that the solid-liquid interfaces remain isotropic (which should be a good approximation for nonfaceted substances), the Gibbs-Thomson conditions, Eq. ( 5), are unchanged. In addition, since we still assume that there is no diffusion in the solid, we do not need to write a local-equilibrium condition along the interphase boundaries. Therefore, the only change that intervenes in the equations is a modification of the local-equilibrium condition at the trijunction, which becomes a Young-Herring equation, that is,

γ αL tαL + γ βL tβL + σ = 0, (16)
with σ given by Eq. (11). Note that the σ vector is not parallel to the interphase boundary (Fig. 3a).

γ αL γ βL α β γ αβ  t  σ ! γ  n z" a) γ αL γ βL α β  σ δ z" θ t" b)
FIG. 3: a) Illustration of the Young-Herring equilibrium condition at the trijunction. n and t: normal and tangent unit vectors of the interphase boundary. b) Definition of the tilt angle θ t and the angle δ between σ and z.

C. Theoretical predictions for the tilt angle

The SP approximation that allows us to obtain a prediction for the growth angle is based on the experimental findings presented in Refs. [20,21]. It was found that, even for lamellae that grow at a large angle with respect to the temperature gradient, the shape of the solidliquid interface is close to the one observed for well-aligned lamellae. More precisely, the "heads" of the lamellae are approximately mirror-symmetric with respect to the mid-plane of the lamellae, which means that the contact angles of the solid-liquid interfaces at the trijunctions are nearly the same on both sides of a lamella (to within the resolution of the experiments). Considering the Young-Herring condition of Eq. ( 16), a strictly symmetric shape is possible only if the vector σ is aligned with the z axis. Using the fact that tαβ = (sin θ, -cos θ) and nαβ = (cosθ, sin θ) in the (x, z) plane, the condition that the x component of σ is zero writes

γ αβ (θ -θ R ) sin θ + γ αβ (θ -θ R ) cos θ = 0. ( 17
)
For a fixed orientation θ R of the eutectic grain, this is a nonlinear equation for the interface orientation θ, which can easily be solved numerically for arbitrary anisotropy functions a c (θ).

As long as the interface stiffness γ αβ (θ) + γ αβ (θ) is positive for all angles, this equation has a unique solution. For negative stiffness, there are ranges of θ R for which there exist three solutions, of which one corresponds to an orientation that is present on the equilibrium shape, one to an unstable orientation, and the third to a metastable orientation with positive stiffness that is missing on the equilibrium shape (belonging to an "ear"). This is the prediction against which we will compare our numerical results. As mentioned, we will note θ sp the value of θ, solution of Eq. ( 17), which is predicted by the SP approximation.

As will be seen below, a steady state with tilted lamellae obtained from our simulations generally does not follow exactly the prediction of Eq. ( 17). The departure from the conjecture can be quantified by the value of the angle between the σ vector and the z axis, which we will denote by δ in the following, as illustrated in Fig. 3b. The SP approximation predicts δ = 0.

III. METHODS

A. Boundary-integral method

For our BI calculations, we have adapted the method developed by Karma and Sarkissian [6]. In the quasistationary approximation, valid for slow growth velocities, the solute diffusion equation, Eq. ( 3) can be replaced by the Laplace equation, ∇ 2 u = 0. The use of Green's function techniques then permits to transform this partial differential equation together with the boundary conditions at the interface into a single integro-differential equation along the solid-liquid front. The numerical procedure to calculate the time evolution of the growth front is as follows. The interfaces are discretized with the help of marker points that are uniformly spaced along the interface. The gradient of the concentration field u, and thus the interface velocity from Eq. ( 4), is obtained from the boundary-integral equation for each interface point except for the trijunction points. After having moved forward the interface points, the new position of the trijunction point is found by solving Eq. ( 8), with the new positions of the first point on each interface taken as input, and the position of the trijunction as unknown. This equation is solved by a relaxation scheme. More details can be found in Ref. [6].

Since, under our hypotheses, the solid-liquid interfaces remain isotropic, the only change that is necessary to incorporate the solid-solid interfacial anisotropy in this model is in the calculation of the trijunction positions. We have replaced Eq. ( 8) by Eq. ( 16); in other words, we have replaced the vector γ αβ tαβ by the vector σ in the routine that calculates the new trijunction positions.

B. Phase-field model

We have used a grand-canonical multi-phase-field model with a multi-obstacle potential, developed and validated recently by Choudhury and Nestler [23]. In this model, accurate simulation results can be obtained with the help of a quantitative thin-interface analysis and an antitrapping current adapted to the obstacle potential. For the solidification of a binary eutectic alloy, we will work with N = 3 phase fields (denoted by φ l , φ α and φ β for liquid, α solid and β solid, respectively) that obey the sum constraint φ l + φ α + φ β = 1, and a diffusion potential µ that is conjugate to the dimensionless concentration field u.

The starting point of the model is the (grand-canonical) free-energy functional

Ω = V ω int (φ, ∇φ) + ω b (T, µ, φ), (18)
where ω int and ω b are the contributions of interfaces and bulk, respectively. The interface part is given by

ω int = a(φ, ∇φ) + 1 W(φ), (19)
where is proportional to the numerical interface thickness. For isotropic interfaces, the gradient energy a(φ, ∇φ) and the multi-obstacle potential W(φ) are given by

a(φ, ∇φ) = N i<j γ ij | q ij | 2 (20)
with q ij = φ i ∇φ j -φ j ∇φ i (a vector that is normal to the ij interface), and

W(φ) =      16 π 2 N i<j γ ij φ i φ j + i<j<k γ ijk φ i φ j φ k if φ ∈ Σ ∞ elsewhere. ( 21
)
The simplex Σ is bounded by φ i ≥ 0 ∀ i and N i=1 φ ν ≤ 1. The first sums in Eqs. ( 20) and ( 21) run over all pairs of distinct phases, with γ ij denoting the isotropic ij surface energy; γ ijk is a third order potential term which prevents the appearance of any "foreign" phases in the binary interfaces.

To define the bulk part of the functional, we start from the Helmholtz free energies of each phase, which we approximate by parabolas as in Ref. [10],

f i (u, T ) = A i u 2 + B i (T )u + E i (T ) (i = α, β, l). ( 22
)
We define the Legendre transforms

ω i (µ, T ) = f i (u, T ) -µu µ = ∂f i /∂u, (23)
which are functions of µ and T . The bulk part of the energy functional is then given by

ω b (µ, T, φ) = N i=1 ω i (µ, T )h i (φ), (24)
where h i (φ) = φ 2 i (3 -2φ i ) + 2 φ i φ j φ k are weight functions which interpolate between phases and satisfy N i=1 h i = 1. For our simulations, we have chosen A α = A β = A l ≡ A, which leads to parallel liquidus and solidus lines [10]. With this choice, the capillary lengths for the α-liquid and β-liquid interfaces are given by

d ν 0 = γ νL ∂ 2 f l /∂u 2 = γ νL 2A . ( 25
)
The equivalence between the two definitions of d ν 0 , Eqs. ( 7) and (25), can be established with the help of a Clausius-Clapeyron relationship for the solid-liquid coexistence line and the definition of u, Eq. ( 2).

Surface tension anisotropy can be incorporated directly into the interface part of this functional. Since both the gradient and potential parts, Eqs. ( 20) and (21), are sums over the different interfaces, it is straightforward to control the anisotropy of each interface independently. Here, we have modified only the contribution of the α-β interface by multiplying the corresponding gradient and/or potential terms with the function a c (θ). We have tested different ways to implement this anisotropy. The technical details, as well as the equa-tions of motion of the model, which are derived from the functional following the lines of Refs. [22,23], can be found in Appendix A.

C. Parameters

We have chosen as a convenient test case a model eutectic alloy with symmetric phase diagram and properties, that is, m α = -m β , u β = -u α , and γ αL = γ βL . In the following, we will therefore drop the phase indices of all parameters (capillary and thermal lengths) for simplicity. We also choose the average value of the solid-solid interface free energy to be equal to the solid-liquid one (γ αβ = γ αL ), which would yield trijunction angles of 120 • for isotropic interfaces. We will specify the simulation parameters using dimensionless ratios, normalizing all lengths with the diffusion length

l D = D V . (26)
We work at the eutectic composition and in the limit of slow velocities and low temperature gradients, which implies d 0 /l D 1, d 0 /l T 1, and the Péclet number

P e = λ/l D 1,
where λ is the lamellar spacing. It may be useful to mention that the value of the Jackson-Hunt minimum undercooling spacing λ m under these conditions is given by

λ m ≈ 3.844 d 0 l D . (27)
As we shall see below, in the regime that we have investigated, the results depend very little on the detailed choice of the parameters.

IV. RESULTS

A. General remarks

Our standard simulation procedure is as follows: we start from a pre-existing pair of lamellae (usually taken from a previous steady-state calculation) and slightly change the rotation angle θ R at a predetermined time step. Immediately after the change in θ R (within a solidification distance of less than one lamellar spacing), the lamellae select a new tilt angle (Figs. 4 and5). With the model Wulff plots that we have implemented, and in the range of control parameters that we have used, the calculations always converged toward a steady-state pattern. We let the system evolve for typically 3-10 lamellar spacings before we measure the steady-state tilt angle θ t of the interphase boundary. In BI simulations, the inclination of the interphase boundary (and those of the tangents to the solid-liquid interfaces at the trijunctions, as well as the σ vector) is automatically given as output data of the simulation. In PF simulations, θ t is measured at a distance from the trijunction of several times the interface thickness, as imposed by the diffuse nature of the trijunction. It was, however, not possible, with our PF method, to extract numerically the contact angles in the moving anisotropic trijunction with a reasonable accuracy. The same is true for the angle δ that quantifies the departure from the SP approximation. between the shape of the anisotropy function and the tilt angles will be detailed below. For the time being, let us comment on some general features of tilted lamellar growth patterns.

First, the selection of a new steady-state tilt angle after a change of the rotation angle θ R occurs very rapidly (over a solidification distance of less than λ), so that the interphase boundary left behind in the solid exhibits a sharp bend, see Figs. 4c and5b. Second, in tilted patterns, the shape ζ(x) of the solid-liquid interfaces most often looks "flatter" than in the fully isotropic case (Fig. 4a,b), in agreement with experimental observations. This results from the fact that the modulus of σ is smaller than the isotropic reference value of the surface free energy, and hence the dihedral angle between the two solid-liquid interfaces at a trijunction is larger than in the isotropic case. Third, the solid-liquid interfaces exhibit a certain degree of asymmetry, a quantitative measure of which can be given by the ratio

[ζ(x) -ζ(-x)
]/ζ 0 , with ζ 0 the amplitude of the solid-liquid cap (and x being centered in the middle of a lamella). In the case of Fig. 4b, that quantity is of about 5%. This is equivalent to stating that (as shown by the BI simulations) the value of δ is not zero. The sign of δ is opposite to that of θ t , which indicates that the diffusion field tends to oppose the lateral drift of the lamellar pattern (as could be expected given the stability of non-tilted lamellar pattern for λ close to λ m in the absence of anisotropy).

Before studying the dynamics of tilted-lamellar solidification patterns by varying the characteristics of the αβ Wulff plot, we have tested the influence of two control parameters, namely, the lamellar spacing λ, and the Péclet number P e, on the tilt angle, for a given crystallographic configuration. Figure 6 shows the variation of θ t as a function of λ/λ m for given values of the other control parameters and for a fixed Wulff plot (BI simulations). In the scanned λ/λ m range, the variations of θ t are very small (a fraction of a degree). Moreover, the angle δ, which characterizes the departure from the SP approximation remains smaller than 1 degree (in absolute value). This indicates that the steady-state tilt angle θ t is smaller than, but close to the SP-approximation value θ sp . It can be seen that for small lamellar spacings, δ tends towards zero and θ t towards θ sp with decreasing λ. The Jackson-Hunt analysis predicts that for λ = λ m , the contributions of the interface curvature and the diffusion field to the front undercooling are of equal magnitude, whereas for smaller spacings capillarity dominates over the composition variations due to diffusion. Since the SP approximation takes into account only capillary phenomena, it is not surprising to see that it becomes more accurate in the small-spacing limit. It should also be mentioned that those small spacings are subject to a lamella-elimination instability for isotropic systems [9]. The value λ ≈ 1.75λ m at which θ t passes through a minimum (and -δ through a maximum) more or less corresponds to a change of the shape of the solid-liquid interface, which, for the larger λ values, exhibits a concave part in the center of each lamella (it may also be noted that, for isotropic interfaces, the lamellar pattern undergoes an oscillatory instability for λ ≈ 2.2λ m [6]).

Figure 7 shows the variation of θ t with the Péclet number P e, calculated with the PF method by varying the diffusion coefficient D in the liquid (this is equivalent to varying the velocity V ), at fixed λ, d 0 , and l T . The graph shows a decreasing, essentially linear variation of θ t with P e. This indicates that the effect of the diffusion field increases when P e increases, and, again, that it is in opposition to the lateral drifting motion of the pattern. A linear fit of the data indicates that the value of θ t extrapolated to P e = 0 remains slightly smaller than the SP-approximation value (θ sp = 9.07 • ). Nevertheless, the variation of θ t remains very small over the (relatively large) range of scanned P e values.

In conclusion, the anisotropy-driven lamellar tilt angle is very little sensitive to the exact values of λ and P e provided that λ is close to λ m and P e 1. Therefore, simulations performed at a fixed λ value close to λ m , and P e significantly less than unity are fully The angle -δ is also plotted (δ is the angle between σ and z). Anisotropy function

a c (θ) = 1 -0.05 cos[2(θ -θ R )] ; θ R = π/3 (θ sp = 5.3 • ).
d 0 /l T = 2.496 × 10 -3 ; a c = 1 + 0.04 cos[4(θ -θ R )], θ R = 15 • (θ sp = 9.07 • ).
representative of the behavior of the system over a large, experimentally relevant range of the main control parameters. In the following, the BI simulations were performed with P e = 0.024, d 0 /l D = 1.9531 × 10 -5 , l T /l D = 4, and the PF simulations with P e = 0.16, d 0 /l D = 7.9 × 10 -3 , l T /l D = 3.167. In both cases, γ α = γ β = γαβ = 1.

B. Anisotropy functions without missing orientations

Figure 8 shows the variation of θ t as a function of the rotation angle θ R , calculated by both BI and PF, for a simple two-fold anisotropy (m = 2) given by Eq. ( 13). Temperature gradient, growth velocity, and lamellar spacing were held constant. The results of both simulation methods follow quite closely the SP approximation, up to differences that do not exceed a degree. As expected, we find symmetric, non-tilted patterns when either the minimum or the maximum of γ is aligned with z. Somewhere in between these orientations, the value of θ t passes thus through a maximum for an orientation that depends on the anisotropy function. We checked that this maximum increases when the anisotropy coefficient is increased. Overall, however, the tilt angle remains much smaller than θ R over a full rotation range. In other words, a weak crystallographic anisotropy, though it has a clearly detectable effect on the dynamics of lamellar eutectic patterns, is not sufficient to induce a significant locking (θ t ≈ θ R ). For convenience, we will call "unlocked patterns" such weakly anisotropic tilted lamellar patterns. Let us mention again that the values of d 0 /l D , l T /l D and P e are actually different for BI and PF ; the results thus demonstrate once more that the influence of all of these parameters is weak in realistic conditions.

θ sp ##### ##} # # θ s # PF# BI "
FIG. 8: Lamellar tilt angle as a function of the rotation angle θ R . Weak, two-fold symmetry anisotropy of the interphase boundary (a c = 1 -0.05 cos[2(θ -θ R )]). The SP-approximation angle θ sp , and the steady-state angle θ t obtained with the BI and PF simulations are both shown. In this graph, as well as in the following, the represented θ R range is limited to [0, π/2] for obvious symmetry reasons.

Next, we have investigated a situation in which a mild lamellar locking takes place. We use an anisotropy function according to Eq. ( 15), with g = 0.05, w g = 0.195, 2 = 0.0854, 4 = 0.0221 (Fig. 9). This function was chosen so as to smoothly reproduce the Wulff plot extracted (assuming the SP approximation) from the experimental data of Figure 7 of Ref.

[21] (see discussion below). The corresponding Wulff shape is an oval with markedly flattened sides (but without straight facets, and without forbidden orientations). The variation of θ sp with θ R is continuous and univalued. The simulation results for θ t nicely follow the SPapproximation curve. Let us analyze the graph in Fig. 9 in more details. The variation of θ t is essentially linear for θ R ranging from 0 to about 20 degrees (we recall that for θ R = 0, the interphase boundary plane of minimum energy is aligned with the temperature gradient). However, the slope is substantially less than 1 [a linear best fit yields a slope of about 0.76 for θ sp , 0.73 for θ t (BI) and 0.72 for θ t (PF)]. In other words, in this regime, the inclination angle of the interphase boundary remains close to, but departs by a measurable amount from the minimum-energy Both BI and PF simulations reproduce the two separate locked and unlocked branches (Fig. 10a). They demonstrate, in particular, the existence of a strong locking effect over a large orientation range, as predicted by the SP approximation. The two methods exhibit, however, a somewhat different behavior in the bistable interval. In the BI simulations, starting from the situation where the facet is aligned with the temperature gradient (θ R = 0), the system closely follows the locked branch upon increasing θ R , up to a limit angle at which it jumps abruptly to the unlocked branch. The jump occurs well before the turning point θ l , at a value of θ R which is close to 45 • . This point can be only attained if θ R is successively increased along the locked branch, and the maximum value of θ t depends on the step size in θ R . It is thus clearly initial-condition dependent and has no obvious connection with the sharp-edge angle θ e . This limit value of θ t in the simulations is also compatible with previous in situ observations of tilted lamellar patterns. Reversely, when θ R is decreased stepwise starting from θ R = π/2, the system describes the whole unlocked branch, within numerical accuracy. This hysteresis is expected for a bistable system. It should be noted that the approach of the limit tilt angle along the locked branch also corresponds to a steep, apparently diverging increase of the δ angle (Fig. 10b). This indicates that, at large tilt angles, the SP approximation, while it still correctly predicts the values of the tilt angle, becomes inaccurate with regards to other aspects of the dynamics. For instance, in contradiction with its basic assumption, the shape of the solid-liquid interface becomes markedly asymmetric at large tilt angles, as can be seen in Fig. 10c.

In contrast with the BI simulations, no hysteresis is observed in the PF simulations. More precisely, no states on the locked branch have been observed for rotation angles larger than about 35 • : when θ R is increased beyond that value, the lamellae switch to the unlocked state. The opposite jump from the unlocked to the locked state was found to occur at the same value of θ R (with an uncertainty of about 1 • ) in runs that were started on the unlocked branch and in which the rotation angle was successively decreased. A refinement of the grid and a decrease in the interface thickness by a factor of two produced a change in this critical value that did not exceed two to three degrees. In contrast, an increase of the amplitude g of the Gaussian in the anisotropy function led to a marked increase of the discontinuous-jump angle.

V. DISCUSSION

A. Reconstructing the anisotropy function

A practical aim of the present study was to give a numerical support to a recent experimental work based on a rotating directional solidification (RDS) method [20,21]. The RDS method uses a standard thin-sample directional solidification setup, and permits, in addition, to rotate the sample at constant angular speed about an axis perpendicular to the (two-dimensional) sample. Under a few, not very restrictive, conditions (zero translation speed, center of rotation placed on the eutectic isotherm, quasi-two-dimensional and quasistatic nature of the front pattern dynamics), this is equivalent to continuously varying the rotation angle θ R of a given eutectic grain with respect to the thermal gradient axis, as we have done in the above calculations. More precisely, the lamellar tilt angle θ rds observed over time during a correctly set-up RDS experiment is equal to the steady-state tilt angle θ t at the current value of θ R . Moreover, in the SP approximation, this entails that the RDS trajectories of the trijunction points are centrosymmetric closed curves homothetic to the section of the Wulff shape of the interphase boundaries by the sample plane, from which a two-dimensional anisotropy function of the interphase interface can be derived (see Ref. [20]).

The anisotropy function used in Fig. 9 was derived from the RDS pattern of a nearly locked grain by this method. As a test for the accuracy of this method, we have reconstructed the Wulff shape from the calculated θ t (θ R ) points under the SP assumption. The two shapes are identical to within experimental error (a few percents of the position vector), as can be seen in Inset 2 of Fig. 9. In conclusion, the errors due to the SP approximation are generally not larger than the experimental uncertainties, which validates the use of the SP approximation in the exploitation of the RDS patterns.

B. Bistability in the numerical simulations

The most important difference between the results of BI and PF simulations is the absence of bistability in the θ t (θ R ) curves obtained with the PF code. As discussed previously, such a bistable behavior should follow from the existence of states of positive stiffness on the "ears" of the Wulff plot. In addition, there is clear evidence of bistabilty in experimentally observed RDS patterns of strongly locked grains, as illustrated in Figure 11 (also see Ref.

[21]). Therefore, the BI method seems to be in better agreement both with the SP theory and with experiments.

Whereas quantitative differences between phase-field and sharp-interface models have often been reported, such a strong qualitative difference between the two methods is a striking finding. Let us first discuss -and, actually, rule out -one possible source for this difference. It is well known that PF models need to be regularized for anisotropy functions that generate missing orientations. Indeed, for orientations with negative stiffness, the evolution equations for the phase fields become ill-posed. This corresponds, in the free-boundary problem, to unstable growth modes with arbitrarily high growth rates in the limit of vanishing wavelengths (absence of stabilization by capillarity). This behavior can be regularized by a convexification of the polar plot of 1/γ [29], or by the addition of higher-order derivatives in the free-energy functional [START_REF] Wheeler | Proc. R. Soc. A[END_REF][START_REF] Torabi | [END_REF]. Since the former method requires a re-calculation of the convexification for each new choice of anisotropy function, we have preferred to use the latter by following Torabi et al. [START_REF] Torabi | [END_REF], who have added to the free energy the square of the mean interface curvature. For our model we used the linearized form, which amounts to adding to the functional of Eq. ( 18) a new term of the form,

γαβ β 2 V ( ∇ • ( q αβ )) 2 , ( 28
)
where β is a length scale over which the sharp corners of the equilibrium shape are smoothed out.

As expected, with the help of this regularization, we have obtained the correct equilibrium shape for a β inclusion inside an α matrix, even for anisotropy functions with missing orientations such as the example used for generating Fig. 10. However, this modification did not appreciably alter the results for lamellar growth: still no bistable behavior was observed in the PF model, and the value of the angle at which the "jump" from the locked to the unlocked branch occurs was not appreciably modified.

We believe that the origin of the difference between BI and PF results is the behavior of the diffuse trijunctions in the PF model. Indeed, in Ref. [10], it was found by direct comparison between PF and BI simulations that the dynamics of diffuse trijunctions deviates from the predictions of sharp-interface models: since the solute diffusivity remains non-zero within the predictions and experimental observation, but is not found in the PF model.

The numerical simulations demonstrate that the real growth shapes are not exactly symmetric, but that this departure from the hypothesis that underlies the SP theory leads only to small differences in the steady-state tilt angle, for various choices of anisotropy functions.

Moreover, these differences are virtually independent of the growth velocity and local spacing, at least in the regime of slow growth (small Peclet numbers) and for spacings close to the minimum-undercooling spacing. Since these conditions are found in many experiments, these results firmly establish the SP theory as a useful tool for predicting the growth angles of lamellar eutectics. In particular, experiments with the rotating directional solidification setup can be used to obtain quantitative information on anisotropy functions, as discussed in Sec. V. It should be noted that in the present numerical study we have limited ourselves to a symmetric phase diagram and to growth at the eutectic composition; however, since the tilt angle is mainly fixed by the (anisotropic) capillary effects, we do not expect major differences for other phase diagrams or compositions.

The important qualitative difference between the predictions of BI and PF simulations is an interesting finding in itself, since correctly designed PF models are usually a faithful representation of the corresponding free-boundary problems. While we have not succeeded in pinpointing the exact origin of the difference between PF and sharp-interface models, we believe that the key point is the dynamic behavior of diffuse anisotropic trijunctions. This is an interesting subject for further studies.

Despite this open problem, in the regime of low anisotropy (without missing orientations), BI and PF simulations show an excellent agreement. This implies that these models can be used to explore the behavior of anisotropic eutectics beyond steady-state growth. In particular, it would be interesting to address the effect of anisotropy on the various instabilities that are known to occur in lamellar eutectics: short-wavelength oscillations [3,6],

spacing diffusion [8,9], and the zig-zag instability [4,11]. Without any doubt, the stability boundaries will depend on the strength of the anisotropy, but it seems difficult to predict in which way without the help of numerical simulations. Finally, it is straightforward to implement anisotropic PF models in three dimensions, which opens the possibility to explore numerically the anisotropy effects on eutectic growth in bulk samples.

Equations (A3) and (A4) can be inverted to express the constants K and H in terms of the surface free energy γ and the length scale = 4W/π:

K = γ , (A5)
H = 16γ π 2 . (A6)
It can easily be verified that this calculation remains valid for a binary interface (an interface between phases i and j along which φ i + φ j = 1 and all other phase fields are zero) in the multi-phase-field setting; therefore, the parameters γ ij in Eqs. ( 20) and ( 21) are directly the surface free energies of the respective interfaces.

In order to generate a surface free energy that depends on the interface orientation according to Eq. ( 9), Eq. (A4) can still be used, but the coefficients K and/or H need to be orientation-dependent. The "standard" procedure (see for example [32]) is to keep H constant and to write K(θ) = Ka 2 c (θ). According to Eq. A3, this creates variations in the interface thickness W , which becomes proportional to a c . This can be avoided by letting K(θ) = Ka c (θ) and H(θ) = Ha c (θ). Finally, it is also possible to keep the gradient energy coefficient constant and to write H(θ) = Ha 2 c (θ), which leads to W ∼ 1/a c (θ).

It should be noted that each of these choices generates different equations of motion for the phase field. The functional derivative, for any functional of the form F = V f (φ, ∇φ, ∇ 2 φ) is explicitly given by

δF δφ = ∂f ∂φ - i=x,y,z ∂ i ∂f ∂(∂ i φ) + i,j=x,y,z ∂ 2 ij ∂f ∂(∂ 2 ij φ) , ( A7
)
where ∂ i φ denote the Cartesian components of ∇φ. Since the interface orientation θ can be expressed as a function of ∇φ, the second term on the right-hand-side acts on any θdependent term in the functional. We have implemented the evolution equations corresponding to the three possibilities outlined above, and have compared the results for a few selected examples. We have found no significant differences. The results presented in the main text are obtained with the anisotropy function in the square gradient term only. In order to treat crystals that are rotated with respect to the temperature gradient (laboratory frame), we use coordinate transformations involving rotation matrices, as discussed in appendix B below.

The equations of motion for all the fields in our multi-phase-field model are obtained from the functional Ω by the procedures detailed in Refs. [22,23] Finally, J at is the antitrapping current that is added to counteract spurious solute trapping effects and guarantee the correct thin-interface limit [33]. Here, we combine the appropriate expression for the double-obstacle potential of Ref. [23] with the interpolation between multiple growing phases developed in Ref. [10].

The time evolution of the phase-fields follows a relaxation equation,

∂φ ν ∂t = - 1 τ δΩ δφ ν -Λ , (A10)
where δΩ/δφ denotes a functional derivative, and Λ is a Lagrange multiplier that is added to ensure N ν=1 φ ν = 1 throughout the system. The relaxation coefficient τ may vary between function is given by a c ( q ) = 1 -4 3 -4

q x 4 + q y 4 |q | 4 , ( B3
)
where 4 is the anisotropy strength; a two-fold anisotropy is generated by

a c ( q ) = 1 -2 q x 2 -q y 2 |q | 2 (B4)
for an arbitrary function a c (θ), we may use that tan θ = q y /q x .

In order to calculate the functional derivative in Eq. (A10), repeated use of the chain rule is made. For instance, we have (B6)

FIG. 1 :

 1 FIG. 1: Lamellar eutectic patterns observed in situ during thin-sample directional solidification (V = 0.5µms -1 ) of a eutectic transparent (CBr 4 -C 2 Cl 6 ) model alloy. a) Symmetric steady-state in a "floating" eutectic grain. b) Tilted lamellae in a "locked" eutectic grain. The growth direction is vertical (liquid on top). Bar: 20 µm.

FIG. 2 :

 2 FIG.2: Schematic repeat units of lamellar eutectic patterns. a) Isotropic system. b) System with an interfacial anisotropy of the interphase boundary in the symmetric-pattern (SP) approximation. α, β: solid phases. L: liquid. z: growth direction parallel to the thermal gradient. x: direction of the isotherms. θ sp : SP-approximation lamellar tilt angle. The lateral drift velocity is given by tan θ sp = V d /V , with V the pulling velocity. Other symbols: see text.

FIG. 4 :FIG. 5 :

 45 FIG. 4: Steady-state lamellar patterns (BI simulations). a) Symmetric pattern without anisotropy. b) Tilted pattern with anisotropic interphase boundaries (anisotropy function given by Eq. (14); g = 0.2 ; w g = 0.1; θ R = 3π/48 ≈ 11.25 • ; P e = 0.024 ; θ t = 11.0 • ). c) Spatio-temporal diagram showing the steady-state dynamics corresponding to a) and b), successively, and the brief transient after the anisotropy was turned on.

FIG. 6 :

 6 FIG.6: Steady-state tilt angle θ t obtained from BI simulations as a function of the reduced lamellar spacing λ/λ m , at constant G and V (d 0 /l D = 1.9531 × 10 -5 , l T /l D = 4).The angle -δ is also plotted (δ is the angle between σ and z). Anisotropy functiona c (θ) = 1 -0.05 cos[2(θ -θ R )] ; θ R = π/3 (θ sp = 5.3 • ).

FIG. 7 :

 7 FIG. 7: Lamellar tilt angle vs Péclet number at constant λ/d 0 = 20.25, d 0 /l T = 2.496 × 10 -3 ; a c = 1 + 0.04 cos[4(θ -θ R )], θ R = 15 • (θ sp = 9.07 • ).

FIG. 9

 9 FIG. 9: a) Lamellar tilt angle as a function of the rotation angle θ R . Mild lamellar locking effect caused by a shallow, smooth local minimum in the Wulff plot of the interphase boundary (see text). Same symbols as in Fig. 8. Dashed line: slope 1. b) Wulff shape (θ R = 0). c) Partial view of the Wulff shape (dash-dotted line) and the shape reconstructed from the θ t data of Fig. 9 under the SP approximation (thin line).

FIG. 10

 10 FIG. 10: a) Lamellar tilt angle as a function of the rotation angle θ R . Strong lamellar locking effect with a sharp local minimum in the Wulff plot (see text). Same symbols as in Fig. 8. Inset: Wulff plot (thin line) and Wulff shape (thick line). Dashed lines: unstable branches. b) Angle δ of σ with z as a function of θ R (BI simulations). c) Shape of the lamellar pattern with the largest θ t value (≈ 48 • ) simulated with the BI code.

FIG. 11 :

 11 FIG. 11: Experimental observation of a coexistence of eutectic-growth domains with two different tilt angles in a single eutectic grain. Rotating directional solidification of a thin (10-µm thick) sample of a eutectic CBr 4 -C 2 Cl 6 alloy. Horizontal dimension: 450 µm.

  . The chemical potential obeys a diffusion equation with source term,χ(φ, µ) ∂µ ∂t = ∇ • M (φ) ∇µ -J at + φ), (A8)where χ(φ, µ) = -∂ 2 ω b (φ, µ)/∂µ 2 is a generalized susceptibility; for our choice of ω b , χ(φ, µ) ≡ 1/2A. Additionally, this choice leads to equal slopes of the liquidus and solidus lines, which renders∂ 2 ω ν ∂µ∂φ T = -u ν ∂h ν (φ) ∂φ and ∂ 2 ω ν ∂µ∂T φ = 0.For an imposed temperature gradient G translated at a given velocity V , ∂T ∂t = -GV from Eq. (1). Furthermore, M (φ) is the atomic mobility of the solute, which is assumed to be zero in the solid and finite in the liquid. Using the standard relationship between mobility and diffusion coefficient,D = M ∂ 2 f /∂u2 , we can write M (φ) = φ l D 2A . (A9)

ACKNOWLEDGMENTS

We thank Alain Karma for sharing with us his original BI code. This work was financially supported by Centre National d' Études Spatiales, France.

plane, which defines a mild locking effect. For rotation angles larger than, say, 45 degrees in Fig. 9, the lamellar growth dynamics escapes the influence of the Gaussian minimum of the Wulff plot. In this unlocked-pattern region, the value of θ t is solely determined by the cos 2θ and cos 4θ terms in the a c function (in this θ R range, we checked that θ t remains equal to 0, as it should, when 2 and 4 are set to 0). In an intermediate θ R interval, there is a steep, but smooth transition between a nearly-locked and an unlocked behavior. In this crossover range, the departure of θ t from the SP approximation is larger than in the rest of the graph. The value of the angle δ as it is given by BI simulations more or less follows the same variation as θ t as a function of θ R (not shown). It reaches a maximum (of about 2 • ) in the crossover region, but this maximum remains small as compared to the corresponding θ t value (≈18 • ).

C. Anisotropy with missing orientations

A strong lamellar-locking effect (Fig. 10a) can be reproduced by using an anisotropy function of the same form as that of Fig. 9, but with a deeper and sharper Gaussian ( g = 0.2 and w g = 0.1 ; see Eq. 14). This modification of the Wulff plot entails the appearance of two (quasi) facets, and four "ears" with long metastable branches and sharp-edge junctions in the Wulff shape of the interphase boundary (see inset in Fig. 10a). Let θ u and θ l , where θ u < θ l , be the tilt angle values at which γ + γ = 0 on the interval [0, π/2] (θ u ≈25 • and θ l ≈ 70 • in the example shown in Fig. 10). The SP approximation predicts three distinct parts in the θ sp vs θ R curve: (i) an essentially linear strongly locked branch with a slope close to 1, which runs from θ R = 0 to θ R = θ l ; (ii) an unlocked, although (weakly) anisotropic, branch for θ R ranging from θ u to π/2; (iii) an intermediate branch, which connects the end points of the locked and unlocked branches, and is, presumably, not observable given that the interface boundary is unstable (γ + γ < 0) along its entire length. In the [θ u , θ l ] interval, two (locked and unlocked) values of the lamellar tilt angle are possible for a given eutecticgrain orientation. In this bistable range, there is a value θ e of the rotation angle (here, approximately 45.3 • ) at which σ has the same value for both branches -this corresponds to a sharp edge in the convex equilibrium shape. diffuse interfaces, the solid-solid interface can actually move close to the trijunction point, contrary to the assumptions made in the one-sided model of solidification, in which the diffusivity becomes zero immediately behind the (sharp) interface. Moreover, a rotation of the trijunction point by a finite angle was observed that persisted even in the sharp-interface limit [10]. Since such effects are present even for isotropic interfaces, one may expect them to be even more important here, where the anisotropic surface energy creates strong "Herring torques". Indeed, we have observed in our simulations that the solid-solid interface tends to slightly change its orientation upon approaching the trijunction region. Therefore, it is possible that the PF model explores a wider range of orientations within the trijunction region than the BI model, which could facilitate the switching between different solution branches.

It should also be noted that the regularization outlined above was developed and validated only for simple interfaces between two phases, and may therefore not be complete in the vicinity of trijunction points. It is possible that the bistable behavior observed in the BI simulations could be recovered by a PF model with correctly regularized trijunction points.

However, since the analytic understanding of trijunction points currently is very limited, it is not clear how such a regularization should be carried out. For all these reasons, a more detailed study of moving diffuse trijunctions certainly is an interesting subject for future work.

VI. CONCLUSIONS AND PERSPECTIVES

We have investigated lamellar eutectic growth with anisotropic solid-solid interphase boundaries, using two different numerical methods: a dynamic boundary-integral and a phase-field code. We have obtained good agreement between our numerical results and the prediction of a recent approximate theory for the growth direction as a function of anisotropy and orientation of the growing bicrystal, which uses the symmetric pattern approximation.

For smooth anisotropy functions (no missing orientations), the two methods are in excellent agreement. If the anisotropy is strong enough to induce missing orientations around a deep minimum in the γ plot, BI exhibits a bistable regime that is in agreement with theoretical Here, we will discuss in more detail several technical points about our phase-field model; these are mostly known facts, but it seems useful to recall them briefly here. We will start with the implementation of the anisotropy, and then give the equations of motion for our model.

Consider a two-phase system characterized by a single phase field φ, with an interfacial free-energy functional of double-obstacle type

where K and H are constants and φ is restricted to the interval [0, 1]. The equilibrium interface profile is given by

for an interface normal to the x direction centered at the origin, with

Furthermore, the standard evaluation of the surface excess free energy (the additional energy created by the presence of an interface) yields

the values τ ij for each interface are chosen such as to make the interface kinetics vanish [23];

the interfaces thus obey the local equilibrium condition of Eq. ( 5).

Appendix B: Rotation matrix representation of Anisotropy

In two dimensions, the interface orientation θ can be obtained for an i-j interface from the Cartesian components of the vector q ij = φ i ∇φ j -φ j ∇φ i . Some simple anisotropy functions a c (θ) can also be directly expressed in terms of the components of q ij instead of θ.

If the growing crystal is rotated, in two dimensions the anisotropy function is changed from a c (θ) to a c (θ -θ R ). We wish to use an implementation that can easily be extended to three dimensions and arbitrary rotations. In order to achieve this, we directly work with the components of the interface normal vector and transform the equations to the reference system of the crystal with the help of the appropriate rotation matrix.

Consider an α-β interface and let the vector q = φ α ∇φ β -φ β ∇φ α (we drop the phase indices for ease of notation), with q x and q y being its Cartesian components in the laboratory frame. We rotate this vector by an angle θ R , which brings it into the reference frame of the crystal,

where R is the rotation matrix

In the reference frame of the crystal, the anisotropy can be expressed in the "primed" vector components, independently of the rotation angle. For example, a simple cubic anisotropy