Michael Fleck

Claas Hüter

Denis Pilipenko

Robert Spatschek

E A Brener
email: e.brener@fz-juelich.de

Pattern formation during diffusion limited transformations in solids

Keywords: diffusion, elasticity, growth, solid-state transformation user supplied): Phase field modeling, Boundary integral

We develop a description of diffusion limited growth in solid-solid transformations, which are strongly influenced by elastic effects. Density differences and structural transformations provoke stresses at interfaces, which affect the phase equilibrium conditions. We formulate equations for the interface kinetics similar to dendritic growth and study the growth of a stable phase from a metastable solid in both a channel geometry and in free space. We perform sharp interface calculations based on Green's function methods and phase field simulations, supplemented by analytical investigations. For pure dilatational transformations we find a single growing finger with symmetry breaking at higher driving forces, whereas for shear transformations the emergence of twin structures can be favorable. We predict the steady state shapes and propagation velocities, which can be higher than in conventional dendritic growth.

Introduction

Transformations between different solid states of a material are essential for many technological and scientific applications, and understanding their kinetics is therefore not only interesting because of the scientific variety and beauty, but also important for tailoring of new materials with specific properties. Often, these processes are accompanied by elastic deformations, which can be e.g. due to density differences or structural changes, which provoke stresses at interfaces between adjacent phases. The influence on the thermodynamics of the transitions between different phases has been thoroughly discussed in the literature [START_REF] Khachaturyan | Theory of Structural Transformation in Solids[END_REF][START_REF] Roitburd | [END_REF], whereas the kinetics of these processes are still far less understood.

There are many important analogies between solidification and diffusion-limited solid-state transformation, and a recent review on similarities is given in [3]. These similarities include the growth of needles and plates similar to Ivantsov paraboloids, i.e. dendrites versus Widmanstätten structures, coupled growth of eutectics and eutectoids, and massive transformations. Differences are, on the other hand, the influence of fluid flow during solidification and elastic effects during solid-state transformations. E.g. Widmanstätten precipitates, which can grow as α needles from β brass, have little tendency for sidebranches and show approximately paraboloidal fronts near the tip region, whereas the side faces are locked into a low energy and mobility zone. They have been modeled using highly anisotropic interfacial free en-ergy [4]. Experiments clearly demonstrate that the growth of Widmanstätten ferrite responds remarkably to externally applied stresses [5], thus emphasizing the role of elastic effects. In real systems many other effects apart from elastic effects can come into play -like inhomogeneous compositions, the Mullins-Sekerka instability, crystal anisotropies, polycrystalline structures, etc. -and also quite successful attempts have been made to model this complex behavior (see for example in [6] and references therein).

From microscopic solvability theory of solidification it is known that the questions concerning existence, shape and growth velocity of steady state patterns in solidification crucially depend on selection mechanisms. In this spirit the influence of many different physical effects on solidification or melting processes has been studied throughout the years [7,8]. For example, the effect of isotropic surface tension was proven not to serve as a selection mechanism for a solution, whereas anisotropic surface tension leads to a unique solution. Here, we mention that the anisotropies of surface tension and interfacial mobility are generally more pronounced in solid-solid transformations, which often leads to faceted interfaces. Even more, high-temperature and bainitic transformation can proceed by the motion of transformation disconnections (defects with ledge/dislocation character), by the motion of growth ledges (defects with pure ledge character or with negligible dislocation content), or by a mixture of these [9,10]. Nevertheless, in the present work we discuss neither the role of anisotropy nor of these (discrete) defect structures but concentrate on the role of elastic effects; the main reason for this simplification is the striking result that steady state dendrite-like growth is even possible without anisotropy of surface tension due to elastic effects, in contrast to diffusion-limited solidification with only rough surfaces [11].

Phase field methods have emerged as a powerful tool to investigate also solidstate transformations; for a general overview see e.g. [12,13]. The microstructural development during precipitation of a coherent orthorhombic phase (O-phase) from an α 2 matrix (DO 19 ) in a Ti-Al-Nb system, which undergoes a hexagonal to orthorhombic transition, has been investigated through computer simulations using the phase-field approach. The orthorhombic phase has three orientation variants with equal energy and formation probability in an isotropic environment. The authors compared the appearing microstructure with experimental observations using transmission electron microscopy [14][15][16][17] and find similar patterns. In particular, they conclude that the application of a moderate external strain does change the essential features and morphology, but only the relative amounts of the phases, as long as the applied strain is small in comparison to the stress-free transformation strain. The underlying time-dependent Ginzburg-Landau equation, however, has not been investigated with respect to a valid sharp interface limit and lacks a detailed understanding of the growth kinetics in the spirit of an interfacial pattern formation process. Recently, Greenwood et al. used phase field methods to investigate two-dimensional dendritic structures in the solid state in the presence of competitive interactions between surface energy and elastic anisotropy [18]. They predict a transition from surface dominated growth along the [10] direction to an elastically driven growth along the [11] direction by changes in the elastic and surface anisotropy and supersaturation. However, there is a still substantial lack of both experimental and theoretical investigations concerning the growth kinetics of the emerging patterns.

The purpose of the present paper therefore is to shed light on the growth kinetics of diffusion-limited solid state transformations, and to highlight differences and similarities to solidification processes with emphasis on elastic effects and their influence on selection principles. The central question addressed here is related to the role of strain energy during solid-state phase transformation processes. We focus on gaining a fundamental understanding of the elastic influence on pattern formation processes, and we restrict our investigations therefore to simple and "clean" model systems. This includes the limitation to isotropic elasticity and surface energy, symmetrical models of heat or component diffusion, local equilibrium conditions at the propagating interfaces and the absence of defects and plastic effects. Even under these simplifying assumptions a very rich solution behavior emerges, providing new insights on the kinetics of diffusion-limited solid-state transformations. Apart from phase field methods, we also use sharp interface methods, that are based on Green's function techniques. More importantly, however, is that we supplement these numerical investigations by analytical calculations which allow to understand these processes and the relevant parameters on a deeper basis.

The formulation of such a model system requires to write down explicit expressions for the boundary conditions and suitable equations of motion at the propagating interfaces, which therefore have to be tracked during the entire dynamical process. In the end, it provides a detailed description of the microstructure evolution on the level of individual grains in a polycrystalline structure. Then, propagation of a front generates or consumes latent heat, which has to diffuse through the system, and we assume that this is the rate-limiting process. This is in contrast to isothermal solid-solid transformations that were discussed in [19], which are driven by local interface kinetics, therefore allowing fast front propagation also on the scale of the sound speed.

The paper is organized as follows: In Section 2 the basic problem of diffusion limited growth in the presence of lattice strain is posed in the spirit of a sharp interface description. In Section 3 different structural transitions are considered that lead, similar to the hexagonal to orthorhombic transition mentioned above, to effectively isotropic elastic media. Section 4 contains the modeling principles for growth in free space using Green's function methods, followed by phase field approaches for channel growth in Section 5. The results for both approaches and geometries are presented in Section 6.

Formulation of the Problem

Specifically, we consider the growth of an equilibrium solid phase β from another metastable solid phase α, as depicted in Fig. 1. At a temperature T eq , both stressfree phases are in bulk equilibrium with each other, whereas below this value the new β phase is thermodynamically favorable. Far ahead of the propagating front, the system is therefore exposed to a lower temperature T ∞ , and we define a dimensionless temperature w = C(T -T ∞ )/L, where L denotes the latent heat of the transition and C the heat capacity. Within the solid phases the evolution of the temperature field obeys a diffusion equation,

D∇ 2 w = ∂w ∂t , ( 1
)
where D is the thermal diffusivity, which we assume to be equal in both phases, as well as the heat capacity (symmetrical model).

The motion of the interfaces follows from energy conservation principles, which account here for the net heat flux and the generation or absorption of latent heat due to the local motion of the interface with normal velocity v n , where n is the interface normal and the superscripts denote the phases. Elastic effects enter into the problem via a modification of the local phase equilibrium condition at the moving interfaces. In contrast to previous investigations [20,21], we now assume the phases to be coherently connected, i.e. we neglect the appearance of any slips, detachments or defect formation at the interfaces, which is reasonable for small misfits. It is important to mention here, that with the assumption of coherency we formulate a completely different model being capable to describe qualitatively different types of elastic effects. For example the wellknown effect of elastic hysteresis [START_REF] Khachaturyan | Theory of Structural Transformation in Solids[END_REF][START_REF] Roitburd | [END_REF], which was proven to disappear without coherency [20], is now included in the model. We point out that the validity of local equilibrium conditions without kinetic corrections at interfaces is probably less obvious for solid-state phase transformations as for solidification in the regime of low undercoolings. Often anisotropy effects and faceting of interfaces can change this behavior, but nevertheless we use the local equilibrium approximation here, together with the restriction to isotropic systems. Also, for diffusion of solute instead of heat, partitioning is usually reached only for fast diffusing elements as carbon in steel, whereas other substitutional elements remain in paraequilibrium at the interface.

v n =Dn • (∇w [β] int -∇w [α] int ), (2)
In a Lagrangian description of linear elasticity, the displacement field is denoted by u, and the coherency condition reads, u [α] = u [β] at the interface. The strain field is then given by the symmetrical spatial derivative, ik = (∂u i /∂x k + ∂u k /∂x i ) /2. This implies that tangential and shear strains τ τ , sτ , ss are continuous at the interface with τ, s being the two tangential directions. Force balance demands the continuity of the normal and shear stresses, σ nn , σ nτ , σ ns , where n is the normal direction. Here, the stresses are defined as the derivative of the free energy with respect to the strains, σ ik = ∂F/∂ ik .

The interface temperature for a planar, stress free interface is given by w = ∆, with ∆ = C(T eq -T ∞ )/L. Capillary corrections induce a curvature dependent term, and for simplicity we discuss only the case of isotropic surface tension γ. Finally, elastic deformations lead to an additional change of the interface temperature. Altogether, this leads to the local equilibrium boundary condition [11]

w int = ∆ -d 0 κ + T eq C L 2 δ Fel , (3)
with d 0 = γT eq C/L 2 being the capillary length and κ the local curvature of the interface, counted positive for a convex phase β (see Fig. 1); Notice that due to the coherency constraint [22,23], the elastic shift of the equilibrium temperature is proportional to the difference of a new potential δ Fel = F

[α] el - F [β]
el , which is related to the elastic free energy F el via the Legendre transformation Fel = F elσ nn nn -2σ nτ nτ -2σ ns ns . From now on, we restrict our considerations to linear isotropic elasticity. Taking the relaxed state of the α phase as reference, the elastic contribution to the free energy there reads

F [α] el = E 2(1 + ν) ν 1 -2ν 2 ii + 2 ik , (4)
where E and ν are Young's elastic modulus and Poisson ratio, respectively. We use here the sum convention for repeated indices for the sake of brevity. In contrast, the elastic contribution to the free energy density of the new phase β reads

F [β] el = E 2(1 + ν) ν 1 -2ν ii -0 ii 2 + ik -0 ik 2 , (5)
which has an other state of zero elastic energy due to the lattice strain 0 ik assigned to the phase transformation (see for example (author?) [START_REF] Khachaturyan | Theory of Structural Transformation in Solids[END_REF][START_REF] Roitburd | [END_REF]). In the following we will call 0 ik the eigenstrain of the transformation. Particular cases of eigenstrains and their interpretation will be discussed in the next section.

The mechanical equilibrium and coherency conditions provide expressions for the discontinuous jumps of the strains at the interface

[β] nn -[α] nn = 0 nn + ν 1 -ν ( 0 ss + 0 τ τ ), ( 6
) [β] nτ -[α] nτ = 0 nτ , [β] ns -[α] ns = 0 ns , (7)
which are only related to the imposed eigenstrain. Defining an eigenstress tensor via Hooke's law

σ 0 ik = E 1 + ν 0 ik + ν 1 -2ν δ ik 0 ll , (8)
we obtain after a few straightforward algebraic manipulations for the elastic contribution to the local equilibrium condition Eq. ( 3)

δ Fel =σ 0 ik [α] ik - E 2(1 -ν 2 ) ( 0 τ τ ) 2 + ( 0 ss ) 2 + 2ν 0 τ τ 0 ss + 2(1 -ν)( 0 τ s ) 2 . (9)
Apparently, this expression depends on the strain state of the α phase at the interface.

As already implicitly assumed above, the elastic degrees of freedom relax fast on diffusive timescales, and therefore the application of static elasticity is legitimate. Therefore, Newton's second law becomes

∂σ ik ∂x k = 0. ( 10
)
The above set of equations now has to be supplemented by boundary conditions at the external system boundaries, and we will discuss different scenarios below. They complete the self-sonsistent description of the moving boundary problem, which require the simultaneous solution of the elastic and diffusion equations, coupled via the boundary conditions at the propagating interfaces, leading to a complicated nonlinear and nonlocal problem.

We note that we will later also discuss more complicated situations with more than one growing phase. Then, also similar boundary conditions and equations have to be set up for the interfaces between them. Typically, then also triple junctions will appear, where the contact angles are given by Young's law.

Models of structural transitions

The simplest case of a characteristic lattice strain is to assume the bond length of the new phase β to be uniformly longer or shorter in all directions in compar- ison to the reference phase, i.e. 0 ik = d ik = δ ik , which implies a density change. Notice, that this kind of transformations, to which we will refer to as dilatational eigenstrain, only changes the volume of an elementary cell, and not its relative dimensions or its scaled shape.

The opposite case are transformations, where the volume of the unit cell is conserved and only the shape is changed. These shear transformations are characterized by a traceless eigenstrain tensor. A particular transition involving shear strain can occur e.g. in hexagonal crystals, where the transition leads to a lowering of the symmetry from C 6 to C 2 . This is the case for example in hexagonal-orthorombic ferroelastics [24], or in the above mentioned Ti-Al-Nb systems [14][15][16][17]. Let the principal axis C 6 be oriented in z direction. For simplicity we neglect all other possible strains with higher (axial) symmetry. By proper choice of the crystal orientation around the main axis in the initial phase, we obtain the new phase in three possible states in the xy plane due to the original hexagonal symmetry (See Fig. 2). Then the non-vanishing components of the eigenstrain tensor 0 ik = s ik are s yy = -s xx = cos 2ϑ and s xy = sin 2ϑ, where the orientation angle ϑ is either ϑ = 0 or ϑ = ±2π/3.

We also consider a more generic mixed mode case, where the self-strain of the transformation is a linear combination of the above two different transformations. We will discuss the following two mixing cases

0 ik = ± ik (η) = η d ik ± (1 -η) s ik ( 11
)
and refer to them as positive or negative mixing.

It is important to mention that switching from positive to negative mixing corresponds to a rotation of the eigenstrain tensor by an angle π/2. Since our model is isotropic in all respects except the choice of the eigenstrain tensor, the only predefined direction is given by the growth velocity. Therefore, we can interpret the negative mixing case to represent physically the same system as the positive mixing case but with a growth direction perpendicular to it.

Free growth

As a first application, we discuss the growth of the stable β phase from the metastable α phase in an infinitely large system, which is exposed to the temperature T ∞ far away from the interface. For simplicity, we discuss an effectively two-dimensional infinite system, assuming translational invariance in z direction. In particular, we assume a plane strain situation from point of view of elasticity. Also, we choose the elastic boundary conditions such that all stresses decay far away from the interface.

Seeking for steady state solutions in free space, we use Green's function methods for solving both the elastic and the diffusional problem, to finally obtain selfconsistently the shape and the velocity of the growing front. While the derivation of the integral representation of the thermal field is well-known, see e.g. [25], we present the calculation of δ Fel in more detail [11]. To this end, we express the elastic problem using the eigenstrain-independent stress tensor σik ; for the α phase it is the same as the nomimal stress, σ

[α] ik = σ [α] ik , whereas in the β phase it is σ[β] ik = σ [β]
ik + σ 0 ik . In terms of the new stress tensor σik , the mechanical equilibrium equation ∂ σik /∂x k = f i introduces a force density which is localized on the interface and vanishes in the bulk,

f i = σ[β] in - σ[α] in = σ 0 in . (12)
This means that the problem of two coherently connected solids with different reference state is equivalent to a monolithic material with a distribution of point forces along the interface. We can therefore use an integral representation of the displacement field,

u i = ds G ik (r -r )f k (r ), (13)
where the integration is performed along the interface, and the free-space Green's tensor G ik (rr ) of the elastic problem is given by [START_REF] Landau | Theory of Elasticity[END_REF] G

ik (r) = 1 + ν 4π(1 -ν)E x i x k r 2 -(3 -4ν)δ ik ln r . ( 14
)
This function is the displacement response at r to a point force located in the origin. As a result, the mismatch provokes long-ranged elastic deformations, and therefore induces a nonlocal character to the problem. The above expressions can then be used to calculate the elastic influence on the interface temperature w int . Far behind the tip the interface profile becomes parabolic and is described asymptotically by the Ivantsov parabola [START_REF] Ivantsov | [END_REF]. In this region, the elastic strains have decayed, and only the constant local contribution to δ Fel (second term in Eq. 9) remains. Consequently, the effective driving force ∆ consists not only of the thermal undercooling ∆, which is controlled by the applied temperature T ∞ , but also a contribution ∆ el from the elasticity,

∆ el = CT eq E[( 0 yy ) 2 + ( 0 zz ) 2 + 2ν 0 yy 0 zz ] 2(1 -ν 2 )L 2 . ( 15
)
We can consider ∆ el as a parameter quantifying the strength of the elastic effects. The elastic effects therefore induce a shift of the equilibrium temperature, and the effective undercooling ∆ is lower than the actually applied thermal driving force,

∆ = ∆ -∆ el , (16)
which reflects the elastic hysteresis [START_REF] Khachaturyan | Theory of Structural Transformation in Solids[END_REF][START_REF] Roitburd | [END_REF]. By elimination of the thermal field in the spirit of [25] and eliminating the elastic fields, with the use of the integral representation for the displacements Eq. ( 13), we can derive the steady state equation for the shape of the solid-solid interface in a closed dimensionless representation, which reads in the co-moving frame of reference [11],

∆ - d 0 κ R + T eq Cδ Fel L 2 = p π ∞ -∞ dx exp[-p(y(x) -y(x ))]K 0 (pµ(x, x )). ( 17
)
Here, p = vR/2D denotes the Peclet number, which is proportional to the product of the steady state velocity v (see Fig. 1) and the radius R of the tip-curvature of the asymptotic Ivantsov parabola.

(x, x ) = [(x -x ) 2 + (y(x) -y(x )) 2 ] 1/2
. Note that all lengths in the above expression were rescaled by the Ivantsov parabola radius R.

Since heat release as well as elasticity depend nonlocally on the interface shape and vice versa, the central equations ( 1) -( 3), together with the expression for the elastic contribution (9), result in the complicated integro-differential equation (17) for the interface shape. It differs from conventional dendritic growth by the appearance of the elastic term. To better understand its influence, we briefly review the solvability conditions for classical dendritic growth [7,8], for situations with pure dilatational and shear strain [11], and cases with a triple junction [28] for the growth along a grain boundary: In classical dendritic growth, any solution for isotropic surface tension with nonzero solvability parameter σ = d 0 /pR must exhibit a finite negative tip kink angle φ. For this reason, dendritic growth is not possible for isotropic surface energy, since no smooth dendrite profiles exist. A dilatational eigenstrain simply shifts the equilibrium interface temperature, and therefore suffers from the same missing selection as the dendrite problem without surface tension anisotropy. The presence of a shear eigenstrain, which will also be discussed in more detail below, leads to the propagation of a bicrystal. In this case, the elastic effects operate as new selection mechanism. Independent of this elastically induced selection mechanisms, the presence of a triple junction alone also works as a selection mechanism. Here, the tip is not smooth with a dihedral angle that is determined by the grain boundary and interfacial free energies.

Concerning the numerical solution of Eq. ( 17) we point out, although it represents the whole problem now also including elastic effects, it is still formally similar to closed representations known from conventional dendritic growth. Therefore, we can solve the corresponding eigenvalue problem numerically for the steady state shape and the stability parameter σ = σ * exactly in the spirit of Ref. [29].

The results are presented in Section 6.

Channel Growth

The growth of the β phase in a channel of finite width W differs significantly from the free growth scenario discussed above. The specific behavior depends of course on the boundary conditions, and we assume fixed grip conditions on the long sides of the channel, i.e. the displacement is constant. Again, we discuss an effectively two-dimensional system by the assumption of translational invariance in z direction (see Fig. 3 for the setup). Furthermore, we assume that the α phase far in front of the growing tip is fully relaxed, which is in agreement with the boundary conditions u x ≡ 0 and u y ≡ 0 along the boundaries of the whole, infinitely long strip. From a thermal point of view, we assume thermal insulation along these boundaries, so no temperature flux through the surfaces, ∂T /∂x = 0 on the outer boundaries. Far ahead of the tip, the temperature is given by the undercooling temperature T ∞ . Due to the phase transition, the release of the latent heat leads to a heating of the system, therefore the temperature far behind the tip T -∞ is higher than in front. Under the given circumstances it is quite obvious that in this tail region a homogeneous state is reached. In particular, since the interfaces do not move there, the temperature becomes indeed constant.

For simplicity, we discuss here only the specific case of a purely diagonal eigenstrain, i.e. ij = 0 for i = j, to demonstrate the analytical procedure. In contrast, in the case of an eigenstrain tensor with nonzero off-diagonal elements, bicrystal patterns emerge, and we will briefly discuss this situation below; it is quite obvious that the analysis can be performed in an analogous way.

Due to the confinement of the system, phase coexistence is possible in a finite range of undercoolings, and we can calculate the asymptotic volume fraction λ of the new phase (see Fig. 3). Notice, that in contrast to the isothermal situation discussed in [19], the temperature in the tail region is not a control parameter but has to be found self-consistently, due to the thermally insulating boundaries. In terms of the asymptotic phase fraction λ, the width of the β phase in the tail region is written as λW . Mechanical compatibility requires (1λ) [α] xx + λ

[β] xx = 0, (18)
since the average lateral strain must vanish due to the boundary conditions for the displacement. Together with the mechanical equilibrium condition of force balance, we can then calculate the elastic fields in the tail region, and consequently the elastic energy. Then, we calculate the energy excess δF, due to the appearance of a finite fraction of the new phase λ. It is defined as the difference between the energy for λ = 0 without any phase β and the energy with an arbitrary but finite fraction,

δF = -W λF [β] el + (1 -λ)F [α] el -2γ + W λ 0 L(T eq -T (λ )) T eq dλ . (19)
Initially, there is only the elastically relaxed but metastable phase α at the constant temperature T ∞ . As soon as we have a finite phase fraction λ, the energy increases by elastic and capillary contributions. Finally, the integral appears due to thermal insulation, since an increase of the amount of the phase β by phase transformation processes generates latent heat, and therefore causes an increase of the temperature.

Using the heat conservation condition we find the relation between this temperature and the phase fraction to be

T (λ) = T ∞ + λ L C . ( 20
)
Far behind the tip the fraction λ of the new phase is constant, which yields a homogeneous elastic state there. While 18), as well as the interface jump conditions Eq. ( 6) and Eq. (7). Consequently, for the non-vanishing components of the stress tensor we obtain σ

[α/β] xx = -σ 0 xx λ, σ [α/β] xy = 0, σ [β] yy = ν(1 -λ)σ 0 xx /(1 -ν) -σ 0 yy and σ [β] zz = ν(1 -λ)σ 0 xx /(1 -ν) -σ 0 zz .
Inserting this into the elastic contribution of the free energy excess, we get

δF el (λ) = -W λ 1 2 σ [β] ik [β] ik -0 ik + 1 2 (1 -λ)σ [α] ik [α] ik = 1 2 λW σ xx 0 xx + σ [β] yy 0 yy + σ [β] zz 0 zz (21)
Finally, using Eq. ( 20) to calculate the thermal energy, the energy excess as a function of the asymptotic phase fraction λ becomes where the parameter ∆ el , as already introduced above in Eq. ( 15), defines the strength of elastic effects and

δF(λ) = W L 2 CT eq (∆ -∆ el )λ - 1 2 (1 + a el ∆ el ) λ 2 - 2d 0 W , (22)
a el = 2 (1 -2ν) (1 -ν) 0 xx + ν 0 yy + 0 zz 2 ( 0 yy ) 2 + ( 0 zz ) 2 + 2ν 0 yy 0 zz , (23)
is a parameter introduced for brevity of presentation, which basically depends on the type of eigenstrain. Note that in general this parameter a el depends also on the off-diagonal component of the eigenstrain tensor 0 xy , but we discus here only the case of vanishing total shear eigenstrain.

The maximum of the energy excess, ∂∆F (λ)/∂λ = 0, provides the equilibrium fraction, which expresses the thermodynamic balance between elastic deformation and free energy release due to the phase transition,

λ = 1 1 + a el ∆ el (∆ -∆ el ) . (24)
If the fraction λ is in the range 0 < λ < 1, we obtain phase coexistence in the asymptotic regime far behind the tip. Then growth of the β phase demands that the energy excess (22) for the equilibrium fraction (24) has to be positive, since the total free energy has to decay in this dissipative process and far behind the tip thermodynamic equilibrium has to be reached. Hence we get the requirement

δF(λ) = W L 2 CT eq 1 2 (1 + a el ∆ el ) λ 2 - 2d 0 W > 0. (25)
This is equivalent to the condition λ > λ crit , where λ crit is given by

λ 2 crit = 4d 0 W 1 1 + a el ∆ el . ( 26
)
Since the fraction of new phase is proportional to the driving force, the growth condition Eq. ( 25) demands that the driving force obeys

(∆ -∆ el ) 2 > W 4d 0 (1 + a el ∆ el ), (27)
in order to make the growth in the channel energetically possible. We note that for the shear transformations discussed above, the most favorable process is the growth of a bicrystal, where two twinned phases β and β with ϑ = ±2π/3 grow together, because this is energetically most favorable [19]. Then also the β/β grain boundary energy appears in the energy balance analogous to Eq. (19). For the bicrystal configuration the analysis can be performed in a similar way, and even the parameter a el remains in the presented form given in Eq. ( 23), since the mean shear eigenstrain of the composite system still vanishes.

Considering the channel geometry, we use a phase field formulation to solve the full free boundary problem numerically. This method is very versatile and allows to study even more complicated scenarios, e.g. the full dynamical evolution and not only the steady state regime, three-dimensionsal situations with full anisotropy of elasticity and surface tension or simulations beyond the symmetrical model, where the two phases have different diffusion constants. Nevertheless, for brevity of presentation, we confine our investigations to the simple, isotropic, two-dimensional and thermally symmetric case that we also used for the Green's function approach in the free space geometry and the above analytical calculations. The phase field method introduces an additional numerical lengthscale ξ, the width of the transition region between the different solid phases. This means that an order parameter, which discriminates between the different phases, has a profile φ = [1 + tanh(n/ξ)]/2 in normal direction through an interface. Although the system size and the tip scale are chosen to be significantly smaller than the interface thickness and the results do not seem to depend significantly on these parameters, they may be improved by using thin interface techniques, as explained in [30]. In the sharp interface limit the description recovers the desired sharp interface equations that were proposed above.

Hence, we introduce the phase field φ having a value φ = 1 in the metastable initial phase, and a value φ = 0 in the new phase. We start from a free energy functional

F [φ, w, u i ] = V (f s + f dw + f ) dV, ( 28
)
where f s (∇φ) = 3γξ(∇φ) 2 /2 is the gradient energy density and f dw (φ) = 6γφ 2 (1φ) 2 /ξ is the double well potential, guaranteeing that the free energy functional has two local minima at φ = 0 and φ = 1 corresponding to the two distinct phases of the system. The form of the double well potential and the gradient energy density are chosen such that the phase field parameter ξ defines the interface width and the parameter γ corresponds to the interface energy of the sharp interface description [31]. The free energy density f (φ, w, u i ) of our phase field description is chosen to be an interpolation between the two distinct free energies of the model,

f (φ, u i , w) =F [α] el h(φ) + F [β] el (1 -h(φ)) + L T -T eq T eq (1 -h(φ)), (29)
where we choose the interpolation function h to be h(φ) = φ 2 (3 -2φ). It is the simplest polynom satisfying the necessary interpolation conditions h(0) = 0 and h(1) = 1, and having a vanishing slope at φ = 0 and φ = 1, in order not to shift the bulk states. The evolution equation of the phase field is given by the variational expression

∂φ ∂t = - M 3γξ δF δφ ui,w . (30)
For large values of the mobility M we recover the desired case of diffusion limited growth. The elastodynamic equations read

ρü i = - δF δu i φ,w = ∂σ ik ∂x k , (31)
which recover static elasticity for slowly propagating interfaces in comparison to the sound speed v s ∼ (E/ρ) 1/2 [19], with ρ being the mass density.

For the temperature field we have the usual thermal diffusion equation with the motion of the phase field or interface being a source of latent heat, where the prime denotes the derivative with respect to φ.

∂w ∂t = D∇ 2 w + h (φ) ∂φ ∂t , (32)
The phase field model presented here is very similar to the model in [32]. However, there only the very special case of melting or crystallization of a confined sphere with elastic effects was discussed.

For the solution of the coupled partial differential equations ( 30)-( 32) we use explicit and parallelized finite difference methods, which are carefully explained in [33]. All simulations have been done on the parallel computer JUGENE operated at the Research Center Jülich. For the simulations we typically use parameters W/d 0 = 100 and W/ξ = 80.

Results

Apart from the theoretical investigations above, we also performed numerical simulations for both the free growth in an infinite system and in a channel or narrow strip. As already mentioned before, we observe the growth of a single new phase if the eigenstrain tensor is purely diagonal, and twin structures for shear transformations with ϑ = ±2π/3. Also, we consider mixed mode situations which are between these two extreme cases.

By construction, the free growth can be captured only by the Green's function method, whereas the channel geometry can be simulated only by the phase field method. The latter is not yet fully quantitative, as the results depend weakly on the interface thickness, but nevertheless we see a very reasonable agreement between the two methods and geometries for rather wide channels.

Single crystal

We first consider the situation of free growth. Here, the simplest case is that of a dilatational eigenstrain where all diagonal elements have the same value, i.e. 0 ik = δ ik ; this case is obtained from Eq. ( 11) by setting η = 1. As already mentioned before, selection of a steady solution in the absence of elastic effects is not possible for isotropic surface tension. It is quite remarkable that also the inclusion of elastic effects due to the dilatational eigenstrain of the transition does not change this situation; the reason is that the elastic effects can be mapped to a change of the undercooling only, and therefore we recover again the missing selection of conventional free dendritic growth [11]. It turns out that also for the volume-preserving shear strain transitions with 0 yy = -0 xx = , where all other components 0 ik vanish, i.e. η = 0 and ϑ = 0, steady state solution do not exist either.

Under these circumstances it is very suprising, that for values of the mixing parameter in between zero and one 0 < η < 1, ϑ = 0 we do find steady state solutions for free growth and negative mixing. This is shown in Fig. 4a, where the selected "stability parameter" σ = d 0 /pR = σ * , which is a dimensionless measure for the growth velocity, is plotted as function of the mixing parameter η. In general, the eigenvalue σ * depends on η, ∆ el and p, but in the regime of small Peclet numbers the elastic driving force and p appear only in the combination ∆ el /p. Fig. 4b shows the dependence of the selected stability parameter on ∆ el /p for two different values of η. For the other type of mixing + ik corresponding to the unprivileged direction of growth, we do not find steady state solutions in free space. We note that from the selected stability parameter σ * one can extract the steady state growth velocity v via the relation

v = 2D d 0 σ * p 2 , (33)
which is obtained from of the above mentioned definition of the stability parameter for fixed Peclet number. In contrast, for the growth in a channel, steady state solutions do also exist for the case of a dilatational eigenstrain, η = 1 (see Fig. 5a). As expressed through Eqs. ( 25) and ( 26), the phase β can not grow for phase fractions λ below the threshold λ crit . Beyond this value the growth of a symmetrical finger (dendrite) is the most favorable solution. We note that in situations where free growth solutions do not exist, the still existing solutions in the channel must have a scale that is determined by the size W of the channel. For higher driving forces an asymmetric finger grows with higher velocities, and is therefore more likely to be observed; typical finger shapes are shown in Fig. 3.

For the situation of η = 0.5, where steady state solutions exist also in the free space, we performed a quantitative comparison of the results of the two methods in Fig. 5b. For this comparison we directly relate the fraction λ to the driving force ∆ for free growth via Eq. ( 24) and Eq. ( 16). Then, in the limit of small Peclet numbers, the latter is related to ∆ according to p ≈ ∆2 /π. Finally, via the calculated function σ * (∆ el /p) and Eq. ( 33) we obtain the growth velocity as a function of λ, which is then compared to the phase field results. The aim of this comparison is to show that the growth velocities obtained from the two different kinds of growth, free growth and channel growth, agree quantitatively in there order of magnitude and therefore are selected by the same underlying physical principle.

The results for the growth velocity agree well in a range of intermediate fractions of the new phase β. In particular, we confirmed that the steady state solutions which were obtained by the Green's function method are indeed stable, as the dynamical phase field method tracks only these stable branches. For low phase fractions, growth is not possible below the threshold λ crit , which tends to zero for infinitely wide channels. For very high driving forces, the growing finger is influenced by the geometrical confinement of the channel width and also affected by kinetic corrections to the sharp interface equations of motion, which therefore lead to a deviation from the free growth result which was obtained by the Green's function method.

It is interesting to compare the situation to a case without elasticity, where the problem reduces to classical growth in a channel [34][35][36][37][38][39]. It is known that there even in wide channels symmetric solutions do not exist below a fraction λ = 1/2, whereas we found that the eigenstrain leads to the existence of solutions at already significantly lower driving forces λ, see Fig. 5a.

Bicrystal growth

In contrast to the cases studied above, the presence of off-diagonal eigenstrain entries lead to the preferential growth of bicrystal patterns. This behavior can easily be understood intuitively, since the eigenstrain in the β phase favors the presence of off-diagonal strain. Then the boundary conditions of fixed displacements require that the surrounding α phase is sheared in the opposite direction. This implies a finite elastic energy which is stored far behind the tip. If, on the other hand, a second finger β of the same size with an off-diagonal eigenstrain with opposite sign appears, the elastic energy contribution that is associated with the off-diagonal elements of the strain tensor vanishes. Consequently, the total elastic energy in the tail region is reduced, and the growth of the bicrystal therefore thermodynamically favored.

An additional energy contribution emerges from this bicrystal configuration, since the boundary between β and β generates a grain boundary energy γ gb . We consider in particular the case of vanishingly small twin boundary energy, γ gb γ. This affects also the geometrical situation at the tip, which is now a triple junction, and the contact angle φ becomes zero, so the overall finger has a smooth contour. Typical shapes of such bicrystal pattern are shown in Fig. 6, which are here related to the growth in a channel, but they look similar for free growth. We note that we simulated only half of the channel using the symmetry of the growing structures [19].

Considering the free space geometry, the first important observation is that in contrast to the single crystal scenario, steady state solution do exist also for pure shear, η = 0, which is already discussed in [11]. New is here that we found the same morphology also in the channel geometry, showing surprisingly good agreement in the comparison of growth velocities, done in the same way as was done for the single crystal described above (see Fig. 7). Again, we would like to point out that the stability of the bicrystal solution obtained from the Green's function technique is confirmed by the phase field simulations.

We also did investigations of the bicrystal free growth considering mixed eigenstrains, and the results are summarized in Fig. 8 and9. In the first graph, the stability parameter σ * , which is again a measure for the dimensionless velocity, is shown as function of the admixture parameter η for fixed ratio ∆ el /p. As before, the elastic driving force and the p appear only the combination ∆ el /p for small Peclet numbers.

As already mentioned above it is quite interesting that the appearance of an additional dilatational strain leads to an effective anisotropy of the system: Whereas in the case of a positive admixture + the velocity grows with increasing η, the opposite happens for negative admixture -. This suggests that in free growth the bicrystal structure will eventually grow in the direction with the higher propagation velocity. For the fast branch, we did not find solutions beyond the point η ≈ 0.25.

In Fig. 10 we show the comparison of the two geometries for the bicrystal case. Again, the phase field results confirm the stability of the solution branch for low driving forces. In this region, we obtain a convincing quantitative agreement between the growth velocities obtained for the free growth scenario using the Green's function method, and the channel geometry with the phase field method. In the left panel the favored direction is shown (positive admixture) for η = 0.1. In this case we again observe dendritic solutions for small driving forces λ (see also Fig. 6a). For the chosen parameters in the channel geometry, all solutions on the upper branch are dendritic. On the lower branch, however, the steady state shapes are double-fingers, which means here the curvature directly at the point of symmetry has changed its sign and we obtain two tips (see Fig. 6b). The velocity of these branched solutions is here always smaller than for the single-tip structures.

For the other direction, i.e. negative admixture, and a mixing parameter η = 0.5, we show the comparison of growth velocities in the right panel of Fig. 10. Here, we obtain a convincing agreement between the different methods and associated geometries at least for small driving forces, indicating dynamical stability also for this type of solution. Obviously, the growth velocities are here substantially lower than for the previous case, and for fully free growth an initial bicrystal seed would therefore preferentially grow in the "fast" direction.

Summary and Conclusion

In summary, we investigated the diffusion-limited kinetics of solid-solid transformations in the presence of lattice strain. We studied in particular the free growth with Green's function methods and the growth in a channel by phase field techniques, exhibiting a comparable behavior. It turns out that the elastic effects have a strong influence on the selection of the arising patterns and the corresponding steady state velocities, which differ crucially from conventional dendritic growth.

In particular, we find that the selected stability parameter and the corresponding growth velocity can be much higher compared to the classical case of dendritic growth, where the solutions are selected by the tiny effect of the anisotropy of surface tension. Different structural transitions, which are described by dilatational and shear eigenstrains, lead to a rich behavior, and we find very different structures already for the simple cases that were discussed here.

Here, we focused on growth situations that are controlled by thermal diffusion, but it is known that the also interesting scenario of chemical diffusion can be mapped to this case.

In comparison to the Green's function method, which is computationally very efficient, the phase field technique has the advantage of higher flexibility. In particular, it can easily be extended to three-dimensional cases, it can be used also for non-steady state situations or anisotropic or non-symmetrical situations with different material parameters in the phases.

To our mind, this work offers a new perspective on diffusion-limited solid-state transformations, but nevertheless it leaves several open questions: What is the role of anisotropy of surface tension, kinetics and elastic effects? How do they influence the general predictions made in the present paper? What happens in particular for materials with rough and facetted interfaces, where the growth may be determined by other mechanisms and selection principles? To which extend is the present theory still applicable? What happens for large deformations, where coherency at the interfaces may be lost and defect emission play a major role, or where irreversible plastic deformation may become substantial? It is tempting to speculate that these additional dissipative processes effectively slow done the kinetics, since they lead to additional strain relaxation. From point of view of modeling, further steps towards a fully quantitative phase field modeling, probably in the direction of [18], seem promising to shed light on these questions. We hope that this work will stimulate further experimental and theoretical investigations, but also re-investigations of earlier experiments in view of these unresolved questions.

This work was supported by the DFG grant SPP 1296 and the German-Israeli Foundation.

F

. Schematic figure of the steady state growth, with constant velocity v, of a thermodynamically favored phase β into a metastable phase α, which is initially at temperature T∞. The motion of the interface is locally described by the interface-normal velocity vn, which is given by the local heat conservation Eq. ( 2). The dashed lines at the two side borders indicate that we will discuss two different types of outer boundary conditions: Free growth, where the strains decay for x, y → ∞, and channel growth, with fixed displacements and thermal insulation at the walls. In the growth direction y the system is assumed to be infinite in both cases. with ∆ el = 0.4 and an undercooling, which is chosen such that according to Eq. ( 24) it would correspond to an asymptotic phase fraction λ = 0.65. While a) shows the shape of the symmetrical finger, b) shows the shape of the sightly faster parity broken finger (see Fig. 5a).

In panel c) we have plotted the shape in the case of negative mixing eigenstrain, η = 0.5, which in that case corresponds to growth in the privileged direction. The parameters are ∆ el = 0.05 and λ = 0.6 (for the resulting velocity see Fig. 5b).

Figure 4. The left panel shows the stability parameter σ * (see Eq. ( 33)) as function of the mixing parameter η for the privileged growth direction with an eigenstrain tensor - ik when ∆ el /p = 1. On the right panel we show the stability parameter σ * (∆ el /p) for two different values of mixing η.

Figure 5. Channel growth with diagonal eigenstrain. The left panel shows the dependence of the growth velocity on the driving force λ, where we choose a dilatational eigenstrain (η = 1, ∆ el = 0.4). For small driving forces we obtain velocities and shapes corresponding to a symmetric solution (see Fig. 3a). For higher driving forces a first order transition towards asymmetrical solutions takes place (see parity broken finger morphology in Fig. 3b). In the right panel the comparison of growth velocities between phase field and the free space calculations is shown, for the case of negative mixing corresponding to the privileged growth direction (η = 0.5, ϑ = 0, ∆ el = 0.05). The agreement is reasonably good, especially for driving forces λ slightly above the critical phase fraction λ crit of zero growth velocity. Figure 6. Bicrystal patterns growing in a finite channel. a) shows the single tip solution and b) the twinned finger morphology, which corresponds to the most favorable solution branch for higher driving forces. The positive mixing chosen for a) and b) leads to growth along the privileged direction and the parameters are η = 0.1, ∆ el = 0.1, λ = 0.6. The panel c) shows the shape of the slow growth direction for η = 0.5, ∆ el = 0.2, λ = 0.6. Figure 7. Bicrystal with pure shear eigenstrain, η = 0: Comparison of growth velocities obtained from the boundary integral technique for the free space geometry [11] with the velocities from phase field simulations for the channel geometry (∆ el = 0.3). Figure 8. Plot of the stability parameter σ * (η) as function of the admixture parameter η, for the two different types of mixing and ∆ el /p = 1.0. The solid line corresponds to the privileged direction of growth, i.e. positive mixing, while the dotted line shows the dependence for the unprivileged negative mixing configuration.

Figure 9. The left panel shows the stability parameter σ * as function of the control parameter ∆ el /p, for positive mixing, whereas the in the right panel σ * (∆ el /p) is plotted for negative mixing. Additionally, both graphs contain the σ * (∆ el /p) dependence for the pure shear configuration η = 0, as a reference to compare both types of mixings, which was taken from [11].

6

 6 Phase transformations, microstructure evolution, and deformation

  [α/β] yy vanishes asymptotically, the strain components [α/β] xx and [α/β] xy are determined by the mechanical compatibility Eq. (

Figure 2 .

 2 Figure 2. Visualisation of the hexagonal to orthorhombic transformation. From the symmetry its obvious, that the orientation angle ϑ of the new orthorhombic phase can have only the three different values 0, ±2π/3.

Figure 3 .

 3 Figure 3. Phase field shapes of single β phases growing with constant velocity v into the metastable α phase in a finite channel of width W . The shapes in a) and b) correspond to the dilatational case (η = 1), with ∆ el = 0.4 and an undercooling, which is chosen such that according to Eq. (24) it would correspond to an asymptotic phase fraction λ = 0.65. While a) shows the shape of the symmetrical finger, b) shows the shape of the sightly faster parity broken finger (see Fig.5a). In panel c) we have plotted the shape in the case of negative mixing eigenstrain, η = 0.5, which in that case corresponds to growth in the privileged direction. The parameters are ∆ el = 0.05 and λ = 0.6 (for the resulting velocity see Fig.5b).

Figure 10 .

 10 Figure 10. Plot of the dimensionless growth velocity as function of the "driving force" λ for bicrystal growth with mixed eigenstrain: comparison of free growth results obtained by Green's function methods to phase field simulations for the growth in a channel. Left panel: positive admixture with η = 0.1. Right panel: negative admixture, η = 0.5.

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters


Francis Corson

Hervé Henry

Mohktar Adda-Bedia

A model for hierarchical patterns under mechanical stresses

Keywords: numerical modelling, pattern dynamics user supplied): mechanical instabilities, phase field

We present a model for mechanically-induced pattern formation in growing biological tissues and discuss its application to the development of leaf venation networks. Drawing an analogy with phase transitions in solids, we use a phase field method to describe the transition between two states of the tissue, e.g. the differentiation of leaf veins, and consider a layered system where mechanical stresses are generated by differential growth. We present analytical and numerical results for one-dimensional systems, showing that a combination of growth and irreversibility gives rise to hierarchical patterns. Two-dimensional simulations suggest that such a mechanism could account for the hierarchical, reticulate structure of leaf venation networks, yet point to the need for a more detailed treatment of the coupling between growth and mechanical stresses.

A broad range of patterns in nature are induced by mechanical forces. While the most well understood examples are found in physical phenomena, e.g. in solids, fracture [START_REF] Freund | Dynamic Fracture Mechanics[END_REF], wrinkling [2], delamination [3], or localized deformation [4], it is increasingly becoming appreciated that mechanical forces are also essential in the morphogenesis of living organisms. In the simplest case, biological tissues can undergo mechanical instabilities as a result of inhomogeneous growth [5,6]. More broadly, the importance of the mechanical and adhesion properties of cells in morphogenesis [7] and the ability of mechanical forces to act as a signal in growing tissues and organisms [8,9,10] suggest that the interplay between mechanical forces and biochemical factors has a central role in development.

An illustration of these ideas can be found in recent investigations of the structure of leaf venation networks [11]. The standard model for vein formation involves the transport of the plant hormone auxin [12]: auxin circulates through the leaf, with the leaf base acting as a sink; as the leaf develops, auxin flow becomes progressively concentrated along preferential paths, which mark the location of future veins. This model has received ample experimental support (see e.g. [13]), yet it fails to explain certain aspects of leaf venation, such as the presence of closed loops. Simulations of auxin transport generally give rise to tree-like structures [14], a typical property of diffusion-limited growth processes [START_REF]Branching in nature[END_REF]. In contrast, it has been noted that the reticulate structure of venation networks is very similar to that of crack patterns, suggesting the possibility that the differentiation of leaf veins could be governed by mechanical stresses [11]. The inner tissues of the leaf, where the veins form, are subjected to compressive stresses by the epidermis, and one of the earliest signs of vein differentiation is the elongation of vascular cells [START_REF] Nelson | [END_REF], which could release the compressive stresses in their vicinity, providing a feedback between differentiation and the mechanical state of the tissue.

In this article, we develop a theoretical approach to study the formation of patterns induced by mechanical stresses in biological tissues. Of particular interest is the fact that these tissues are growing as the patterns form. We show how growth can be the driving force that causes the formation of new structures, and how the history of their formation can be reflected in their organization. Let us emphasize that this study mainly aims at demonstrating possible effects of mechanical stresses, and that the approach developed here would have to be integrated with biochemical mechanisms in a more complete model.

Stated in general terms, the patterning mechanism we consider is the following: compressive stresses imparted on a tissue cause it to switch to a "collapsed" state occupying a smaller volume, whereby the compression is released. In the case of leaf venation, this transition describes the differentiation of veins, and the collapsed state refers to the elongated vascular cells. An analogy can be made with the behavior of solid foams, which exhibit localized deformation upon compression [4]. Our model is not explicitly constructed according to this analogy. Instead, we consider the initial and collapsed states of the tissue as two phases of a continuous medium, which we describe using a phase field method. However, the response of this two-phase medium is analogous to that of a nonlinear elastic material, and can be analyzed with reference to the behavior of such materials 1 .

In what follows, we first present our model in a one-dimensional setting. We define a minimal model, amenable to analytical treatment, that implements the patterning mechanism described above. We consider a layered system, in which the tissue is coupled to a rigid substrate, which could represent a stiffer tissue such as the epidermis of plant leaves. Stresses are induced by a mismatch between the two layers, which can occur as a result of differential growth. If the transition between the two states of the tissue is reversible, regular patterns are obtained, in which domains of collapsed and uncollapsed tissue alternate. If, instead, the transition to the collapsed state is made irreversible to model When turning to two-dimensional systems, the tensorial nature of mechanical fields comes into play, and a much greater variety of behaviors are possible. Rather than to attempt a systematic exploration, we present a few illustrative examples of the patterns that can be obtained. In particular, we show that anisotropic collapse, which could describe the elongation of vascular cells in plant leaves, produces reticulate patterns comprised of interconnected stripes of collapsed tissue. Our simulation results suggest that hierarchical, reticulate patterns might be obtained when irreversibility is introduced. However, in the simple model presented here, these patterns are disrupted by strong residual stresses that develop in the vicinity of growing collapsed regions. We suggest that a more satisfactory model would require a more detailed description of the interplay between mechanical stresses and growth. We also contrast our model with a similar approach presented in [17].

One-dimensional model

In this section, we present our modeling approach in the context of a onedimensional system. We consider a tissue that can exist in two states, an "initial" state and a "collapsed" state that has a smaller volume. Transition to the collapsed state is triggered by compressive stresses, and releases the compression. While, in general, the two states of the tissue could differ in their stiffness (the collapse could be due to the transition to a softer state), we assume that they only differ in their rest configurations (the collapsed state has a smaller volume at rest). In the case of leaf veins, the initial and collapsed states of the tissue represent the undifferentiated and vascular tissues, respectively, and the change in rest configuration corresponds to the elongation of vascular cells as they differentiate. The tissue is described as a continuous medium, and its two states as two "phases", using a phase field model.

The phase field approach, used here in the context of mechanically induced pattern formation [17], was first introduced to describe the growth of a solid in a liquid [18,19,[START_REF] Fix | Free Boundary Problems: Theory and Applications Pitman[END_REF][START_REF] Langer | Model of pattern formation in first order phase transitions[END_REF] (reviewed in [22,[START_REF] González-Cinca | Phase field models in interfacial pattern formation out of equilibrium[END_REF]) and has later proved to be a powerful tool to describe free boundary problems in fluid flows [START_REF] Folch | [END_REF]25] and solid-solid phase transitions in alloys [26,27]. This is especially the case in 3D, where interface tracking methods are extremely difficult to implement [28,29]. Our model is similar to those developed for fracture in [30,31,32,33,34] and for solid-solid phase transitions in [26,27], which comprise phases having different mechanical properties.

A continuously varying auxiliary field φ is introduced to indicate the local state of the tissue, with the convention that φ = 0 corresponds to the initial state, and φ = 1 to collapsed regions. The behavior of the system, i.e. both its deformations and the transitions between the two states, is assumed to be governed by an energy functional E{φ, u}, where u is the displacement field. While this is natural for a physical system, it may seem rather arbitrary for a biological one. However, we will see that the resulting behavior is sufficiently general. The evolution of the system is assumed to be quasi-static, i.e. at any given time it is in a state of equilibrium that minimizes E, and satisfies the equilibrium equations

δE δu = 0 (1) δE δφ = 0. ( 2
)
For a one-dimensional system, the energy E is defined by

E = D 2 (∂ x φ) 2 + e(φ, ) + f (φ) dx. (3)
The first term in the integral is a squared gradient term, which penalizes sharp variations in the phase field. e is the density of elastic energy, which depends on the phase and on the strain = ∂ x u. f (φ) is a potential that, together with the coupling between φ and through e, governs the behavior of the phase field. Assuming that both phases (treated as linear elastic) have the same elastic modulus µ and differ only in their rest configuration, the energy defined by eq. 3 takes the form

E = D 2 (∂ x φ) 2 + µ 2 ( -r (φ)) 2 + f (φ) dx, (4)
where the rest configuration is given by

r (φ) = h(φ) r1 . (5)
The function h (defined below) satisfies h(0) = 0 and h(1) = 1, and r1 is the rest configuration of the collapsed phase. The potential f is defined by

f (φ) = αh(φ) 2 + βh(φ). (6)
h(φ) = 3φ 2 -2φ 3 is chosen such that h (0) = h (1) = 0, so that uniform phases (φ = 0 or φ = 1) are stationary states. Qualitatively, it can seen be from eq. 4 that under sufficient compression, the system will tend to switch to the collapsed state to lower its energy. The response of the system can be described by examining its uniform states of equilibrium (where φ and are uniform). In the range of parameters considered here 2 , φ has a unique equilibrium value for every value of the strain . Defining the critical strains

c0 = β µ r1 , ( 7
) c1 = µ 2 r1 + 2α + β µ r1 , (8)
2

The parameters are chosen such that the critical strains defined below satisfy c1 < c0 . φ is equal to 0 when > c0 (note that since we are considering compressive strains, and σ are generally negative), varies gradually from 0 to 1 as goes from c0 to c1 , and is equal to 1 when < c1 (for strong compressive strains). Since φ is uniquely determined by the strain , so is the stress σ = µ(r (φ)).

Our two-phase model thus behaves as a non-linear elastic material with a certain stress-strain response.

As shown on fig. 1(a), this response is non-monotonic, and the phase transitions in our model are equivalent to mechanical instabilities in an elastic material with a non-monotonic stress-strain curve [35]. The states for which the slope is negative, which correspond to values of the phase intermediate between 0 and 1, are unstable. If the system is subjected to increasing strain, it deforms uniformly until = c0 , then phase-separates. When the two phases coexist, the stress has a fixed "plateau" value σ p , as do the strains in the two phases (see fig. 1(a)). In a biological context, the plateau stress can be understood as a "homeostatic stress" restored by the formation of a collapsed region. If the strain is further increased, the proportion of the collapsed phase increases until the system is entirely collapsed, after which the strain is uniform again. Note that states intermediate between the plateau stress and the instability threshold are metastable, and in actual physical systems, the distinction between the two can be blurred by disorder [4]. This can be reproduced in our model by adding a random prefactor r(x) to the terms e(φ, ) + f (φ) in the expression of the energy (eq. 3). This creates local weak spots, so that the collapsed phase nucleates before the critical stress for an ideal, homogeneous system (compare figs. 1(a) and (b)).

As illustrated by the simple, piecewise linear stress-strain curve shown in fig. 1(a), our model provides a minimal description of a mechanism for the re- lease of mechanical stresses by localized deformation. There are just enough parameters to allow independent adjustment of the properties that will determine its pattern-forming behavior, i.e. the instability threshold, the amount of deformation associated with the collapse, and the equilibrium stress at the interface between the two phases.

In what follows, we consider the system shown in fig. 2, in which a layer of tissue described by the above model is coupled to a rigid substrate. If there is any mismatch between the two layers, mechanical stresses are induced in the tissue. Such a layered arrangement is commonly used to model various mechanical instabilities, e.g. the fragmentation of thin films [36,37], and can be likened with the layered structure of plant leaves [11], with the rigid substrate representing the epidermis.

The coupling between the two layers is described by a quadratic potential, which included into the energy (eq. 4) yields

E = D 2 (∂ x φ) 2 + µ 2 ( -r (φ)) 2 + f (φ) + k 2 (u -ηx) 2 dx. (9)
The parameter k represents the strength of the coupling, and η is a measure of the mismatch between the two layers, defined by L s = (1 + η)L, where the L s and L are the lengths of the substrate and of the tissue layer in an unstressed state.

To understand how the phase transition can be triggered, it is useful to compute the strain field when the tissue is in its initial state (φ = 0) 3 . The equation of mechanical equilibrium (eq. 1) yields

µ∂ xx u -k(u -ηx) = 0, (10)
which is linear in u and admits solutions of the form

u = ηx + Ae x λ + Be -x λ . (11)
In this equation, λ = µ/k is a characteristic elastic length scale of the system, and A and B are constants that are determined by the boundary conditions.

Here, the edges of the tissue are free of stresses, and ∂ x u = 0 at x = ±L s .

Computing A and B and deriving eq. 11 with respect to x, we find that the strain in the tissue is given by With these values, the elastic length scale is λ = .25, the instability threshold is c0 = -0.25 and the plateau stress is σ p = 0. Accordingly the strains in the two phases when they coexist are equal to the rest strains in their rest configurations, i.e. p0 = 0 and p1 = -1.

= η 1 - cosh(x/λ) cosh(L s /λ) . (12)
If the system is large compared to the elastic length scale λ, the strain is close to the strain imposed by the substrate everywhere except at the edges, and the stresses exhibit a broad plateau. If the system size is comparable to or smaller than the elastic length scale, the stresses have a more localized maximum at the center of the system, and the value of this maximum increases with the size of the system.

One way instabilities can be induced is by increasing the mismatch between the two layers (this is analogous to the situation of a shrinking film on a substrate). From the above, it appears that instabilities can also be induced by growth with the mismatch between the two layers being fixed, that is as result of uniform growth. The latter regime will always be considered when simulating patterning in growing tissues 4 .

We now turn to a numerical analysis of the evolution of the system. So far, we specified only its energy, which determines its states of equilibrium and their stability. To determine the new state of equilibrium reached after an instability occurs, its dynamics must be specified. A reasonable assumption is that mechanical relaxation is much faster than the transition between the two states of the tissue, so that at any time the system is at mechanical equilibrium:

0 = δE δu = -µ∂ x (∂ x u -r (φ)) + k(u -ηx) (13)
Dynamics need thus be specified only for the phase field. We assume the latter

∂ t φ ∼ - δE δφ = D∆φ -∂ φ e -∂ φ f. (14)
At each time step of the overall evolution of the system, i.e. of its growth, this relaxation equation is solved assuming the system is at mechanical equilibrium, until an equilibrium is reached 5 . And this state of equilibrium satisfies the equations 6 :

0 = µ∂ x (∂ x u -r (φ)) -k(u -ηx) 0 = D∆φ -∂ φ e -∂ φ f. (15)
Most of the time, no instability occurs, and the new state of equilibrium is very close to the previous one. If an instability occurs, the phase field evolves to a new state that contains a different number of collapsed regions. Fig. 3 shows a typical equilibrium configuration.

The results obtained in various conditions are shown in fig. 4. In the first example (fig. 4(a)), the tissue has a fixed size and the mismatch between the two layers is progressively increased. A limited number of collapsed regions are formed and grow to occupy the entire tissue. In the second example (fig. 4(b)), the system grows with the mismatch remaining constant (L = e t τ L 0 and L s = e t τ L s0 ). In that case, instabilities occur with time in an approximately selfsimilar cascade. Note that in both of these examples, the existing collapsed zones readjust after each new instability, and the system relaxes to a regular pattern that minimizes its energy. In the process, existing collapsed zones may change size or position, which is possible because the transition between the phases is reversible. In the case of a growing system (fig. 4(b)), this means that the pattern does not retain any trace of the history of its formation. Now, if the transition to the collapsed state represents tissue differentiation, it should be irreversible, precluding the above rearrangements. A simple way of implementing irreversibility is to make the evolution of the phase field φ irreversible. Assuming that the transition must progress to a certain point before it becomes irreversible, we allow φ to increase or decrease when it is smaller than a certain threshold (1/2), and only to increase when it is larger 7 . As shown in fig. 4(c), this yields hierarchical patterns, which reflect the history of the system. The collapsed zones grow along with the system, so that earlier-formed zones are larger than more recent ones.

So far, we assumed that the tissue was perfectly homogeneous, and it may be of interest to examine how disorder, which is inherent in biological tissues, 5 The equations of mechanical equilibrium (a set of linear equations) are solved at each step of the relaxation. The relaxation is stopped when the rate of change of the phase field falls below a certain threshold. 6 We are thus assuming that there are three well separated time scales τ mech τ diff τ growth , where τ mech , τ diff , and τ growth are the characteristic time scales corresponding to mechanical relaxation, cell differentiation, and growth, respectively.

7 Such an ad hoc rule would clearly be unsatisfactory to describe a physical system, but in the case of biological systems, it can be seen as a necessary simplification of the complex phenomena involved. would affect the outcome of the model. As described earlier, disorder can be incorporated as a multiplicative noise term in the energy of the system. Instabilities then no longer occur exactly at the locations where the stresses are largest. Fig. 4(d) shows an example of the resulting irregular patterns. Consistent with theoretical arguments and experimental results concerning crack patterns [37,38], we find that the effect of disorder depends crucially on how the distance between collapsed regions compares with the elastic length scale λ. If this distance is much larger than λ, then the stresses exhibit a broad plateau and the location of new collapsed zones is very variable, being governed essentially by the weak spots due to disorder. On the other hand, if the distance between collapsed regions is comparable to or smaller than λ, then the stresses have a more localized maximum and the location of new collapsed zones is close to the location of this maximum. The latter regime is thus more favorable for obtaining regular and reproducible patterns. It is worth noting that in growing plant leaves, the distance between veins is of the same order as the thickness of the leaf [START_REF] Nelson | [END_REF] (which sets the characteristic length scale for elastic coupling between veins).

Two-dimensional model

We now turn to the extension of the model to two-dimensional systems. In this context, the one-dimensional patterns described in the previous section can be interpreted as parallel stripes of localized deformation. The more general patterns that we expect to describe would take the form of networks of stripes of different orientations, corresponding to different principal directions of deformation. It is thus clear that a scalar parameter is no longer sufficient to describe the collapsed state. Since the difference between the two phases lies only in their rest configurations, we choose here to identify the phase field φ and r , which in 2D is a tensor describing the local rest configuration. Likewise, the mismatch between the two layers is now described by a tensor η. The 2D analog of eq. 9 thus takes the form

E = D 2 |∇ r | 2 + λ 2 Tr( -r ) 2 + µ| -r | 2 + f ( r ) + k(u -ηx) 2 dS. ( 16
)
While η may be anisotropic, we assume that the tissue has no preferential direction. Accordingly, f can be expressed as a function of the invariants of r , e.g. its trace and its norm.

To obtain a behavior similar to that of the one-dimensional model, we consider potentials that have a minimum for r = 0 and a minimum of the same depth for values of r corresponding to a smaller natural volume (i.e. Tr r < 0). Among all possible potentials satisfying this constraint, we have chosen the following three illustrative examples. In the first, the secondary minimum is reached for all deformations having a given amplitude (norm), regardless of the shear. The second one specifically favors uniaxial deformations. In contrast, the third one has a secondary minimum for isotropic compression. More specifically, the first potential is defined by

f ( r ) = -(| r | -1/2) 2 + (| r | -1/2) 4 , (17)
which has minima for r = 0 and for | r | = 1 (see fig. 5(a)). The second and third examples are defined as polynomials of the invariants of r , i.e.

f ( r ) = c ij (Tr r ) i | r | 2j . ( 18
)
The two sets of coefficients used are c 20 = . In contrast with the one-dimensional system show in fig. 2, in which the edges were free of stresses, most of the 2D simulations were done with periodic boundary conditions, which allow the mechanical equilibrium equation to be efficiently solved in Fourier space. In this case, the phase field is initially chosen to be random, to break the symmetry of the system. The initial conditions (mismatch or size for a growing system) are chosen such than an instability immediately develops, preventing relaxation to a uniform state.

As in the one-dimensional case, we first consider reversible evolution in a non-growing system (fig. 6. Here, the mismatch between the two layers is directly set to its final value, and the system is allowed to relax to equilibrium. The first two potentials (shown in figs. 5(a) and (b)), which allow or favor uniaxial compression, yield very similar patterns of interconnected stripes, which could be likened to crack networks in directional growth systems such as basalt columns [39]. In contrast, the potential shown in fig. 5(c), which favors isotropic compression, produces islands of localized deformation arranged in an hexagonal lattice, which are similar to some structures observed in solid-solid phase transitions [27]. The characteristic length of the above patterns, as in 1D, is set by the elastic length scale λ = µ/k. The difference between the two potentials that produce reticulate patterns becomes more apparent when considering a growing system (and a reversible evolution), which makes it possible to observe individual nucleation events leading to the formation of new collapsed regions. Indeed, in the case of the potential of fig. 5(a), new collapsed regions take the form of trijunctions that split an uncollapsed region into three 9 , while in the case of the potential of fig. 5(b), they occur as single stripes that divide an uncollapsed region across its longest extension 10 . This qualitative difference can be explained by the fact that the core of a trijunction is compressed in both directions, a state that is disfavored by the potential of fig. 5(b).

Finally, we consider irreversible evolution in a growing system, which we 9 A typical example in the case of an irreversible evolution can be seen in fig. 8(a) 10 A typical example can be seen in the case of an irreversible evolution can be seen in fig. 8(c found to lead to hierarchical patterns in 1D. In that case, irreversibility was implemented by forcing the evolution of the scalar phase field to be monotonic once a certain threshold was exceeded. A straightforward extension to two dimensions would be to apply the same rule to the trace of the rest configuration tensor r . However, this would allow the principal directions of the rest configuration in the collapsed phase to vary (and simulations indicate that this can indeed occur). Instead, we wish the establishment of these special directions to be definitive, as is the transition to the collapsed state (think of the elongation of vascular cells in plant leaves). Accordingly, we impose that, once the compression in some direction has exceeded a certain threshold, it can no longer decrease 11 .

As in 1D, we consider a growing system with a fixed mismatch, and find that hierarchical structures are obtained. As shown on fig. 7, the shape of the patterns depends both on the potential used and on the shape of the growing domain (for the earliest structures that are formed). However, we find that the simple model presented here produces unexpected artifacts, such as regions of lower compression within the collapsed zones, and open-ended collapsed regions that fail to reconnect with previously formed structures.

As illustrated by fig. 8, these effects are caused by tensions that accumulate in the collapsed structures as they grow. When the transition to the collapsed state is reversible, the collapsed regions tend to an equilibrium width such that compressive stresses are relaxed in their neighborhood (i.e. for an ideal, infinitely long and straight stripe, the stress at the interface has a fixed value equivalent to the plateau stress σ p of the one-dimensional model). In contrast, in the irreversible case, the width of each collapsed region grows along with the system and exceeds its equilibrium width. Instead of just releasing the compression in their vicinity, these growing collapsed regions give rise to tensile stresses. As can be seen on fig. 7, these tensile stresses both destabilize the collapsed regions, and prevent the progression of new collapsed regions towards them. It must be noted that such tensile stresses are also generated in the 1D model. However, they do not have such dramatic effects in that case, because new collapsed structures always form away from existing ones.

Discussion

In this article, we explored a possible mechanism for biological patterning by mechanical stresses. This mechanism, which involves the stress-mediated transition of a tissue between two states, was shown to yield a response similar to 11 To this end, we define a "maximum deformation tensor" m that keeps track of the maximum past compression in any direction, i.e. the tensor r (t)m(t ) is positive for all times t < t . From this maximum deformation tensor, we define an "irreversible deformation tensor" i . In a basis that diagonalizes m, i = diag(f (λ 1 ), f (λ 2 )), where m = diag(λ 1 , λ 2 ). The function f is defined by f (x) = 0.9x if x < -0.5 and f (x) = ∞ if x > 0.5, so that compressions in excess of .5 are stored in i . Irreversibility is implemented by imposing that the tensor ri is negative. The factor .9 in the definition of f was introduced to prevent numerical instability. that of a non-linear elastic material, justifying an analogy with mechanical instabilities. The formation of patterns was driven by coupling the tissue to a rigid substrate, which could represent a stiffer tissue such as the epidermis of plant leaves. As in the mechanical instabilities of layered physical systems, such as the wrinkling of a film bound to a substrate, this elastic coupling introduces a characteristic length scale that governs the size of the patterns.

We first analyzed a one-dimensional system. In this case, the model involves a small number of parameters, and it is possible to describe the equilibrium patterns analytically. Numerical simulations showed that regular patterns, which are largely independent of the history of the system, are obtained when the transition between the two states of the tissue is reversible. In contrast, when this transition is made irreversible to represent tissue differentiation, the history of the system is retained in its final state. In a growing system, new structures keep forming over time, yielding hierarchical patterns. This is very similar to the development of leaf venation networks, in which veins of different orders form successively over the course of leaf growth [START_REF] Nelson | [END_REF].

In two dimensions, because of the tensorial nature of elastic fields, the potential that characterizes the preferred states of the tissue can have many different forms. While we did not carry out a systematic analysis, the different examples considered suggest that reticulate patterns can readily be obtained when the potential allows or favors uniaxial deformations of the tissue. In contrast, islands of deformed tissue are obtained with a potential that favors isotropic deformations. This suggests that a mechanical model of leaf venation patterning could account for both wild-type patterns and the disconnected vasculature observed in certain mutants [40].

A very similar model of leaf vein patterning was recently proposed in [17]. In that approach, however, the different states of the tissue are characterized by a single scalar order parameter. The transition to the collapsed state is associated with an isotropic deformation, and the collapsed tissue is assumed to have a reduced shear modulus to allow the formation of reticulate patterns. The elongation of vascular cells is described as a change of their rest volume combined with a strong elastic shear. One benefit of this approach is that it involves a much smaller number of parameters. However, it is unclear that it gives a realistic description of the mechanical state of vascular cells. While it comes at the cost of a greater complexity, we find it more appropriate to describe the elongation of vascular cells as a change in their rest shapes. This also makes it possible to address the effect of variations in these rest shapes on the overall pattern. For instance, as mentioned above, our results suggests that failure to elongate properly can lead to disconnected patterns such as are observed in some mutants.

From the preceding results, we expected that the introduction of irreversibility in a growing two-dimensional system would lead to hierarchical, reticulate patterns. However, we found that strong tensile stresses develop in collapsed regions as they grow beyond their initial, equilibrium size, and dramatically affect the evolution of the system. This is a consequence of the very simple growth laws assumed in our model. The rest configuration of the tissue changes when it switches between its two states, but growth in each of the two states is uniform and proceeds at a constant rate, allowing residual stresses to accumulate. In reality, the growth of biological tissues depends on the mechanical stresses to which they are subjected [41,42], which can be expected to limit the buildup of residual stresses. Incorporating this dependence would be essential to further investigate biological patterning by mechanical forces. Another perspective would be to integrate the mechanical approach developed here with the biochemical factors of tissue differentiation. For instance, the differentiation of vascular cells in plant leaves involves the expression and polar localization of auxin carrier proteins [13], and it would be of interest to investigate how these processes are connected with the elongation of vascular cells.

  differentiation, hierarchical patterns are obtained.

4 PageFigure 1 :

 41 Figure 1: (a) Sketch of the stress-strain curve of the model. The dashed segment corresponds to unstable states, which phase-separate. The dotted horizontal line indicates the plateau stress observed when the two phases coexist. (b) Stressstrain response of a simulated inhomogeneous system.

5 PageFigure 2 :

 52 Figure 2: Diagram of the system. A layer of tissue (top) is coupled to a rigid substrate (bottom).

Figure 3 :

 3 Figure 3: A simulated pattern with three collapsed zones. The solid line and the dashed line represent the strain and the stress, respectively. The parameters of the model are µ = 1, r1 = -1, α = -0.25, β = 0.25, D = 4 × 10 -4 , k = 16, L = 1, and η = -.5.With these values, the elastic length scale is λ = .25, the instability threshold is c0 = -0.25 and the plateau stress is σ p = 0. Accordingly the strains in the two phases when they coexist are equal to the rest strains in their rest configurations, i.e. p0 = 0 and p1 = -1.

Figure 4 : 10 Page

 410 Figure 4: Patterns obtained in a one-dimensional system. These plots show the evolution of the phase field, with the collapsed regions appearing in black. Model parameters are given in fig. (3). (a) Increasing mismatch between the two layers, with the size of the tissue being fixed (L = 1). (b-d) Growing system with a fixed mismatch (η = -0.275). The initial size of the tissue is L 0 = 0.25. (b) Reversible evolution. Note that the collapsed regions appear to shrink only because the total size of the system increases. (c) Irreversible evolution. (d) Same as (c) with disorder.

Figure 5 :

 5 Figure 5: Gray scale plots of the different potentials f ( r ) used. 1 and 2 are the eigenvalues of r . Darker regions correspond to lower values. (a) The potential defined by eq. 17 favors compressive strains of a given amplitude, regardless of the shear. (b) The potential defined by eq. 18 and the first set of coefficients given in the text favors uniaxial compression. (c) The potential defined by eq. 18 and the second set of coefficients favors isotropic compression.

  4, c 01 = .2, c 30 = .7, c 11 = .5, c 40 = .4, c 21 = .1, c 02 = .1, which corresponds to minima for uniaxial deformations (see fig. 5(b)) and c 20 = .25, c 01 = .5, c 30 = 1, c 11 = 0, c 40 = .375, c 21 = 0, c 02 = .5, which leads to a minimum for isotropic compression (see fig. 5(c)) 8 .

Figure 6 :

 6 Figure 6: Typical patterns observed in a layered two-dimensional system with the different potentials. Periodic boundary conditions were used and the initial state was randomly perturbed. Gray levels correspond to the trace of the rest configuration r , with darker areas corresponding to collapsed regions. Potentials that allow or favor uniaxial compression (a and b, corresponding to the potentials of figs. 5(a) and 5(b), respectively) yield reticulate patterns, while the potential that favors isotropic compression produces discontinuous domains arranged in a hexagonal lattice (c, corresponding to the potential of fig. 5(c)).

) 12 Page

 12

Figure 7 :

 7 Figure 7: Patterns generated by irreversible evolution in growing systems with different potentials and boundary conditions (periodic boundary conditions or circular domains). (a) and (c) were obtained with the potential of fig. 5(a), (b) and (d) with the potential of fig. 5(b).As expected, hierarchical structures are formed. However, they are disrupted by the tensions that accumulate within and around the collapsed regions as they grow (see text and fig.8). Notice the lighter spots or stripes within the collapsed regions.

14 PageFigure 8 :

 148 Figure 8: (Color online) Effect of irreversibility on the mechanical state of the system and pattern formation. All figures represent growing systems with periodic boundary conditions. (a) and (b) were obtained with the potential of fig. 5(a), (c) and (d) with the potential of fig. 5(b). Evolution was reversible in (a) and (c), irreversible in (b) and (d). Collapsed regions (defined by thresholding the value of Tr( r )) are shown in blue.Gray levels indicate compressive stresses (lighter regions are more strongly compressed), and regions appearing in black are under tension. In the reversible case, the stresses decrease near the collapsed regions, yet generally remain compressive. In the irreversible case, this is also true for recently formed structures. In contrast, the larger, older structures generate tensile stresses. These tensile stresses hinder the progression of newly formed structures, yielding disconnected patterns.

Conversely, localized deformation in non-linear elastic materials can be interpreted as a first-order phase transition.

Similar calculations could be carried out for the states of equilibrium in which the two phases coexist if the width of the interfaces is neglected. 6 Page 6 of 18 http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

If the two layers were instead assumed to grow at different rates, then the mismatch would become arbitrarily large with time.

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Note that the potential defined by eq. 17 cannot be written in this form.

http://mc.manuscriptcentral.com/pm-pml

Acknowledgements

This work was supported by EC NEST project MechPlant. The authors would like to thank one anonymous referee for detailed remarks and helpful suggestion on a previous version of this manuscript.