Jihène Ghmadh

Jean-Marc Debierre

Marc Georgelin

Rahma Guérin

Alain Pocheau

Physical mechanisms of the oscillatory 2λ-O mode in directional solidification

Keywords: numbers: 81.10.Aj, 68.70.+w, 81.30.Fb

The 2λ-O oscillatory mode of cellular solidification patterns is studied in thin samples of a succinonitrile-acrylonitrile transparent alloy. The origin and the nature of oscillations are revisited and discussed by combining experiment with 3D phase-field numerical simulations. The existence domain of 2λ-O oscillations and the evolution of their period with growth velocity are determined and compared. Simulations evidence transversal solute fluxes between neighbor cells as an essential feature of cell dynamics. A solute balance model in which transversal fluxes are crucial for oscillations recovers the emergence of a 2λ-O mode and its period-velocity relationship. It thus confirms the fundamental role of transversal fluxes and provides a complete coherent description of the physical mechanisms of 2λ-O oscillations. Parametric excitations are finally used to force 2λ-O oscillations beyond their stability domain and highlight the nature of the underlying oscillator, especially its non-linearity responsible for intermittent oscillations and complex behavior in the resonance band.

I. INTRODUCTION

In physical, chemical or even biological systems, a canonical mechanism for the emergence of oscillatory patterns stands in the destabilization of a fixed point into a limit cycle [1,[START_REF] Misbah | Dynamique complexe et morphogénèse[END_REF]. In directional solidification, spectacular examples of such oscillatory growth modes are provided in either the isotropic phase of a cholesteric liquid crystal [START_REF] Cladis | [END_REF], eutectic alloys [4][5][6], or diluted binary alloys [7][8][9][10].

In the latter case, 2λ-O oscillations have been observed in thin samples of a diluted succinonitrile-based alloy [7]. Adjacent cells were found to form pairs (spacing 2λ) that oscillate (O) in phase opposition while remaining left-right symmetric. The stability of cellular arrays with respect to various modes including 2λ-O oscillations has been determined by a Floquet-Bloch analysis and found to critically depend on surface tension [11]. Related vacillating-breathing modes have also been observed in numerical studies based on amplitude equations [12,13] or on a phase-field model [14].

In the above experiment [7], a large scan of velocity and cell spacing provided 2λ-O oscillations on a bounded velocity range 1.0 < V /V c < 4.5 close to the critical velocity V c . As several instabilities compete in this velocity domain [11,[START_REF] Georgelin | Étude expérimentale des instabilités secondaires en croissance cristalline directionnelle[END_REF], saturated oscillations lasting ten or more periods were seldom observed. In particular, the large cell deformations displayed during an oscillation period offered the solidification front the opportunity to evolve toward asymmetric structures like doublons [START_REF] Jamgotchian | [END_REF]. Meanwhile, both the spatial and the temporal symmetries broke since the system turned from an oscillatory to a stationary growth mode [7]. This corresponds to a rare experimental example of such a doubly-broken symmetry. Following this destabilization to doublons, the 2λ-O mode is usually impermanent in practice, so that observing it requires some patience, both in experiments and in simulations. In this context, combining both experiment and numerical simulation could therefore pro-vide valuable insights into its mechanism.

More recently, experiments were performed under microgravity conditions in extended samples of a different succinonitrile-based alloy. They provided globally disordered arrays of oscillating cells displaying locally correlated oscillations with a π phase-shift in the case of local square ordering or a 2π/3 phase-shift in the case of local hexagonal ordering [8][9][10]. Here too, these spatiotemporal patterns refer to generic modes [START_REF] Misbah | Dynamique complexe et morphogénèse[END_REF]. On the other hand, numerical simulations were performed in extended systems, first by using amplitude equations [17], later by using the two-sided phase-field model for a model system [18]. Parameterizing the one-sided phase-field model with the actual physical parameters of microgravity experiments, π and 2π/3 oscillatory modes were recently recovered with oscillation periods comparable to the experimental ones [8][9][10]. This strongly suggests that quantitative numerical results could also be obtained for the 2λ-O oscillating mode.

For thin samples, a pioneering study using the onesided phase field model [19] reported the onset of an oscillatory single structure in an intermediate domain of spacings but, to our knowledge, no simulations were performed so far to specifically study the 2λ-O mode. We report here on new experimental results and numerical simulations relying on a phase-field code designed to incorporate the most salient experimental characteristics. A very similar code was recently validated [20] by means of quantitative comparisons with the orientational response of inclined solidification structures growing in comparable experimental systems [21][22][23]. The good agreement found there motivates us to adopt here a similar strategy for studying the 2λ-O mode by combining experiments with numerical simulations.

Altogether, the experimental and numerical studies performed here provide new insights into both the existence domain of 2λ-O oscillations and the relation between the oscillation period T and the pushing velocity V . In addition, a parametric excitation of the 2λ-O mode, performed both in experiment and simulation, reveals an intermittent occurrence of 2λ-O oscillations in the resonance band which emphasizes the intrinsic nonlinearity of this mode. On the other hand, an analysis of the dynamics of the concentration field revealed by simulation points out the importance of transversal fluxes between cells. On this basis, a simple model of cell dynamics based on concentration balance leads the 2λ-O oscillations to be recovered. Its underlying physical mechanisms emphasize the role of cell coupling by transversal fluxes in the emergence of oscillations of a cell pair.

This paper is organized as follows. Section II describes the experimental procedure and the 2λ-O mode, while section III reports the basic ingredients of the phase-field code used in the simulations. Section IV is devoted to the stability diagram and to the period-velocity relation. The data extracted from the concentration field are analyzed in section V and the essential role of transversal fluxes between cells is identified. Their implication is analyzed in section VI within a model of concentration balance and of cell evolution which succeeds in recovering 2λ-O oscillations. Section VII reports the parametric excitation of the 2λ-O mode worked out in experiment and simulation and emphasizes the role of non-linearity in the intermittent occurrence of 2λ-O oscillations. Finally, a discussion completed by our main conclusions is provided in section VIII.

II. EXPERIMENTAL

The experimental setup is designed to achieve directional solidification of a layer of cells or dendrites in homogeneous and controlled conditions capable of allowing the real-time visualization of microstructures dynamics.

Details of the classical directional solidification setup are described elsewhere [22,24]. A thin sample filled with the solidifying material is pushed in a controlled uniform thermal gradient towards the cold zone at a given velocity V . The growing liquid-solid interface is then observed through an appropriate optical stage. The sample is made of two parallel glass plates glued on their sides to define a 45mm×100mm×100µm thin cavity. Sample thickness is selected by using mylar sheets as calibrated spacers, thin enough to avoid the emergence of a second layer and thick enough to ensure a 3D behavior of microstructures as opposed to the 2D ribbon-like behavior displayed for excessive squeezing [19]. A thermal gradient of 110 K.cm -1 is provided by two heaters and two coolers which sandwich the sample, all being electronically regulated with an accuracy better than 0.1 • C. When addressing the periods of the 2λ-O mode, it has also been extended to 140 K.cm -1 . The material used in experiments and modeled in numerical simulations is a transparent plastic crystal, succinonitrile, with a small amount of acrylonitrile as dilute solute. Care has been taken to obtain a single crystal without grain boundaries in the whole sample, the selected orientation involving a [0 0 1] direction perpendicular to the sample plane and a [1 0 0] direction parallel to the thermal gradient (and thus normal to isothermal lines and planar fronts). Physical and control parameters of experiments and simulations are given in Tab. I. On a given sample, the solidification velocity is increased up to a constant value V at which observations are performed. The velocity growth rate influences the mean cell spacing, narrow spacings being displayed at large rates and wide spacing at small ones. In any case, sufficient time is given to spacing inhomogeneities to relax by diffusion, thereby providing a quasi homogeneous spacing over the observation window which typically spans over a dozen of cells. Instabilities, as the 2λ-O instability studied here, occur naturally by development of fluctuations usually on steady patterns but sometimes also on dynamical ones, following the perturbations induced by a previous instability.

Figure 1 displays a typical sequence of oscillation evidenced on a cycle of the 2λ-O mode at V = 8µm.s -1 . (See the corresponding movie in [25] as well as the movie for G = 140K.cm -1 with a larger number of cells.) The spatial period of oscillation extends over two cell spacings with neighboring cells in phase opposition. This actually corresponds to a 2λ instability. Oscillations involve both the cell shape and the cell position (or undercooling). As displayed in figure 1-a, when cell positions are the same, the shape differences are the largest. This means that their oscillations are nearly in phase quadrature. Typical values of the oscillation periods are 1 or 2 min. This stands as an unusually long characteristic time for cellular instabilities since, for instance dendritic oscillations only last a few seconds [24].

Figure 2 displays a cycle provided by oscillating cells over many periods. To enhance the signal, the difference between the tip positions δz t (t) and the cell widths δλ(t) of two adjacent cells is considered. On these oscillating cells, the cell widths are measured at the position of their tip curvature center.

At a given V , a scan of the 40 mm wide front displays both steady and oscillating cells. As the cell spacing is quite uniform even on such a large distance, it appears that at the same values of V and λ, steady and oscillating cells can coexist [7]. This makes a difference with other cell instabilities such as cell elimination, tip splitting or dendritic emission which develop as soon as the required conditions in (V , λ) are satisfied. Here, the coexistence of steady and oscillating cells suggests a sub-critical nature of the 2λ-O instability. This is supported by the response to a sudden velocity pulse which makes previously steady cells develop 2λ-O oscillations.

In practice, the fact that some cells remain steady whereas other similar cells have undergone the instability means that the former have not yet encountered the perturbations that made the latter go unstable. The reason is that the perturbations yielding the instability are rare, for instance because they are intense or spatially coherent. Intense fluctuations are actually required in sub-critical instabilities to leave the basin of attraction of a stable fixed point. They then involve a threshold in intensity that makes their probability of occurrence low. Here, in addition, one may suspect that some level of spatial coherence over a couple of cells is required to make oscillations develop. In particular, no oscillation of a single cell has ever been observed and the model of section VI will find them damped. Accordingly, one may guess that fluctuations involving a characteristic length of 2λ are required, with a positive perturbation on a cell and a negative one on its neighbors. Here again, this specificity reduces the probability of occurrence.

In this context, two kinds of procedure have been used to document the 2λ-O mode. The first one consists in applying a slow ramp or small increments of velocities and in observing a given domain of a dozen of cells over a long time to detect and follow the birth and the development of the oscillations. As this is time and sample consuming, this cannot be used systematically. The second procedure consists in scanning the whole front at a given velocity V so as to detect yet developed oscillations. It has been used to fill in the existence diagram in the (V , λ) space (Sect.IV), the spacing λ being determined from the distance between the cell axes. Both procedures have their own merit : the first one enables to detect 2λ-O oscillations from their birth to their longterm development ; the second one enables to determine the existence diagram in a reasonable time.

For velocities close to the upper limit of the existence domain, oscillations succeed in keeping the same amplitude over about ten periods (Fig. 2). However, for velocities below, i.e. in most of the existence domain, 2λ oscillations display an increasing amplitude which yields one of the oscillating cells to either tip-split or get eliminated depending on whether its initial spacing is large or small. In both cases, neighboring cells relax toward an asymmetric state corresponding to a so-called doublon that is stable with respect to the 2λ-O mode [7]. In this respect, the 2λ-O instability stands as a mediator for the formation of doublons. This emphasizes its long-term role on the cellular pattern dynamics.

III. PHASE-FIELD

The thin interface phase-field model (TIPM) was introduced by Alain Karma and co-workers to simulate the solidification of a pure substance [26,27]. This model was extended later to the case of a dilute binary alloy [28,29]. Recently, direct comparisons of the TIPM with 3D experiments of alloy solidification performed in confined samples [19,20] and in extended ones [8][9][10] became feasible, and quantitative agreement was found in each case.

Here, as in a number of recent studies [8-10, 20, 30, 31], we replace the usual phase-field ϕ ∈ [-1, 1] by the preconditioned phase-field which is a signed distance from the interface. For a given numerical accuracy, this transformation allows one to use coarser numerical grids, which significantly reduces the numerical effort [32]. Graphics Processing Unit (GPU) parallel programming is used to further increase the code efficiency [8][9][10]33]. We only give the main lines of the model here. For more details, we refer the reader to our recent study of inclined solidification structures under very similar experimental conditions [20]. Numerical simulations are performed in a thin parallelepipedic domain, x lying along the sample width, y along the sample thickness, and z along the temperature gradient.

ψ = √ 2 tanh -1 (ϕ), (1)

A. Evolution equations

In the TIPM, the solid-liquid interface is given a thin but finite width W 0 which is used as the unit length here [26,27]. For rough materials like succinonitrile, attachment kinetics can be neglected. We thus impose a zero kinetic coefficient here by setting the characteristic time of the phase-field model to τ 0 = a 0 W 3 0 /(Dd 0 ), where d 0 is the capillary length, D the solute diffusion constant in the liquid phase, and a 0 0.5539 [26,27]. In order to obtain a nondimensional version of the evolution equation for ψ, we follow the usual convention where lengths and times are respectively divided by W 0 and τ 0 . The resulting equation is

(1 -β k z * )a 2 s ∂ψ ∂t = √ 2 ϕ -C(1 -ϕ 2 )(U + z * ) +a 2 s ∇ 2 ψ - √ 2ϕ(∇ψ) 2 + 2a s ∇a s • ∇ψ + √ 2 (1 -ϕ 2 ) ∇ • A, (2)
with β k = (1 -k), k being the partition coefficient. Here, the nondimensional concentration field U is related to the physical concentration field c by

c = c ∞ k 1 + β k U α k -β k ϕ 2 , (3)
where α k = (1 + k), c ∞ is the nominal solute concentration, and the additional term β k z * on the left hand side of Eq. ( 2) is introduced to reduce higher order corrections that appear in the asymptotic expansion of the phase-field equations [29]. The nondimensional variable z * is defined as

z * = [z + z s (t) -V t]/l T , (4)
where

l T = mc ∞ (k -1)/(kG) (5)
is the thermal length. The quantity V t-z s is a small vertical offset that becomes constant (up to one mesh size) at large times. It is the difference between the physical pushing velocity term V t and the numerical term z s (t) which adds up the vertical shifts performed during the simulation (as explained in subsection B).

For our rough material, the constant that couples the nondimensional concentration field U to the phase-field is C = 75D * /47, with D * = Dτ 0 /W 2 0 [27]. Both a s and A depend on the crystal anisotropy. In the present case, the anisotropy is cubic and the crystal axes [1 0 0] and [0 0 1] are oriented as in the experiment (see section II). The corresponding anisotropy function reads then [34]

a s = (1 -3 4 ) + 4 4 (n 4 x + n 4 y + n 4 z ), (6)
where n x , n y , n z are the components of the unit vector n along the normal to the solid-liquid interface and 4 is the anisotropy strength. Moreover, the three components of the anisotropy vector A are given by

A µ = 16 4 (1 -ϕ 2 ) √ 2 |∇ψ|a s n µ (n 4 x + n 4 y + n 4 z ) -n 2 µ , (7)
with µ = x, y, z.

We use here the one-sided model that neglects solute diffusion in the solid. Following references [28,29], a corrective solute current

j at = W 0 2 √ 2 c ∞ k β k γn, (8)
is used to avoid artificial solute trapping and other spurious corrections due to the finite interface thickness W 0 . Here

γ = 1 -ϕ 2 √ 2 (1 + β k U ) ∂ψ ∂t . (9)
The evolution equation for the nondimensional concentration field U is then [28,29]

α k -β k ϕ ∂U ∂t = (1 -ϕ)D * ∇ 2 U - (1 -ϕ 2 ) √ 2 D * ∇ψ • ∇U - 1 √ 2 n • ∇γ + γ 1 - ∇ • n √ 2 . (10)

B. Parameters and conditions imposed

Along the y direction (sample thickness), no-flux boundary conditions are imposed on U while a contact angle θ c is prescribed by imposing ∂ψ/∂y at both boundaries (see for instance [35] for a detailed description of such limit conditions). We use here the same value as in [20], ∂ψ/∂y = ±0.8333, the positive (negative) value corresponding to the front (rear) boundary, so that θ c 33.5 • . We verified that for such small contact angles the actual value of θ c has a weak influence on the simulation results [START_REF] Ghmadh | Étude par la méthode du champ de phase à trois dimensions de la solidification dirigée dans des lames minces[END_REF]. Periodic boundary conditions are preferred for the x direction (sample width). Along the z direction (thermal gradient), the boundary conditions are no-flux at the bottom of the domain and U = -1.0 at its top. In addition, the fields ψ and U are shifted when needed to keep the interface roughly at the same vertical location in the simulation domain. In the stationary regime, the cumulated shift z s (t) = V t + z 0 , where z 0 is a constant offset due to the fact that the growth velocity varies during the initial transient. The initial conditions used here are either a randomly perturbed flat interface or a 2λ-O mode obtained from a previous simulation on which a dilation or a contraction of space is applied in order to vary the spacing λ. Both lead to the same features for the developed steady or oscillating state. The convergence parameter ξ = W 0 /d 0 has been tested on steady cells until reaching constant cell features. It has also been taken sufficiently small that no spurious side-branching occurs on the cells and that the oscillation period T remains constant within a few percent if ξ is further decreased by about ten percents. Altogether, a satisfactory convergence is obtained for ξ = 90.0.

A spatiotemporal noise of small amplitude F [8-10, 20, 31, 37] is introduced during the early stages of the simulations (t ≤ 100s). Afterwards, noise is turned off (F = 0) so that only the intrinsic numerical noise remains : all the results presented here exclusively correspond to this F = 0 regime. Although the growth time of the 2λ-O oscillations may depend on the details of the initial noise, the features of the developed oscillating states do not. In contrast, the simulations of parametric oscillations considered in section VII are always performed without spatiotemporal noise from the beginning (F = 0) to avoid introducing a bias in the forcing.

The numerical parameters used in the simulations are gathered in table II. Based on our previous study [20], we use 4 = 0.011 [38] for the crystal anisotropy and 35µm for the sample thickness. For this thickness, the cell tip has a curvature radius of the order of ten micrometers in the yz plane perpendicular to the sample (Fig. 3). This strongly curved interface suggests a truly 3D behavior. A very different behavior is expected for the 2D ribbon-shaped interface that is obtained in thinner samples [19]. To evaluate the differences between the two cases, we also performed fully 2D simulations [START_REF] Ghmadh | Étude par la méthode du champ de phase à trois dimensions de la solidification dirigée dans des lames minces[END_REF]. They showed that comparable oscillations also appear in 2D but their amplitude is weaker, they rapidly damp out and the scaling exponent of relation ( 11) is about -3/4 instead of -3/2.

35µm

FIG. 3: Side view of the numerical domain (35µm in thickness) showing one of the two oscillating fingers intersected by the yz plane that passes through its tip. The tip radius of curvature is roughly equal to 11.5µm here.

IV. EXISTENCE DIAGRAM AND OSCILLATION PERIOD

Varying both the pushing velocity V and the spacing λ in the simulations allowed us to observe a number of occurrences of the 2λ-O mode. (See the movie in [25] for an example of oscillations simulated at V = 9.0µm.s -1 and λ = 77.5µm.) All the occurrences of the 2λ-O mode are reported in the existence diagram of figure 4, together with experimental points. Whereas an excellent agreement is obtained regarding the velocity range, the experimental modes span over a wider range of spacing with, in the central part of this range, a gap in which only stable cells are observed.

To better compare experiments and numerics, we stress that numerics provide, for a given set of physical parameters, values of λ smaller than in experiment, as previously noticed in [8] or in [20]. This corresponds in Fig. 4 to an overall shift of the numerical data to the left following which the numerical domain appears as a shift of the large spacings experimental domain. However, no similar shift of the experimental domain for small spacing is displayed on the numerical data of Fig. 4. This may be traced back to the impossibility of simulating small spacings owing to the much shallower grooves in numerics than in experiment (compare Figs. 7 and8 with Fig. 1). As studied in details for the one-sided phase-field model in [39], this groove feature is found to induce a rapid cell elimination that preempts any possible emergence of durable oscillations. Lastly, we also notice that, following the sensitivity of stability bands to crystalline anisotropy [19], its slight variation between experiment and numerics may yield noticeable implication on the stability diagram. Considering these differences between experiment and simulation, we may finally view Fig. 4 as providing a global agreement between them regarding the existence domain of the 2λ-O mode.

To provide a reliable measurement of their period, oscillations must last several cycles with nearly constant amplitude, no transition to doublet and no perturbation from neighbors. For these reasons, only a part of the experimental events displayed in figure 4 provided data on oscillation periods. Similar restrictions prevented to measure the period of the numerical occurrence at V = 4µm/s. In ref. [7], a large set of experimental results obtained for temperature gradient G = 110 or 140K/cm yielded the relation

T = AV -3/2 , ( 11
)
where A is a constant. As this relationship does not depend on the thermal gradient in this short range, we have addressed its relevance by completing experimental data with few periods obtained at G = 140K/cm. As can be seen in Fig. 5a, both experimental and numerical data clearly follow this power law relationship within numerical accuracy. Fitting them to Eq. 11, one obtains A exp = (2.46 ± 0.08) × 10 3 µm 3/2 s -1/2 and A num = (2.73 ± 0.08) × 10 3 µm 3/2 s -1/2 . The scaling law (11) agrees with that provided by the 1λ-O instability [40] both regarding the exponent and the prefactor which, for a relative spacing of 1/2, amounts to

A = 2π(2k) 1/2 (1 + k) -1 D/V 1/2 c
, i.e., 2.90 × 10 3 µm 3/2 s -1/2 here. We also note that the same relationship between period and velocity, with a rather similar value of coefficient A, was also obtained for 2π/3 oscillation modes in extended 3D samples, although both the oscillation mode and the solute were different there [8][9][10].

To address the dependence of the oscillation period on the cell spacing λ, we have looked for a relation T (V, λ) in non-dimensional variables based on the following diffusive characteristic variables : the diffusion time τ D = D/V 2 , the diffusion velocity V D = D/λ and the diffu- (11) extends to :

sion length l D = D/V . With T /τ D = b(λ/l D ) α (V /V D ) β , relation
T /τ D = b Pe 0.5 , (12)
Pe denoting the Péclet number Pe = λV /D and b being a non-dimensional prefactor. As can be seen in Fig. 5b, experimental data agree with this extended scaling law with a best fitting prefactor b = 8.1 ± 0.4, but numerical data do not extend over a sufficient range of Péclet number to be conclusive.

V. CONCENTRATION FIELD

Intuitively, one may expect the 2λ-O mode to result from a coupling between neighboring cells but a part of the information required to elaborate a complete description of the underlying mechanism is not easy to get from experiment. The situation is somewhat better in numerics because, in addition to purely geometric characteristics of the cell tip like position z t , velocity żt , and radius of curvature ρ, phase-field simulations also provide the concentration field (see Figs. 6,7, and 8). Thus, in complement to the geometrical information provided by the microstructures, one may expect the missing information to be gained from the spatiotemporal evolution of the concentration field.

A concept that proved useful in quantifying the dynamics of the growth modes is constitutional supercooling [41]. It is commonly defined as a temperature difference but may also be expressed in terms of concentrations,

δc(r, t) = c L (z, t) -c(r, t), (13)
where c(r, t) is the actual concentration at location r and time t and c L (z, t) is the liquidus concentration at the corresponding temperature T (z, t),

c L (z, t) = c ∞ k + G m (z + z s -V t) = c ∞ k + T G m z * , ( 14
)
where z * is defined in (4). Using (3) and ( 5) finally yields the nondimensional undercooling

Ω = - k c ∞ δc = β k z * + 1 + β k U α k -β k ϕ 2 -1. (15)
When, for a given velocity, the thermal gradient G is less than a critical value, the planar front destabilizes and 6: (Color online) Time evolution of the fluctuations (denoted by a tilde) of four different parameters (schematic representation). The fluctuations of a given parameter are defined as its instant value minus its average value and they are normalized to one here. The parameters represented are the solute concentration at a fixed distance above the cell tip, ca, and three geometric parameters of the cell tip: position zt, velocity żt, and radius of curvature ρ. The oscillation time period is denoted by T and labels a to i refer to those used in Figs. 7 and8 cells form. A zone of negative undercooling then builds up in the liquid and Ω goes through a minimum Ω a < 0 at some distance ∆z a ahead of the cell tips, as seen in Fig. 7. Here and in the following, we want to focus on the physical limit where the interface becomes very sharp. To do so, everywhere in the liquid Ω is obtained by evaluating the right hand side of Eq. ( 15) for the far-field value ϕ = -1.0. As α k + β k = 2, one then gets

Ω = β k (z * + U ) (in the liquid). ( 16
)
In the solid, we uniformly set the undercooling to an arbitrary negative constant

Ω = cst < 0 (in the solid). ( 17
)
Figure 7 represents Ω in the xz plane that intersects the cell tips. The simulation is performed for a single pair of cells with periodic boundary conditions along x. The resulting data are then duplicated along x to ease visual inspection by the reader. As detailed below, the main characteristics observed in the experiments (Fig. 1) are recovered here.

The fluctuations of four cell parameters are schematically represented as functions of time in Fig. 6. The three geometrical parameters were introduced above, the fourth, c a is the solute concentration at the altitude ∆z a above the cell tip where Ω is minimal. One can follow their evolutions by focussing attention on the central cell in the successive panels of Fig. 7. At initial time (panel a), żt is maximum, ca and zt are zero and increasing and ρ increases after having reached its minimum about a picture earlier (panel h). Taking the tip velocity fluctuation żt as the phase reference, the initial phase shifts remain constant in time, that is roughly π/2 for ca and zt , and 3π/4 for ρ.

At time t = T /2 (panel e) the tip velocity is minimum, so that the vertical solute flux rejected at the tip is lower than average. As a result, at later times, the solute concentration c a ahead of the tip decreases, as shown in panels f, g, h, and i in Fig. 7, producing the large black area at a distance ∆z a above the cell tip (panels f and g). As ∆z a 40µm, the vertical diffusion time ∆z 2 a /D 1.2s is rather short as compared to the oscillation period T 108s. The effects of żt on the dynamics of c a are thus almost immediate.

The counteracting mechanism that re-equilibrates the ∂Ω/∂x concentration in the dark stripe around altitude ∆z a above the cell tips comes from horizontal solute diffusion currents. In the far-field limit of Eqs. ( 16) and ( 17), we obtain

∂Ω ∂x = β k ∂U ∂x (in the liquid), (18)
and

∂Ω ∂x = 0 (in the solid). ( 19
)
Note that this quantity is simply proportional to the horizontal component of the concentration gradient ∂c(r, t)/∂x, or to the opposite of the horizontal concentration current j x = -D(∂c/∂x). Figure 8 represents the time evolution of ∂Ω/∂x over an oscillation period. When c a is minimum (panel g) the horizontal solute currents coming from the two neighboring cells reach maximum extensions and converge to increase the solute concentration above the tip of the central cell (panel g). This mechanism definitely supersedes the coupling between first-neighbor cells but, since the central cell is not totally screened, it still survives for the next period. One can see that the two converging horizontal currents cooperate to pinch the tip of the central cell. As a consequence, ρ decreases to reach a minimum value on the next panel (panel h). This increases diffusion currents at the tip and thus its velocity żt . In the meantime, the situation is the exact opposite for the two cells neighboring the central one: the tip radius increases, so that the tip velocity decreases. This completes the explanation of the way the central cell manages to survive.

As compared to the time analysis just detailed that starts from t = T /2 (panels e to i), the sequence is precisely reversed half a period earlier, as shown by the panels starting from t = 0 (panels a to e). Thus a unique mechanism accounts well for the periodic nature of the 2λ-O mode, with a doubling of the spacing. The following section is devoted to model it.

VI. EMERGENCE OF OSCILLATIONS

The specific features of the 2λ-O instability are its oscillatory character and its wavelength, twice the cell spacing λ. The former feature calls for understanding why a cell that is on the way to be eliminated succeeds in reversing its dynamics to get back to the mean cell position. The latter feature suggests that the same mech-anisms that render unstable the oscillatory modulations at the scale of a cell pair leave those at a single cell scale stable. Our objective here is to identify the main mechanisms responsible for these features and to model their interplay so as to uncover the emergence of oscillatory states.

In terms of temporal eigenvalues of perturbation modes, σ = σ r + iσ i , this requires to find out the mechanisms yielding σ r (V, λ) to change sign or σ i (V, λ) to become non-zero. The former event deals with stability and the latter with the occurrence of oscillations (damped or amplified). Here we shall focus on the latter one since our concern points more to the emergence of oscillations than to stability.

To date, the neutral stability of the 2λ-O mode has been determined in a symmetric model by Floquet analysis at given surface tension and anisotropy [11] and a similar instability called vacillating-breathing instability has been evidenced at large solidification velocity by simulation of a symmetric model with constant miscibility gap [12]. However, to our knowledge, no theoretical modeling of the instability mechanism of the 2λ-O mode has been proposed so far. This contrasts with the case of oscillations of individual cells, i.e., with the 1λ-O mode, where various theoretical [40] or analytical and numerical results [42,43] have been obtained, with controversial conclusions however. In comparison, the major difference that will be brought about by the 2λ-O mode here will be the coupling between neighbor cells.

In this context, the ambition of our modeling will be more qualitative than quantitative : identifying the physical mechanisms responsible for the oscillations, modeling them and their interplay so as to uncover cell oscillations and their main features. Our analysis will rely on the major outcomes pointed out by simulation : the coupling between cells and the interplay between cell form and cell position. It will then aim at expressing them on a rational basis, up to modeling approximations and thus up to quantitative inaccuracies.

In practice, our model will assume quasi-steady states, Laplacian dynamics close to cell tips, a 2D space, negligible surface tension, and a given family of cell forms with infinite grooves. It will explicitly consider the transversal fluxes between cells as the major factor of the 2λ-O instability and will provide both their expression and their implication on cell form and cell tip position.

The model will succeed in uncovering the emergence of oscillations, i.e., of non-zero σ i . However, the velocity at which oscillations occur as well as the stability of the state on which they emerge will differ from those determined in experiment or simulation. In particular, oscillations will appear on an unstable state rather than on a stable one, so that their occurrence will correspond to a transition between 2λ-S and 2λ-O instabilities [11]. Beyond these differences with observations, the mechanisms that generate non-zero σ i and thus oscillations may be expected to be responsible for those of the 2λ-O mode. In particular, the eigenvalue σ i will yield a rela-tionship between period and velocity in close agreement with that evidenced in experiment and simulation.

We now build the model on the following ingredients revealed by experiment and simulation :

(i) two dynamical cell variables : the cell tip position z t and the cell tip curvature radius ρ (Figs. 1 and6).

(ii) a major dynamical factor : the transversal solute flux between cells. It is correlated to the non-dimensional undercooling Ω (15) and drives the modulations of both z t and ρ (Figs. 78).

We respectively label the interface and the cell tips by the indexes I and t. We consider a central cell and we index its left and right neighbors by "-" and "+" respectively. The origin of coordinates (x, z) is placed in the middle of the central cell and on the solidus. Time derivation will be denoted by a dot.

A. Regimes and assumptions

Following the experimental and numerical results, the velocities that will be considered are weak enough for making the Péclet number Pe = λV /D significantly smaller than unity : Pe << 1. For instance, for a cell spacing λ = 80µm, Pe = 1 is reached at a velocity V D = D/λ ≈ 17µm/s, well above the studied range V ≤ 10µm/s (Fig. 4).

We then consider the tip region defined as the cell domain closer to the tip than half a cell spacing λ. There, as variations of solute concentration involve the length scale λ, the low Pe condition yields the advection term V.∇c to be negligible compared to the diffusion term D∇ 2 c. In steady states, the dynamics of solute concentration then yields the concentration field to be quasi-Laplacian, λ 2 ∇ 2 c/c ∞ = O(Pe) and the variation of solute concentration over the cell width to be negligible in comparison to the variation c t -c ∞ undergone over the liquid phase [40].

Regarding the dynamics of the 2λ-O mode, we notice that its period T is actually large compared to the diffusion time τ D = D/V 2 : (T /τ D ) 2 = V /V T with V T = D 2 /A 2 = 0.25µm/s. As the instability domain starts at V ≈ 4µm/s (Fig. 4), this yields T /τ D > 4 following which a quasi-steady approximation for the dynamical states may be assumed. According to it, we may invoke an exponential relaxation of the concentration field in the z direction over a diffusion length lD = D/V I based on the instantaneous growth velocity

V I = V I • e z = V + żt .
Considering a two-dimensional growth in the plane (x, z), the above approximations then allow the concentration field to be approximated as :

c(x, z, t) = c ∞ + (c t -c ∞ ) exp[-(z -z t )/ lD ] + o(c t -c ∞ ).
(20) For simplicity, we shall work within a prescribed family of forms involving a circular tip of radius ρ, a constant width l = 2ρ and infinite grooves (Fig. 9). In the oscillatory state, the cell form will then drift within this family with variable ρ or . This is analogous to the cell evolution considered in the analysis of the 1λ-O mode [40] except that the family of forms was that of Saffman-Taylor fingers and no coupling between cells was in order. In particular, as in [40] and in agreement with the quasisteady approximation, we shall not consider the effects of the groove modulation implied by the variation of cell width λ(t).

The last assumption will be to neglect the effect of surface tension on cell tip undercooling. In particular, in the same set-up, a study of tip position and of tip temperature evidenced a linear evolution of the tip undercooling = (c t -c ∞ )/∆c with 1/V , where ∆c = c ∞ (1 -k)/k denotes the miscibility gap [44]. Assuming a zero undercooling at infinite V then yields a relationship, = V c /V , that is similar to the Bower-Brody-Flemings (BBF) criterion [45]. On this basis, the capillary correction (|m|∆c d 0 /ρ < 3.10 -3 K) to the drop of melting temperature (|m|c t > 3.5K) appears to be of order of 10 -3 . In addition, as the oscillation amplitudes of tip curvature radius ρ and of tip position z t are of the same order, the oscillations of tip temperature induced by ρ appear to be less than a percent of those induced by z t .

These estimates thus legitimize our assumption to neglect capillary corrections to the cell tip undercooling in both steady and oscillatory states. Accordingly, the Gibbs-Thomson relationship at the cell tips reduces to the liquidus relationship :

c t = c ∞ k -(z t -z p ) G |m| , (21)
where z p denotes the position of a steadily growing planar front. For the sake of simplicity, hereafter the origin of the z-axis will be placed there : z p ≡ 0.

B. Dynamics of tip position

We now consider a box of width λ and height L larger than l D , extending from the tip of the central cell towards the liquid phase. We wish to follow the mean solute concentration in it (Fig. 9). Its evolution is linked by solute conservation to the concentration fluxes over the box boundary:

dC dt = d dt V (c -c ∞ ) dx dz = - Γ j • n dl, (22)
where C denotes the net excess of solute in the box beyond the concentration c ∞ , V and Γ the volume and the contour of the box, n its outward normals and j the concentration fluxes.

On the contour Γ, the fluxes to consider are at z = z t + L, the advective flux j ∞ = -c ∞ V I e z , the diffusion fluxes being negligible; at z = z t : the advective flux j t and the diffusion flux j d induced by solute rejection at the interface; at x = ±λ/2, the lateral horizontal fluxes j -and j + generated by the neighboring cells. As the concentration field is quasi-Laplacian in the tip region, the divergence theorem may be invoked to express the integral of the diffusion flux j d = -D∇c on the line z = z t from that on the interface I where, by solute conservation, -D∇c • n I = (1 -k)c I V I • n I . One obtains [28]:

- z=zt j d • n dx = D G |m| (λ -) + (1 -k)c t V I , (23)
where denotes the cell width and n = -e z the outward normal to the contour Γ at z = z t . Following the exponential relaxation (20), one may approximate the excess of solute concentration in the box by C ≈ λ lD (c t -c ∞ ). Similarly, the advection flux j t at z t may be approximated by by its value j t = -c t V I e z at the cell tip and the integrals of the lateral fluxes over the lateral sides of the box as lD times their value j ± (z t ) at z t . These fluxes are driven by the differences of concentration δc ± = c T± -c M± between the tip of a neighbor cell (point T ± in Fig. 9) and the center of the central cell at the same height z t± (point M ± in Fig. 9). Fluxes may thus be approximated from the resulting mean concentration gradient δc ± /λ as : j ± (z t ) ≈ ∓D δc ± /λ e x .

Following the liquidus relationship (21), the cell tip concentrations c T± equal the liquidus concentration c L (z t± , t) at the altitude z t± : c T± = c L (z t± , t). Evaluating c M± from (20) then yields, at first order in (z t± -z t ) and żt ,

δc ± = [c t -c ∞ -l D G/|m|](z t± -z t )/l D or δc ± = ∆c ( -V c /V ) (z t± -z t )/l D .
Accordingly, the concentration differences δc ± , and thus the resulting fluxes j ± (z t ), depend on the deviation ( -V c /V ) to the BBF criterion which is known to be small here [44]. Our simulations provide however the opportunity of determining this deviation since, following c T± = c L (z t± , t), it appears that δc ± = δc(0, z t± , t) as defined in (13) or δc ± = -Ω(0, z t± , t) c ∞ /k according to [START_REF] Georgelin | Étude expérimentale des instabilités secondaires en croissance cristalline directionnelle[END_REF]. In particular, noticing that Ω(0, z t , t) = 0 and linearizing Ω with respect to the shift (z t± -z t ) of cell tip positions yields

δc ± = -(c ∞ /k) (∂Ω/∂z)| zt (z t± -z t )
and, finally :

j ± (z t ) = ± D λ c ∞ k ∂Ω ∂z zt (z t± -z t ) e x . ( 24
)
Knowing the concentration fluxes, we may now express relation ( 22 As, in our velocity range, |δ| stands below 1/4, we shall neglect its contribution in the following. This turns out neglecting the variations of the diffusion length lD and thus the internal dynamics of the solute layer. One then obtains from (22):

żt = τ -1 D [ (1 -/λ)(l T -z t ) -k z t /λ] -(1 -/λ)V I + D 2 λ 2 V c (1 -k) -1 ∂Ω ∂z zt (z t+ + z t--2z t ), ( 25
)
where τD = D/V 2 I . This relation shows the effect on the tip velocity of both the cell form (i.e., the cell width or the tip radius ρ = /2) and the lateral fluxes induced by a shift (z t± -z t ) of tip position with respect to neighbor cells. It calls for completing the dynamical system with a dynamical equation for the tip curvature radius ρ.

C. Dynamics of cell form

For isolated dendrites growing in a homogeneous melt, a solvability condition prescribes the tip curvature radius ρ in steady and quasi-steady growth states : ρ 2 V ∝ γ/c t where γ denotes surface tension [46]. This relation can be extended to dendritic arrays provided that dendrites are separated enough for not being noticeably influenced by their neighbors [47,48]. This restriction requires their spacing λ to be large compared to the diffusion length l D = D/V , and thus their Péclet number Pe = λV /D to be large : Pe >> 1.

Here, the cellular pattern stands in the opposite regime, Pe = O (1), where cells are close enough to significantly interact by diffusion. Then, their steady form is largely dependent on the cell spacing λ [49]. It then does not refer to a solvability condition which, in particular, involves no dependence on λ. For small Péclet numbers, Pe 1, and a uniform cellular pattern, steady cells correspond to Saffman-Taylor fingers [49] whose relative width l/λ is selected by surface tension [50,51]. However, this kind of form ceases to be valid when lateral fluxes of solute are induced by neighboring cells. We shall then model the dynamics of form evolution by considering two factors : (i) the direct effect of transversal fluxes on the interface position ; (ii) the relaxation of perturbed forms towards the steady form.

(i) Lateral fluxes make the solid phase melt, as modeled in (25) for the tip position dynamics. However, this melting may be differential, yielding an evolution of form, i.e. of ρ. We model this effect by considering the evolution of concentration on the side of a cell (Fig. 9). As we restrain attention to the tip vicinity, fluxes may be taken as constant and equal to j ± (z t ). We then look at an altitude z I for solute conservation in a band of liquid of height δz and length λ/2 -|x I | where x I denotes the interface abscissa. Assuming homogeneous solute concentration, we obtain the evolution of concentration ċI = ∓j ± (z t )/(λ/2 -|x I |). Invoking the liquidus relationship ( 21) at the interface,

c I = c ∞ /k-z I G/|m|, yields the evolution of interface po- sition, żI = ±j ± (z t )/(λ/2 -|x I |) (|m|/G). Using ρ -1 = -∂ 2 z I /∂x 2 I | x I =0
, we now deduce from it the evolution of the tip curvature radius : ρ = ±16 (|m|/G) j ± (z t )ρ 2 /λ 3 .

(ii) The basic stability of the cellular form in absence of transversal fluxes goes together with a natural relaxation. As this relaxation is driven by diffusion, it involves a diffusive time scale τ ρ driven by the relevant length scale, the steady state curvature radius ρ : τ ρ ∝ ρ2 /D. For simplicity, we shall assume that the relaxation time scale reduces to the diffusive time scale, thus yielding the following contribution to the tip curvature radius dynamics: ρ = -D/ρ 2 (ρ -ρ).

Using (24), the two above contributions altogether yield the following evolution equation for the tip curvature radius:

ρ = 16 D 2 λ 2 V c ρ 2 λ 2 (1 -k) -1 ∂Ω ∂z zt (z t± -z t ) - D ρ2 (ρ -ρ). ( 26
)
D. Coupled dynamics and 2λ-0 oscillations Altogether, relations ( 25) and ( 26) provide, with = 2ρ, a coupled dynamical system for z t and ρ. Labeling by a tilde the fluctuations zt = z t -zt , ρ = ρ-ρ where zt and ρ denote average values, their linearized dynamics may be straightforwardly derived (see Appendix). It reads

żt = az t + bρ, ( 27
) ρ = cz t + dρ, (28)
where, in case of phase opposition between a cell and its neighbors, zt± = -zt , and

a = (λ/ ¯ + k -1) ( V 2 -V 2 D ), ( 29
) b = -2k (λ/ ¯ ) 2 λ/ ¯ + k -1 V λ ( V V c -1), ( 30
) c = -32 D 2 λ 2 V c ρ2 λ 2 (1 -k) -1 ∂Ω ∂z zt , (31)
d = - D ρ2 , (32)
with

V 2 = -4 λ/ ¯ λ/ ¯ + k -1 D 3 λ 2 V c (1 -k) -1 ∂Ω ∂z zt . ( 33
)
Transverse fluxes strenghten the V 2 term and parameter c. We stress that, as ∂Ω/∂z(z t ) is negative (Fig. 7), both of them are intrinsically positive. On the other hand, a is positive for V < V and negative above, while b is negative above onset V > V c and d is always negative.

The eigenvalues σ 1 and σ 2 of the linear system ( 27 8 : λ = 77.5µm, ρ = 27µm, ¯ = 2ρ, ∂Ω/∂z(z t ) = -3.75 × 10 -3 µm -1 that we assume independent of V , and V = 75.6µm/s. It shows the occurrence of a nonzero imaginary part of eigenvalues for V > V O . Below V O , eigenvalues are real with one of them positive, so that the cellular state is unstable with respect to the 2λ-S instability [11]. This mode, which drives an elimination of one cell over two, is reminiscent of the spatial perioddoubling instability of dendrites [47,52].

The facts that the cellular pattern on which the oscillatory instability develops is actually stable and that V O = 18.0µm/s is large compared to observations, indicate that the modeling is too crude to grasp the correct behavior of the real part of eigenvalues. Fixing this caveat stands beyond the scope of this study. This would require a more refined model that better addresses cell geometry, cell undercooling, transversal fluxes and the pulsation of the boundary layer implied by cell oscillations. In particular, the latter is an important ingredient of the 1λ-O instability [40] that has been overlooked here by neglecting the factor δ in the derivation of (25). On the other hand, determining more accurately the relaxation time scale τ ρ and its dependency on velocity would be important for dealing with stability. Finally, as stressed by Kessler and Levine on the 1λ-O mode [42], the interaction between tail and tip, which is overlooked here, might be relevant to cell stability. Beyond these shortcomings, the occurrence of complex eigenvalues (σ i = 0) in the actual modeling succeeds in pointing out how transversal solute fluxes driven by differential cell positions can trigger oscillations.

It is worth highlighting the essential role of cell coupling by transversal fluxes on both oscillation and stability. Regarding oscillations, these fluxes sustain the parameter c that is essential for having negative D and hence oscillations. In particular, with no fluxes, c would vanish and yield positive D and thus real eigenvalues. This situation refers in particular to 1λ-modes where, as all cells undergo the same evolution (z t± = z t ), transversal fluxes vanish by symmetry. The model is thus unable to capture 1λ-O oscillations, presumably because of its neglect of the boundary layer pulsation. Regarding stability, fluxes provide a positive contribution to parameter a from V 2 , which supports instability. In particular, we note that, for a single oscillating cell surrounded by steady neighbors, i.e., for zt± = 0, the flux amplitudes would be divided by two. Then the same modeling would turn unstable oscillations into damped oscillations.

From a more physical point of view, the mechanisms yielding oscillations may be described this way. Consider a positive modulation of curvature radius at same tip positions : ρ > 0, zt = 0, Fig. 7e. As the flatter the tip, the weaker the diffusion flux around it, the resulting depletion of diffusion on the central cell yields a negative tip velocity modulation : żt < 0, zt < 0, ρ > 0, Fig. 7f. The difference of tip position between the cell and its neighbors then drives lateral fluxes which make the cell melt at its tip and on its sides. Tip melting raises the strength of the negative modulation of velocity żt < 0, possibly yielding cell elimination by the 2λ-S instability. Sides melting yields the tip radius to decrease, in addition to its natural tendency to relax towards its mean value. When this effect has enough time to develop, it may yield the modulation of curvature radius to reverse sign : ρ < 0, zt < 0, Fig. 7g. Meanwhile, thanks to a resulting stronger diffusion at the cell tip, the tip velocity raises. When this effect is stronger than that induced by tip melting, it may even reverse the velocity modulation, żt > 0, zt < 0, ρ < 0, Fig. 7h, and yield the cell to recover the neutral position but with a positive velocity and a thin curvature radius : zt = 0, żt > 0, ρ < 0, Fig. 7i. The remaining half cycle may then resume and complete the oscillation.

Accordingly, transversal fluxes between cells generate both tip melting and side melting. Whereas cell elimination is induced by tip melting, oscillations result from side melting which can make a cell drifting down the grooves succeed to recover. The dynamical competition between both these effects either yield the 2λ-S instability, the 2λ-O instability or cell stability (Fig. 10).

FIG. 11: Evolution with the pushing velocity V of the period Ti = 2π/|σi| of the 2λ-O mode at λ = 77.5µm, following the imaginary part of the eigenvalue σi provided by the system (27) (28). To compare with the experimental and numerical trend (11), the period Ti is multiplied by V 3/2 . The nearly constant value reached shortly after the occurrence of oscillations at VO = 18µm.s -1 agrees with the observed relationship (11). The insert shows data in log-log coordinates and a comparison with the trend Ti ∝ V -3/2 of relation (11).

On a quantitative ground, the period T i of oscillations is provided by the imaginary part σ i of eigenvalues :

T i = 2π/|σ i |.
Its evolution with V at λ = 77.5µm, is compared in Fig. 11 to the relationship (11) found in experiment and simulation. Apart from a sudden decrease at the onset of oscillations V O , the combination T i V 3/2 quickly reaches a nearly constant value of about A ≈ 0.4 × 10 3 µm 3/2 s -1/2 , of the same order as the values found in experiment (A exp = 2.46 × 10 3 µm 3/2 s -1/2 ) or simulation (A num = 2.73 × 10 3 µm 3/2 s -1/2 ). This agreement on both the type of power law and the order of magnitude of its prefactor gives additional credit on the relevance of our model regarding oscillations.

VII. PARAMETRIC EXCITATION

The 2λ-O instability corresponds to a spontaneous destabilization of the cellular state towards a limit cycle and thus to the natural emergence of an oscillating pattern by a Hopf bifurcation. A canonical way to explore the vicinity of this bifurcation and reveal the main features of the underlying oscillator consists in forcing oscillations by an external, resonant, mean. Below, we use this mean to explore some features of the 2λ-O oscillator by addressing, both experimentally and numerically, its response to a parametric excitation.

Extrinsic excitations are induced here by modulations of the pushing velocity. This corresponds to uniform modulation whose spatial signature thus differs from that of the 2λ-O mode. It may yield an additive or a multiplicative forcing depending on whether the modulation affects the absolute or the relative value of the main variable, here the solute concentration. We show below that solute advection by this velocity forcing yields a multiplicative forcing which is known to possibly excite oscillatory modes when parametric resonance is in order.

We thus introduce, both in experiment and simulations, a sinusoidal modulation of the pushing velocity:

V f (t) = V 0 + V 1 (t) = V 0 + V 1 cos(ω f t), ( 34
)
where V 0 is the average pushing velocity, V 1 the excitation amplitude, ω f the forcing pulsation and T f = 2π/ω f the forcing period. Regarding the experiment, we stress that the forcing period will be larger than the relaxation time of the thermal field, τ = g 2 /κ = 50s, g = 5mm being the gap between heaters and coolers and κ = 5.10 5 µm 2 .s -1 the thermal diffusivity of the sample plates. Accordingly, we may consider that the velocity modulation negligibly affects the thermal field [24]. However it actually modulates the advection of solute concentration in the thermal field frame, V f (t)∂ z c, with noticeable consequences. The sinusoidal forcing of solute advection, V 1 (t)∂ z c, splits into the advection of the time-independent solute field c and that of the remaining fluctuating part c. The former term, V 1 (t)∂ z c, yields an additive sinusoidal modulation to the dynamics of concentration which simply drives an homogeneous oscillation of solute concentration and thus of the whole cellular pattern, at the same period T f . The latter term, V 1 (t)∂ z c, corresponds to a multiplicative forcing term whose time-dependence depends on that of both ∂ z c and V 1 (t). It then drives a parametric excitation of the system from which the oscillatory 2λ-O mode may emerge, in case of resonance.

Looking at a cell as an oscillator, damped or amplified, of angular pulsation ω = 2π/T , T being the period of the 2λ-O mode, it is instructive to determine the modulation of ω induced by the velocity modulation (34). For this, we notice that Eq. 11 still applies since the modulation periods T f that will be considered hereafter are long as compared to the diffusive and advective relaxation times λ 2 /D < 9s and τ D = D/V 2 ≈ 28s. In terms of ω, Eq. 11 then reads ω 2 = BV 3 with B = 4π 2 /A 2 a constant. Plugging the forced velocity V = V f (34) in this last relation yields ω(t) 2 = B V 3 0 [1 + ν cos(ω f t)] 3 with ν = V 1 /V 0 and, as ν will be taken much smaller than unity here,

ω(t) 2 ω 2 [1 + 3ν cos(ω f t)] ( 35
)
up to first order in ν. This clearly evidences, in the framework of oscillators, the parametric excitation of cells induced by V 1 (t) in terms of a modulation of the pulsation of the 2λ-O mode.

Given an oscillator of natural period T , the condition for parametric resonance at the lowest forcing period T f is T f ≈ T /2. Applying this to the velocity modulation (34) will first yield a spatially homogeneous oscillation of both the concentration field and the interface position at the forcing period T f , in response to the additive modulation brought about by concentration advection. Moreover, if parametric resonance occurs, cell oscillations should appear at period 2T f ≈ T , with phase relationships between cells to be determined.

To better highlight parametric resonance, velocity modulations are induced, both in experiment and in simulation, on steady states that stand within the stability domain of the 2λ-O mode but that are close to the instability frontier. As expected, these modulations are found to induce first global oscillations of cells at the forcing period T f and, for T f ≈ T /2, cell oscillations similar to those of the 2λ-O mode with a period close to 2T f ≈ T and phase opposition between neighbor cells (Figs. 12,13,14). While the period doubling between excitation at T f and oscillations at 2T f results from parametric instability, we notice that the breaking of spatial uniformity when passing from the homogeneous modulation to the phase opposition between oscillating cells does not result from the forcing. It then reveals the spatial nature of the underlying oscillator, the 2λ-O mode.

Regarding experiment, parametrically induced cell oscillations are displayed over a cycle in Figs. 12 and13. Figure 12 shows phase opposition between neighbor cells through spatial alternation of larger and thinner cells on an isotherm. Over a cycle, each cell alternates between thin and fat and upward or downward motion. As oscillation amplitudes are faint, we complete observations by figure 13 that shows a span of cells over their whole groove. On each snapshot, the cell widths are almost the same at cell tips but significantly differ at cell tails. They thus encode the oscillation of cell tips over the grooves, i.e. here, over a time interval of about T /3.

Regarding simulations, Fig. 14 displays, for a given basic state and various forcing periods T f , the oscillations of the lateral widths of two adjacent cells. At the beginning of the parametric modulation, one notices, for all forcing periods, in-phase oscillations of cells at T f , in response to the additive forcing induced by the velocity modulation. However, some periods later, it also appears oscillations at the double period 2T f that are in phase opposition between neighbor cells. These oscillations correspond to a parametrically excited 2λ-O mode. We signal them by double arrows in figure 14.

Occurrence of parametric oscillations is especially manifest for T f = 115s (Fig. 14a) where, at t ≈ 2700s, the cell signals grow out of phase to reach robust oscillatory states displaying a clear phase opposition. This evolution gathers both the temporal feature of a parametric resonance (period doubling) and the spatial feature inherent to the 2λ-O mode (phase opposition between adjacent cells). It thus consists of a mix between parametric resonance and 2λ-O instability. In phase space, this corresponds to a basic state which, after destabilization by parametric excitation, restabilizes onto a limit cycle akin to the 2λ-O oscillatory mode.

On the other hand, at a smaller forcing period T f = 110s (Fig. 14b), the maxima of the signals show clear oscillations at 2T f , but their minima involve local maxima. This indicates that oscillations at T f are still noticeably present. When these local maxima raise at t ≈ 1400s, t ≈ 2800s and 3800 < t < 4200s, they then make the oscillation at T f dominant, meaning that global cell oscillations have overwhelmed the 2λ-O like oscillations. Oscillations at 2T f of both the maxima and the minima of the signals are then only noticeable on limited time intervals. Parametric resonance thus still occurs, but by intermittency.

At still smaller forcing periods, T f = 90s (Fig. 14c) and T f = 80s (Fig. 14d), parametric resonance fails to clearly appear. However, oscillations at 2T f of the signal maxima are still noticeable for limited time intervals of several periods in Fig. 14c and of even reduced duration in Fig. 14d.

Altogether, these dynamics reveal both a robust parametrically excited 2λ-O mode for T f = 115s (Fig. 14a) and an intermittent evolution towards this limit cycle, with an intermittent character that grows when T f decreases (Figs. 14b,c,d). The latter feature signals an increased complexity of the phase space which still includes a destabilized basic state by parametric excitation and an attracting limit cycle, but also additional structures responsible for this intermittency.

Table III sums up the experimental and the numerical conditions of the observed parametric oscillations. Here, one considers parametric oscillations in a broad sense, i.e. not only for robust 2λ-O modes as in figure 14a but also for those that show up the characteristic features of the 2λ-O mode intermittently (Figs. 14 b,c,d). The upper and lower bounds of the corresponding forcing periods T f determine a resonance band roughly centered around half of the natural period T of the unforced 2λ-O mode at the same spacing and velocity.

Figure 15 reports the way the oscillation amplitudes vary with the forcing period T f in simulations. A scan of T f has been performed from 70s to 160s but parametric oscillations have only been found in the range 90 -125s. As the oscillation amplitudes are tiny, of the order of few micrometers, error bars are relatively large. In experiment, the oscillation amplitudes are of the same order but the resolution is even lower, which prevents to obtain relevant amplitude variations.

In both experiments and simulations, parametric oscillations have displayed a band of resonance [T f -, T f + ] that is reported in Table III together with its bandwidth ∆T f = T f + -T f -and the resonance bandwidth predicted for a linear oscillator, δT f = 3νT /4 [START_REF] Landau | Mechanics Course of Theoretical Physics[END_REF]. For both experiments and numerics, the former is about ten times larger than the latter. We aim below to relate this discrepancy to the intermittency noticed on the cell signals and to the expected non-linear nature of the 2λ-O oscillator.

To model the parametric excitation of the 2λ-O mode with minimal complexity, one may overlook its spatial features to focus attention on a dynamical system involving parametric modulation around a Hopf bifurcation. The Hopf bifurcation signals that the eigenvalues λ = σ ±iω of the system linearized around the fixed point (0, 0) turn their real part σ positive at some threshold value of a control parameter, here the pushing velocity V . This yields the following linear modeling in the vicinity of the basic state : z-2σ ż+(σ 2 +ω 2 )z = 0. Parametric excitation is introduced as a temporal modulation ω(t) of the pulsation ω following which the linear dynamical system is then represented by a Mathieu's equation. Its standard analysis shows that parametric resonance occurs for forcing pulsations ω f belonging to a band centered on 2ω, whose width extends over δω f = 3νω when no damping is in order and on a more restricted and shifted domain otherwise [START_REF] Landau | Mechanics Course of Theoretical Physics[END_REF]. This Arnold tongue corresponds to forcing periods T f centered on T /2 and extended over a bandwidth δT f = 3νT /4. Accordingly, the linear modeling recovers parametric resonance but on a bandwidth much smaller than ∆T f . It thus explains the emergence of a robust 2λ-O mode (i.e. Fig. 14a, T f = 115s) with a shift (T f -T /2) linked to damping and anharmonicity of the oscillator. However, it fails to uncover the intermittent occurrences of this mode beyond (Figs 14 b,c,d). Non-linearities must thus be considered for explaining them.

Introducing non-linearities reveal a far more complicated behavior involving, for damped oscillations (σ < 0), a parametric resonance and, for amplified oscillations (σ > 0), a combination of parametric resonance and of Hopf-induced oscillations [START_REF] Rand | [END_REF]. In particular, in a kind of Poincaré map, the phase space shows in the former case a fixed point destabilization towards non-trivial stable fixed points corresponding to parametric oscillations and, in the latter case, a more complex structure involving stable and unstable fixed points together with stable and unstable limit cycles induced by the Hopf bifurcation [START_REF] Rand | [END_REF]55]. These limit cycles then modulate the parametric oscillations yielding quasiperiodic evolutions.

Here, the basic state of parametric forcing was initially taken in the stable domain of the 2λ-O mode, close to the stability limit. However the additive forcing induces oscillations of the growth velocity that make it enter the unstable domain on some part of the cycle. For this reason, the damping coefficient σ is presumably not a small negative constant but a slightly oscillating parameter which enters the unstable domain of positive values σ > 0. Therefore, not only robust parametric oscillations but also quasi-periodic evolutions resulting from a combination of 2λ-O mode oscillations at period T and of parametric oscillations at 2T f may be expected. In addition, a sensitivity to perturbations may result from the existence of unstable fixed points and unstable limit cycles in the Poincaré map of the model.

These dynamical features actually reproduce those found in numerical simulations for various forcing period T f (Fig. 14). They highlight the difference between the parametric oscillations in a broad sense that includes intermittency and those in the strict sense observed on the robust oscillations at T f = 115s (Fig. 14a). In particular, the small value of the bandwidth δT f of a linear oscillator compared to the observed bandwidths ∆T f shows the practical relevance of non-linearities for uncovering the intermittent emergence of a parametric 2λ-O mode. We have revisited the emergence and the nature of the 2λ-O oscillatory mode of solidification by combining experimental and 3D phase field numerical studies.

Both agree qualitatively on the behavior of oscillatory cells and on the possibility to parametrically excite the 2λ-O mode. They also both conclude quantitatively to the same bounded range of velocity for the occurrence of oscillations and on the same relationship between period and velocity. As phase field simulations involved the physical parameters of the experiment, this agreement confirms the capability of actual 3D phase field simulations to provide realistic dynamics of solidification interfaces.

Analysis of the concentration fields displayed in simulation revealed noticeable transversal fluxes of solute concentration between a cell and its neighbors. To explore whether they could be an important part of the oscillation mechanism, we have built a simple model of concentration balance in which transversal fluxes led cell tips to melt with implications on both tip position and tip curvature radius. These fluxes depend not only on the temperature difference between cell tips but also on the deviation of the cell tip undercooling from the Bower-Brody-Flemings criterion [45]. Emergence of oscillations were then recovered, the transversal fluxes being mandatory for them to occur. In addition, the model recovered the relationship between period and velocity with the appropriate order of magnitude, thus giving credibility to its relevance. Analysis of the causal links between the phenomena involved in the model then provides a complete description of the mechanism of oscillations. However, some discrepancies regarding the transition from damped to amplified oscillations call for extending the model to a more realistic context.

To explore the nature of the oscillator underlying the 2λ-O instability, we have induced parametric excitations from modulations of the pushing velocity. This homogeneous forcing led, both in experiment and simulation, 2λ-O oscillations to emerge, thus exhibiting a breaking of homogeneity. This parametric excitation enabled the 2λ-O instability to be studied on a broader range extending over its stable domain. It then appeared a complex behavior mixing 2λ-O oscillations with global oscillations. Its occurrence has been linked to the nonlinearity of the underlying oscillator, following which parametric excitations around the corresponding Hopf bifurcation generate a complex phase space and intermittency.

Altogether, these results deepen the actual knowledge on the 2λ-O mode while attesting the relevance of 3D phase field simulations to accurately uncover cell dynamics. More fundamentally, the essential role of cell coupling by transversal fluxes offers a useful insight into the physical origin of multi-cell dynamics and opens the way to further model them on various instances.

As = 2ρ, one gets ¯ = 2ρ and ˜ = 2ρ. Using (A1) (A2) in (A3) with VI = V and ṼI = żt then yields :

żt = - V 2 D (λ/ ¯ + k -1) zt (A5) -4 λ/ ¯ D 2 λ 2 V c (1 -k) -1 ∂Ω ∂z zt zt -2k (λ/ ¯ ) 2 λ/ ¯ + k -1 V λ ( V V c -1) ρ.
Relations (A4) (A5) correspond to the linearized dynamical system (27) (28) with the coefficients (a, b, c, d) reported in section VI D.

FIG. 1 :

 1 FIG. 1: Experimental evidence of a 2λ-O mode. From (a) to (f), six successive snapshots of the same front are shown at time intervals of 20s over a cycle. The pushing velocity is V = 8 µm/s, the image width is 280 µm, the oscillation period 100s and the thermal gradient G = 110K/cm. Snapshots (a) and (f) thus refer to the same phase of the cycle. Oscillations of tip position are tiny but those of cell width and of tip curvature are noticeable. First neighbor cells are in phase opposition and arrows in snapshot (a) indicate the direction of cell velocities in the frame of the thermal gradient.

FIG. 2 :

 2 FIG. 2: Limit cycle of oscillating cells obtained from the differences of cell tip positions δzt and cell widths δλ between two neighbor cells. V = 13µm.s -1 , λ = 45µm, G = 140K/cm.

FIG. 4 :

 4 FIG. 4: Existence diagram of the 2λ-O mode obtained at G = 110K/cm : different symbols correspond to the occurrences found experimentally (full circles) and numerically (empty triangles).

FIG. 5 :

 5 FIG.5: Log-log plot of (a) the oscillation period T vs pushing velocity V and (b) T /τD vs Pe. Full circles for experimental data, empty squares for numerical data. Data are fitted with respect to relation(11) in (a) and to relation (12) in (b). Continuous or dotted line for a fit to experimental or numerical data.

  FIG.

FIG. 7 :

 7 FIG. 7: (Color online) Time evolution of the 2λ-O mode obtained numerically. The pushing velocity is V = 9µm/s and the image width is 310µm. Solid cells appear in black in the lower part of the panels and white arrows indicate the cell tip velocities. In the liquid phase above them, isosurfaces are drawn for the constitutional undercooling Ω defined in the text. A complete oscillation period, T 108s is swept, the frames respectively corresponding to times (a) 0; (b) T /8; (c) T /4; (d) 3T /8; (e) T /2; (f) 5T /8; (g) 3T /4; (h) 7T /8; (i) T .

FIG. 8 :

 8 FIG. 8: (Color online) Time evolution of the field ∂Ω/∂x for the same parameters as those of Fig. 7. The two arrows in panel g indicate the orientations of the horizontal solute currents, proportional to -∂Ω/∂x, that converge toward the central cell. At the bottom and the top of the frames, the zero current areas respectively represent the solid and the far-field liquid.

FIG. 9 :

 9 FIG.9: Sketch of the domain used for modeling. Cells geometries are taken as circular tips of radius ρ = /2 followed by a constant width . The coordinates of a cell interface point are (xI , zI ). The tips of the first neighbor cells are designed by T-and T+. The points in the middle of the central cell that stand at the same height as the neighbor cell tips are called M±. A control volume V of width λ and height L > lD is defined in which solute evolution is determined from the solute fluxes across its sides. In addition to the usual advection and diffusion fluxes, transversal fluxes are induced here by the shift of tip position between a cell and its neighbors. The contour of the box is labelled Γ and its outward normals n. At a finer scale, the evolution of solute concentration in elementary bands δz around the cell tip enables the implication of transversal fluxes on the evolution of tip curvature radius ρ to be determined.

  ) as a dynamical equation for the tip position. Regarding the time derivative of C ≈ λ lD (c t -c ∞ ), we note that dc t /dt = -żt G/|m| and that d lD /dt = -lD zt /V I so that dC/dt = -λ lD G/|m| żt [1 + δ] with δ = (c t -c ∞ )/V I |m|/G zt / żt . The BBF criterion and the amplitude estimate | zt | ≈ | żt |/T then yield |δ| ≈ τ D /T .

  ) (28) are monitored by the trace Σ = a + d = σ 1 + σ 2 , the determinant Π = ad-bc = σ 1 σ 2 of the corresponding Jacobian matrix, and by the discriminant D = Σ 2 -4Π = (a -d) 2 + 4bc of its characteristic polynomial. In particular, emergence of oscillations corresponds to D turning negative. Here D(V ) is positive at V c since b(V c ) = 0. It then decreases with V as a(V ) and b(V ) do and eventually turns negative at some velocity V O at which an oscillatory dynamics emerges.

FIG. 10 :

 10 FIG. 10: (Color online) Evolution of the eigenvalue spectrum with V for modulations at two cell spacings. Dashed lines : real eigenvalues σ± with σ-< σ+ ; full lines : complex eigenvalues with σr and σi the real and imaginary parts. As V increases, a transition from 2λ-S to 2λ-O instability is encountered at VO = 18µm.s -1 . At larger velocity, the 2λ-O instability restabilizes as the real part σr turns negative.

Figure ( 10 )

 10 Figure(10) displays the evolution with V of eigenvalues for the following specific values of parameters corresponding to the simulation inFigs 7 and

FIG. 12 :

 12 FIG.12: Experimental evidence of a parametrically excited 2λ-O mode over a cycle. G = 140K/cm, V0 = 7µm/s, V1 = 0.5µm/s and T f = 141s. The snapshots show cell oscillations in the vicinity of cell tips over an oscillation period T = 282s that is twice the forcing period T f . The delay between them is a quarter of period. Arrows denote the direction of cell velocities and the mean cell spacing is 85µm. Each snapshot shows phase opposition between neighbor cells. These features correspond to those a 2λ-O mode parametrically excited at a forcing period T /2.

FIG. 14 :

 14 FIG. 14: (Color online) Numerical simulations of parametric excitations with a forcing amplitude V1 = 0.2µm.s -1 . The basic state refers to V0 = 7µm.s -1 , λ = 65µm, and the forcing period T f is (a) 115 s, (b) 110 s, (c) 90 s, (d) 80 s. The variables λi refer to cell width of neighbor cells. Double arrows indicate oscillations at the double period 2T f , with a phase opposition between neighbor cells. They correspond to a parametrically forced 2λ-O mode. Figure (a) shows the growth and saturation of a robust parametrically excited mode involving, as the 2λ-O mode, a period 2T f for both λi and phase opposition between them. The remaining figures show that the parametric mode appears more and more intermittently as T f decreases below its value in (a).

2 FIG. 15 :

 215 FIG.15: Peak to peak amplitude of the driven 2λ-O mode as a function of the driving period T f . Average velocity V0 = 7µm/s, velocity amplitude V1 = 0.2µm/s, spacing 2λ = 130µm. Amplitude is maximum for T f = 115s close to T /2 = 100.5s. The shift (T f -T /2) refers to damping and anharmonicity of the oscillator.

TABLE I :

 I Physical parameters of the succinonitrileacrylonitrile alloy of interest and control parameters imposed in the experiments and the numerical simulations. In the first entry, m is the liquidus slope (negative here) and c∞ is the acrylonitrile concentration

	Physical parameter	Symbol Magnitude Unit
	liquidus temperature drop	-mc∞	2.0	K
	solute diffusion coefficient (liq.)	D	1350	µm 2 /s
	capillary length	d0	0.01296	µm
	partition coefficient	k	0.286	-
	anisotropy strength	4	0.011	-
	Control parameter	Symbol Magnitude Unit
	thermal gradient	G	0.011	K/µm
	pushing velocity	V	3.0 -13.5 µm/s
	critical velocity	Vc	2.97	µm/s

TABLE II :

 II Numerical parameters used in the simulations.

	Parameter	Expression Value
	Convergence	W0/d0	90.0
	Grid spacing	δs/W0	1.2
	Noise amplitude	F	0.05

TABLE III :

 III Experimental and numerical parameters for the parametric excitation. The fifth entry gives the range [T f -, T f + ] of forcing periods for which the 2λ-O mode is parametrically excited either robustly or intermittently. ∆T f = T f + -T f -provides the bandwidth of this resonance band and δT f its predicted value for a linear oscillator.

	Parameter	Experiments	Numerics	Unit
	λ	100.0	65.5	µm
	V0	7.0	7.0	µm/s
	V1	0.5	0.2	µm/s
	T /2	70.5	100.5	s
	[T f -, T f + ]	[50.0, 120.0]	[90.0, 125.0]	s
	∆T f	70	35	s
	δT f = 3νT /4	7.6	4.3	s
	VIII. SUMMARY AND CONCLUSION

Acknowledgment We would like to thank Laurent Raymond, Damien Tourret and Florent Perez (supported through A*MIDEX grant number ANR-11-IDEX-0001-02 funded by the Investissements d'Avenir programme of the French Government) for their help to port our codes to GPU-based computers. Région PACA provided financial support through the open project ENEMS ( Étude Numérique et Expérimentale des Microstructures de Solidification).

APPENDIX A: LINEARIZED SYSTEM

We seek to linearize the dynamical system (25) (26) by expanding it at first order in perturbation amplitudes. We denote by a bar v the average values of variables v and by a tilde ṽ their fluctuating part v -v.

We restrict here the analysis to a 2λ mode for which zt± = -z t . As zt± = zt , we obtain :

• at order zero in fluctuations :

which yields the average tip position :

• at first order in fluctuations :