# On the Short-Time Compositional Stability of Periodic Multilayers

## I Introduction

Modern techniques of thin film deposition permit the preparation of multilayers with nearly arbitrary concentration profiles. The knowledge of the stability of artificial multilayers at elevated temperatures is of great practical interest. With increasing atomic mobility, compositional changes occur where the diffusion distance is determined by the mobility and the available time. In the present paper, the case of coherent periodic multilayers consisting of two immiscible components A and B is considered. Accordingly, competing driving forces for compositional changes in these multilayers are the reduction of the energy of mixing of the two components and the reduction of the interfacial energy of interphase boundaries.

To evaluate the thermal stability of multilayers, we investigate in the following the early stage of compositional changes perpendicular to the individual layers within a one-dimensional model ignoring lateral perturbations of the layered structure as well as boundary effects at the multilayer surface and the interface to the substrate. The neglect of boundary effects is justified as long as the characteristic length of atomic diffusion for the considered time scale is small compared to the total multilayer thickness. The compositional evolution is studied within the framework of the nonlinear Cahn-Hilliard diffusion equation . Correspondingly, the multilayer is described by a smooth concentration profile where individual layers are separated by diffuse interphase boundaries. The continuum description seems questionable in the limiting case of very thin layers of a few monolayers only. However, many predictions of the continuum approach agree qualitatively with those of a detailed analysis of an appropriate lattice model by Hillert .

The composition in multilayers evolves quite differently, depending on the initial concentration profile. This is illustrated by the three examples in Fig. 7. The as-prepared multilayers in Fig. 7 were assumed as periodic A/B layer stacks of pure individual layers, with the exception of case (c) which exhibits a small thickness perturbation. The curves show the concentration of component B as a function of position. They were obtained as numerical solution of the Cahn-Hilliard diffusion equation. Depending on the initial layer thicknesses, different cases can be distinguished: (a) dissolution of very thin layers despite the immiscibility of the components (phase separation can occur only on a larger length scale), (b) relaxation of a periodic structure to a stationary state with smooth concentration profile, and (c) rapid smoothening of the profile followed by slow thinning of thinner layers up to their complete dissolution. The main driving force for the composition changes in Fig. 1 is the minimisation of interface energy. This happens even at the expense of volume energy so that the minimal and maximal values of the concentration profile generally differ from the concentrations of the corresponding bulk phases. A series of experimental findings  reported recently seems to be related to these peculiarities.

During mechanical alloying, the formation of nonequilibrium supersaturated phases of immiscible components as e. g. Ag–Cu  and Co–Cu  was observed. Due to the small layer thicknesses of the lamellar structure which arises during ball milling, a partial intermixing of the two components could be energetically favourable. Similarly, an enhanced solution of carbon in nickel layers was detected in Ni/C multilayers prepared by pulsed laser evaporation with individual layer thicknesses of only a few nanometers . Another surprising observation is the formation of a mixed Co–Cu phase during annealing of Co/Cu multilayers despite the immiscibility of Co and Cu . The metastable mixed phase formed obviously due to the large excess of interface energy in the multilayer. A strong intermixing was also found during deposition of a few monolayers Ni onto Au . In this case, besides interface energy, the large elastic energy of the strained Ni layer is an additional driving force for intermixing. For a better understanding of these experimental findings, the present theoretical work deals with the effect of the high portion of interface energy on the early composition evolution in nanoscale multilayers.

The simulations in Fig. 7 as well as numerical and analytical investigations by other authors  reveal that concentration profiles evolve in characteristic stages: (i) the relaxation to layered quasi-stationary states or the complete vanishing (dissolution) of very thin layers takes place comparatively rapidly; (ii) at a longer time scale, a slow ripening process involving diffusion between distant layers occurs, i. e. a thinning of the thinnest layers and a corresponding thickening of the thicker ones. In the following, the slow ripening process is referred to as long-time composition evolution, whereas the relaxation of an arbitrary periodic concentration profile to a stationary solution of the Cahn-Hilliard diffusion equation is referred to as short-time evolution (Fig. 7b). Also, the rapid dissolution of layers as shown in Fig. 7a is considered as a particular case of the short-time evolution.

The aim of the present work is to analyse the conditions under which either a rapid dissolution of thin layers occurs initially (Fig. 7a) or a stationary periodic concentration profile evolves (Fig. 7b). To this end, we focus on one-dimensional stationary solutions of the Cahn-Hilliard diffusion equation which are characterised by the multilayer period $`d`$ and the average composition $`\overline{c}`$ of the multilayer. These two parameters are usually controlled by the multilayer preparation and are conserved during relaxation towards quasi-stationary states. From the following analysis, a $`\overline{c}`$-$`d`$ diagram results which shows the existence region as well as the stability properties of the stationary solutions. In the present paper, a periodic solution is called (globally) stable if its energy is smaller than those of all other concentration profiles with the same periodicity $`d`$ and average composition $`\overline{c}`$. All other stationary solutions are metastable or unstable with respect to perturbations conserving the periodicity and average composition. As outlined by Langer , stationary periodic solutions are always unstable against perturbations of the periodicity (cf. also Fig. 7c). However, as mentioned above, these periodicity perturbations develop in general on a longer time scale and are not considered in this work. Also, small lateral composition perturbations as well as a roughening of interfaces between strained inidvidual layers  are not analysed here because it is expected that such perturbations develop on a longer time scale.

The paper is organised as follows. After the description of the model in Sect. II, stationary solutions of the Cahn-Hilliard diffusion equation are considered in Sect. III. Sect. IV deals with a specific model for the Gibbs free energy which allows an analytical calculation of stationary concentration profiles. In Sect. V, the stability of these concentration profiles is investigated. Finally, the results are discussed and summarised.

## II Model

The Gibbs free energy density $`f(c)`$ of a binary A-B system as a function of the uniform concentration $`c`$ (mole fraction) of component B is to exhibit two minima. The equilibrium concentrations, $`\alpha `$ and $`\beta `$, of large coexisting phase regions (strictly two half-spaces) are determined by the common tangent construction (Fig. 7). In the present work, the case of multilayers with planar interfaces is considered. Following Cahn and Hilliard , the free energy per unit area of a system with concentration $`c(x)`$ varying in one dimension is described by

$$F[c]=𝑑x\left[f(c)+\kappa \left(\frac{dc}{dx}\right)^2\right].$$
(1)   
The second term on the right-hand-side of (1) represents the energy contribution due to a concentration gradient where $`\kappa `$ is the gradient energy coefficient.

The interdiffusion flux in the system is given by $`j(x)=\stackrel{~}{M}(\delta F/\delta c)/x`$ where $`\stackrel{~}{M}`$ is the atomic mobility. Together with the continuity equation $`c/t+\mathrm{\Omega }j/x=0`$, the following nonlinear diffusion equation results 

$$\frac{c}{t}=M\frac{^2}{x^2}\left(\frac{\delta F[c]}{\delta c(x)}\right)=M\frac{^2}{x^2}\left[f^{}(c)2\kappa \frac{^2c}{x^2}\right]$$
(2)   
where $`M\mathrm{\Omega }\stackrel{~}{M}`$ with $`\mathrm{\Omega }`$ the atomic volume. For simplicity, the atomic volume of the two components has been assumed to be equal and the composition dependence of $`M`$ has been omitted. Starting from any initial concentration profile, the further evolution can be calculated numerically from (2) (see e. g. Fig. 1). However, in view of the great practical importance of quasi-stationary concentration profiles, we will analyse the stationary solutions of the Cahn-Hilliard diffusion equation (2) more systematically in the following.

## III Stationary solutions

Equilibrium concentration profiles are determined by the extrema of the free energy under the constraint of particle conservation. This leads to the variational problem

$$\frac{\delta }{\delta c(x)}\left(F[c]\mu 𝑑xc(x)\right)=0$$
(3)   
with the result

$$f^{}(c)2\kappa \frac{d^2c}{dx^2}\mu =0$$
(4)   
(the prime denotes the derivative with respect to $`c`$). Comparison of (4) and (2) reveals that solutions of (4) are also stationary solutions of the diffusion equation (2). The Lagrange multiplier $`\mu `$ is identified as interdiffusion potential. When $`\mu =\delta F/\delta c`$ is uniform, the particle flux vanishes. Integration of (4) leads to

$$\kappa \left(\frac{dc}{dx}\right)^2=f(c)\mu c+KD(c)$$
(5)   
where $`K`$ is an integration constant. The last equality in (5) defines the function $`D(c)`$ used in the following. In general, the physically relevant solutions of (5) are periodic concentration profiles $`c(x)`$ oscillating between a minimal value $`a`$ and a maximal value $`b`$ (cf. Figs. 2 and 3). The extrema $`a`$ and $`b`$ are related to the parameters $`\mu `$ and $`K`$ by the conditions $`D(a)=D(b)=0`$ which lead to

$$\mu =\frac{f(b)f(a)}{ba},K=\frac{f(b)af(a)b}{ba}.$$
(6)   
Further integration of (5) yields the inverse function of the concentration profile

$$x(c)=_{c_0}^c𝑑cI(c)$$
(7)   
with $`I(c)=\sqrt{\kappa /D(c)}=dx/dc`$. The integration bounds have been chosen in such a way that the origin of the $`x`$-coordinate is located at an interphase boundary defined by $`c=c_0(a+b)/2`$. Thus, equation (7) represents the concentration profile in half a period from $`c=a`$ at $`x=d_a/2`$ to $`c=b`$ at $`x=d_b/2`$ (Fig. 7). The individual layer thicknesses $`d_a`$ and $`d_b`$ of the two phase regions, briefly called phase ’a’ and phase ’b’, are given by

$$d_a=2_a^{c_0}𝑑cI(c),d_b=2_{c_0}^b𝑑cI(c).$$
(8)   
Similarly, the multilayer period length $`d=d_a+d_b`$ and the mean composition $`\overline{c}`$ are given by

$$d=2_a^b𝑑cI(c),\overline{c}=\frac{2}{d}_{d_a/2}^{d_b/2}𝑑xc(x)=\frac{2}{d}_a^b𝑑cI(c)c.$$
(9)   
For given concentrations $`a`$ and $`b`$, the concentration profile $`c(x)`$ can be calculated directly from (7). However, from the experimental point of view, $`a`$ and $`b`$ are not known a priori. Usually, components A and B are deposited consecutively with fixed individual layer thicknesses. During the early annealing stage, a smooth concentration profile develops which is similar to the stationary profiles derived here. The arising concentrations $`a`$ and $`b`$ are determined by equations (9) where the mean composition $`\overline{c}`$ and period length $`d`$ are given. The solution of equations (9) for $`a`$ and $`b`$ is, however, not always unique. If there is more than one solution, one has to compare the free energies of different solutions in order to find that with the lowest one. From (1) and (5), the free energy of one multilayer period results as

$$F_p=4\sqrt{\kappa }_a^b𝑑c\sqrt{D(c)}+\left[\frac{b\overline{c}}{ba}f(a)+\frac{\overline{c}a}{ba}f(b)\right]d.$$
(10)   
An important intrinsic length of the present problem is the width of the interphase boundary $`\xi `$ defined by

$$\xi =(ba)\frac{dx(c)}{dc}|_{c=c_0}=(ba)\sqrt{\frac{\kappa }{D(c_0)}}.$$
(11)   
The second equality follows from (5). In the limiting case of spatially extended phases ($`d_a,d_b\xi `$; i. e. $`a\alpha `$, $`b\beta `$), the interface width becomes

$$\xi =(\beta \alpha )\sqrt{\kappa /f_0}(\beta \alpha )l_0$$
(12)   
with $`f_0D((\alpha +\beta )/2)`$. $`f_0`$ characterises the height of the free energy wall of the $`f(c)`$ curve referred to the common tangent (cf. Fig. 7). The last equality in (12) defines the length unit $`l_0`$ used in the following.

For very thin individual layers, comparable with the width of the interphase boundary, the concentrations in the middle of the layers, $`a`$ and $`b`$, differ from those of the corresponding bulk phases $`\alpha `$ and $`\beta `$ ($`\alpha <a<b<\beta `$, Fig. 7) because the common tangent construction does not apply to thin layers. The concentrations $`a`$ and $`b`$ define a secant with the $`f(c)`$ curve as shown in Fig. 7. An estimate of the difference between the concentrations $`\beta `$ and $`b`$ is given by

$$\frac{\beta b}{\beta \alpha }=\rho _b\mathrm{exp}(d_b/2\xi _b)$$
(13)   
(see Appendix A), where $`\xi _b\sqrt{2\kappa /f^{\prime \prime }(\beta )}`$ and $`\rho _b`$ is a numerical factor of the order of unity. An analogous formula applies to the difference $`a\alpha `$. In the limiting case $`d_b\xi _b`$, a factor of $`\rho _b=2`$ has been calculated for the special composition dependence of the Gibbs free energy considered in the following section. The estimate (13) clearly reveals that concentrations $`\beta `$ and $`b`$ differ significantly when the layer thickness $`d_b`$ approaches the characteristic length $`\xi _b`$.

## IV Special case: $`c^4`$–model

To simplify the calculation of the concentration profile (7), we consider the case where the free energy as a function of concentration can be represented as a polynomial of the fourth power

$$f(c)=A(c\alpha )^2(c\beta )^2+B(c\alpha )+C.$$
(14)   
The parameters $`B`$ and $`C`$ turn out to be unimportant for the composition evolution. The characteristic energy unit $`f_0`$ results as $`f_0=A(\beta \alpha )^4/16`$. The parameters $`\alpha `$ and $`\beta `$ in (14) coincide with the equilibrium concentrations of spatially extended phases corresponding to the common tangent construction (Fig. 7). The characteristic length $`\xi _b`$ in (13) is obtained as $`\xi _b=(\beta \alpha )l_0/4`$.

For the case (14), briefly called ’$`c^4`$-model’, the inverse concentration profile (7), the period length, and the mean composition (9) can be expressed by elliptic integrals

$$x(c)=\omega \sqrt{\frac{\kappa }{A}}[F_e(\varphi (c),m)F_e(\varphi (c_0),m)],$$
(15)   
$$d=2\omega \sqrt{\frac{\kappa }{A}}K(m),$$
(16)   
$$\overline{c}=a+\frac{2}{\omega }Z(\varphi _Z,m)$$
(17)   
with $`\omega =2/\sqrt{(b_1a)(ba_1)}`$, $`m=(ba)(b_1a_1)/((ba_1)(b_1a))`$, $`\varphi _Z=\mathrm{arcsin}\sqrt{(b_1a)/(b_1a_1)}`$, and $`\varphi (c)=\mathrm{arcsin}\sqrt{(ba_1)(ca)/((ba)(ca_1))}`$. $`K(m)`$ and $`F_e(\varphi ,m)`$ are the complete and incomplete elliptic integrals of first kind, and $`Z(\varphi ,m)`$ is the Jacobi zeta function . $`a`$ and $`b`$ are the concentrations in the middle of the individual layers, and $`a_1`$ and $`b_1`$ are the further two intersections of the $`f(c)`$-curve with the secant shown in Fig. 7. $`a_1<a<b<b_1`$ are also the four zeros of the function $`D(c)`$ defined in (5). Choosing $`a`$ and $`b`$, the other two zeros are given by

$$a_1=\alpha +\beta c_0w,b_1=\alpha +\beta c_0+w$$
(18)   
with $`w\sqrt{2[c_0(\alpha +\beta )\alpha \beta ]+ab3c_0^2}`$ and $`c_0(a+b)/2`$. Analytical expressions for the stationary solutions of the Cahn-Hilliard diffusion equation in the case of the $`c^4`$-model have been given previously in terms of Jacobian Elliptic functions by Tsakalos  for the symmetric case $`d_a=d_b`$ and by Novick-Cohen and Segel for the general case .

Fig. 7 shows a series of interesting quantities as a function of the mean concentration $`\overline{c}`$ for fixed period length $`d`$, calculated within the $`c^4`$-model. In a certain $`\overline{c}`$$`d`$ region, two stationary periodic solutions have been found. The corresponding branches in Fig. 7 are denoted by #1 and #2, respectively. The free energy $`F_p`$ of the periodic concentration profile with smaller concentration variation $`ba`$ (solution 2) is higher than that of solution 1 (Fig. 7a). For small values of $`\overline{c}`$, the free energy of the homogeneous concentration, $`F_h=f(\overline{c})d`$, is lower than that of the periodic solution 1. The values of the concentrations $`a`$ and $`b`$ are shown in Fig. 7b. A striking feature of these plots is the existence of a minimal value $`\overline{c}_a^M`$ of the mean concentration. At the minimum $`\overline{c}_a^M`$, solutions 1 and 2 merge. The corresponding layer thickness of phase ’b’ is denoted by $`d_b^M`$ (Fig. 7c). It correlates roughly with the minimal value of $`d_b`$ in this case. As discussed in the next section, solution 2 does not appear for multilayer period lengths $`d`$ below a critical value.

## V Stability of stationary solutions

In the following, let us consider the stability of periodic concentration profiles $`c(x)`$ which are obtained as solutions of equation (5). These concentration profiles are also stationary solutions of the Cahn-Hilliard diffusion equation (2). They are stable against any infinitesimal perturbation $`\delta c(x)`$ if the second variation of the free energy

$$\delta ^2F=𝑑x\left[f^{\prime \prime }(c)(\delta c)^2+2\kappa \left(\frac{d\delta c}{dx}\right)^2\right]$$
(19)   
is positive definite.

At first, the stability of the homogeneous concentration $`c(x)=\overline{c}`$ is considered. Without any restriction, the perturbation $`\delta c`$ can be represented as a Fourier series. From (19), it is evident that the stabilising influence of the gradient term (second term on the right-hand-side of (19)) is stronger the shorter the wavelength of the perturbation is. Considering a periodic perturbation with period $`d`$, i. e. $`\delta c\mathrm{sin}(2\pi x/d)`$, we find stability $`\delta ^2F0`$ for period lengths

$$d<d^S(\overline{c})=2\pi \sqrt{\frac{2\kappa }{f^{\prime \prime }(\overline{c})}}.$$
(20)   
$`d^S(\overline{c})`$ is equal to the smallest wavelength of spinodal decomposition obtained from a stability analysis of diffusion equation (2) . The inverse function of $`d^S(\overline{c})`$ exhibits two branches which are denoted by $`\overline{c}_a^S(d)`$ and $`\overline{c}_b^S(d)`$ (Fig. 7). For the $`c^4`$-model, one obtains

$$\overline{c}_{b/a}^S(d)=\frac{\alpha +\beta }{2}\pm \frac{\beta \alpha }{2\sqrt{3}}\left[1\frac{\pi ^2}{2}(\beta \alpha )^2\left(\frac{l_0}{d}\right)^2\right]^{1/2}.$$
(21)   
Within the interval $`\overline{c}_a^S<\overline{c}<\overline{c}_b^S`$, the homogeneous solution is unstable and spinodal decomposition takes place. According to (20), the maximum of $`f^{\prime \prime }(\overline{c})`$ yields a minimal period length $`d_{min}`$ below which the homogeneous solution is stable for all values of $`\overline{c}`$. For the $`c^4`$-model, $`d_{min}=\pi (\beta \alpha )l_0/\sqrt{2}`$ follows.

Since within the region $`\overline{c}_a^S<\overline{c}<\overline{c}_b^S`$ the homogeneous solution is unstable, there must be a stable periodic solution. The investigation of the existence region and stability of periodic solutions is more complicated than that of the homogeneous one. In the following, the results of a numerical calculation of the free energies (10) belonging to the solutions (15) of the $`c^4`$-model are summarised (Fig. 7). It is expected that the qualitative features are the same for other double-well potentials $`f(c)`$. Similar results have been obtained by Hillert  for a lattice model.

The behaviour of periodic solutions changes qualitatively in dependence on the multilayer period length $`d`$. The stability of concentration profiles can be discussed in a convenient way by including non-stationary states with nonequilibrium amplitudes $`ba`$ in the consideration. Although the parameter $`ba`$ does not give a complete characterisation of non-stationary states, it is an appropriate quantity to illustrate the stability behaviour of periodic concentration profiles. In Fig. 7, the free energy $`F`$ of concentration profiles is sketched as a function of their amplitude for different mean compositions and for two multilayer period lengths. Extrema of $`F`$ correspond to stationary solutions. In Fig. 7, these extrema are marked by filled and open circles denoting stable and unstable stationary solutions, respectively.

For period lengths above $`d_{min}`$, but below a certain critical value $`d_c`$, the situation is represented by the curves 1 and 2 in Fig. 7a: (i) Within the region $`\overline{c}_a^S<\overline{c}<\overline{c}_b^S`$, there is a stable periodic solution corresponding to the minimum of $`F`$ (curve 1). The homogeneous solution ($`ba=0`$) belongs to a maximum of $`F`$ (more precisely a saddle point) and is therefore unstable. (ii) Outside the region $`\overline{c}_a^S<\overline{c}<\overline{c}_b^S`$, there is no stationary periodic solution (curve 2). The minimum of the free energy $`F`$ is given by the homogeneous solution which, therefore, is stable. At the concentrations $`\overline{c}=\overline{c}_{a,b}^S`$, a continuous transition between the homogeneous and the periodic stationary solution occurs. This transition will be referred to as second order ’phase transition’ in the following despite the fact that the resulting ’phases’ are only stable in the sense discussed in the Introduction.

For larger period lengths $`d>d_c`$, the behaviour is more complex as illustrated in Fig. 7b. Within the region $`\overline{c}_a^S<\overline{c}<\overline{c}_b^S`$, the situation is the same as for $`d<d_c`$ (curve 1). Outside certain limiting values $`\overline{c}_{a,b}^M`$ of the mean concentration, there is no stationary periodic solution (curve 4). However, within the regions $`\overline{c}_a^M<\overline{c}<\overline{c}_a^S`$ and $`\overline{c}_b^S<\overline{c}<\overline{c}_b^M`$, two stationary periodic solutions exist (curves 2 and 3; cf. also Fig. 7). One solution (solution 1) corresponds to a minimum of the free energy $`F`$ and the other one (solution 2) to a maximum. Thus, solution 2 is unstable. Choosing solution 2 as initial condition for the concentration profile, the numerical solution of the diffusion equation (2) reveals a rapid change of the profile. Depending on the initial fluctuations, it evolves either into the periodic solution 1 or into the homogeneous solution. The latter one corresponds also to a minimum of $`F`$. Whether the periodic solution 1 or the homogeneous one is more stable depends on the corresponding values of the free energy $`F_p^{(1)}`$ and $`F_h=f(\overline{c})d`$. At a certain concentration $`\overline{c}=\overline{c}_{a,b}^T(d)`$, for which $`F_p^{(1)}=F_h`$, a first order ’phase transition’ between the periodic and the homogeneous states takes place (cf. Fig. 7). Between $`\overline{c}^S(d)`$ and $`\overline{c}^T(d)`$ the homogeneous solution is metastable ($`F_p^{(1)}<F_h`$, curve 2 in Fig. 7b), whereas between $`\overline{c}^T(d)`$ and $`\overline{c}^M(d)`$ the periodic solution is metastable (curve 3 in Fig. 7b).

At $`\overline{c}=\overline{c}^S(d)`$, solution 2 disappears by merging into the homogeneous solution. With increasing distance of $`\overline{c}`$ from $`\overline{c}^S`$, the difference in the free energies of solutions 1 and 2 decreases and at the concentrations $`\overline{c}^M`$ they coincide (cf. also Fig. 7). This implies that the second variation of the free energy functional vanishes (see Appendix B) and, consequently, the periodic solutions become marginally stable at $`\overline{c}^M`$. Outside the region $`\overline{c}_a^M<\overline{c}<\overline{c}_b^M`$, no inhomogeneous stationary solution exists. Thus, $`\overline{c}_a^M(d)`$ and $`\overline{c}_b^M(d)`$ represent the minimal and maximal mean concentrations for the existence of metastable or globally stable periodic structures with given multilayer period length.

The critical value of the multilayer period $`d_c`$, which separates the regions of first and second order ’phase transitions’ in Fig. 7, can be derived analytically by means of a third-order perturbation analysis. The corresponding critical concentrations $`\overline{c}_c=c^S(d_c)`$ are determined by the equation $`(f^{\prime \prime \prime })^2+\mathrm{\hspace{0.17em}3}f^{\prime \prime }f^{\prime \prime \prime \prime }=0`$, where $`f^{\prime \prime }`$, $`f^{\prime \prime \prime }`$ and $`f^{\prime \prime \prime \prime }`$ denote the second to fourth derivative of $`f(c)`$ at $`c=\overline{c}_c`$. According to (20), the critical period length is then given by $`d_c=2\pi \sqrt{\mathrm{\hspace{0.17em}2}\kappa /(f^{\prime \prime }(\overline{c}_c))}`$. For the $`c^4`$-model, one obtains $`\overline{c}_c=[\alpha +\beta \pm (\beta \alpha )/\sqrt{5}]/2`$ and $`d_c=\sqrt{5/2}d_{min}`$.

Let us recall that the above picture on the stability of periodic solutions was obtained for the special $`c^4`$-dependence of the Gibbs free energy. Further analysis of other dependencies $`f(c)`$ is desirable to confirm the present stability diagram qualitatively and to study quantitative changes. Numerical simulations of the composition evolution in the present work were performed by means of a finite difference method described in ref.. This method is based on a semi-implicit finite difference scheme coupled with a fast Fourier transformation. For numerical details, the reader is referred to Copetti and Elliott . The spatial grid spacing in our calculations was less than 10% of the interphase boundary width of stationary states and for the time integration an adaptive step size control was applied. For example, 1024 grid points were used for the calculations in Fig. 1.

## VI Discussion

Multilayers for practical applications are usually deposited as layer stacks of nearly pure layers of components A and B with thicknesses $`d_A`$ and $`d_B`$. During annealing, a considerable interdiffusion between individual layers can occur depending on the values of the equilibrium bulk concentrations $`\alpha `$ and $`\beta `$, which change with temperature. Consider for definiteness the case of thin layers of component B between thick layers of component A ($`d_B<d_A`$). For arbitrary multilayer period lengths, the B-layers dissolve at elevated temperatures if the mean composition $`\overline{c}=d_B/d`$ (supposing $`\mathrm{\Omega }_A=\mathrm{\Omega }_B`$) is smaller than $`\alpha `$, where $`\alpha `$ increases with increasing temperature. However, layer dissolution can occur also for $`\overline{c}>\alpha `$, if the multilayer period is small enough.

The composition evolution in multilayers is controlled by the competition between volume free energy and interface energy. For mean compositions outside the region $`\overline{c}_a^M(d)<\overline{c}<\overline{c}_b^M(d)`$, no stationary periodic solutions of the Cahn-Hilliard diffusion equation have been found. This implies that individual layers of such multilayers dissolve in the early stage of annealing because the free energy gain by phase separation is too small to overcompensate the interface energy (Fig. 1a). The thickness of the ’b’-layers $`d_b^Md_b(\overline{c}=\overline{c}_a^M)`$, which corresponds to the minimal mean concentration $`\overline{c}_a^M(d)`$, is shown in Fig. 7 as a function of the period length $`d`$. For large period lengths ($`d>7l_0`$ for the example in Figs. 7 and 7), $`d_b^M(d)`$ represents the minimal layer thickness of phase ’b’ of stable concentration profiles. For smaller $`d`$ (but larger than $`d_c`$), the minimal thickness $`d_b`$ appeared in the present case at a mean concentration slightly larger than $`\overline{c}_a^M`$ (cf. Fig. 7c).

Despite the large change of the multilayer period length in Fig. 7, the characteristic layer thickness $`d_b^M(d)`$ changes only slightly and is always of the order of the characteristic length $`l_0=\sqrt{\kappa /f_0}`$. The gradient energy coefficient $`\kappa `$ can be estimated within the framework of the regular solution model  as $`\kappa \mathrm{\Delta }U/r_0`$ where $`r_0`$ is the interatomic distance and $`\mathrm{\Delta }U`$ is the energy of mixing (per atom) of an equiatomic solution ($`\overline{c}=0.5`$); typically, $`\kappa =10^{11}`$ to $`10^{10}`$ J/m. The value of $`l_0`$ considerably exceeds the interatomic distance when $`f_0\mathrm{\Delta }U/r_0^3`$, where $`\mathrm{\Delta }U/r_0^3`$ is typically in the range of $`10^8`$ to $`10^9`$ J/m<sup>3</sup>. Choosing for example $`\kappa =310^{11}`$ J/m and $`f_0=10^8`$ J/m<sup>3</sup>, one obtains $`l_00.5`$ nm which is about twice the interatomic distance.

The gradient energy coefficient changes only very weakly with temperature, whereas the parameter $`f_0`$ (cf. Fig. 2) decreases with increasing temperature, approaching zero at the critical phase separation temperature $`T_c`$. As a consequence, the length $`l_0`$ and correspondingly the thickness $`d_b^M`$, characterising the onset of layer dissolution, diverge as $`|TT_c|^{1/2}`$ at the critical temperature. Comparatively small critical temperatures in technologically important regions of a few $`100^{}`$C are found, for example, for systems with large lattice mismatch because mechanical stresses in coherent layers cause a considerable lowering of the critical temperature compared to incoherent phases .

The layer dissolution in multilayers with very small period length, which is driven by the reduction of interface energy, leads to the formation of supersaturated phases. Subsequent phase separation occurs with a larger period length on a longer time scale and can be kinetically hindered by rapid quenching. As mentioned in the introduction, such a situation could be present during mechanical alloying , when a nanoscale lamellar structure develops in the course of ball-milling, or in the case of multilayer deposition  with subsequent short-time annealing at moderate temperatures . The nonequilibrium phase formation observed in those experiments could be related to the layer dissolution discussed here. Based on an estimate of the chemical energy of coherent phase boundaries, such a mechanism has been suggested by Gente et al.  to explain the observed solid solution formation of immiscible elements due to long-time ball-milling.

## VII Summary

Stationary solutions of the one-dimensional Cahn-Hilliard diffusion equation have been analysed in dependence on the mean composition $`\overline{c}`$ of a binary system and the multilayer period length $`d`$. A $`\overline{c}`$-$`d`$ diagram has been established showing the regions of existence, metastability and global stability of stationary periodic concentration profiles, as well as of the homogeneous concentration (Fig. 7). The diagram was derived under the constraint of fixed period length of the concentration profile. Actually, periodic solutions are unstable against thickness fluctuations  which leads to layer thickness coarsening in the course of annealing. However, this happens on a longer time scale than the relaxation to quasi-stationary concentration profiles as long as the individual layer thicknesses are significantly larger than the interphase boundary width.

The present analysis revealed that very thin individual layers (typically a few monolayers for small mutual solubility) dissolve during annealing if the mean composition of artificial multilayers is lower than a critical value $`\overline{c}_a^M(d)`$. The individual layer thickness $`d_b^M`$ corresponding to the minimal value $`\overline{c}_a^M(d)`$ increases slightly with increasing multilayer period length. The layer dissolution is driven by two mechanisms: (i) interdiffusion between pure individual layers to establish the equilibrium phase concentrations and (ii) reduction of interface energy in the case of very thin layers of the order of the interphase boundary width. In a certain $`\overline{c}`$-$`d`$ region, the layered structure can exist as metastable state although the free energy of the homogeneous concentration is lower. According to the Cahn-Hilliard theory, the equilibrium composition in thin layers differs from that of the corresponding bulk phase if the layer thickness becomes comparable with the interphase boundary width .

Consideration of the present conclusions in the design of layered structures could help to improve their thermal stability or, on the other hand, to prepare new metastable phases by controlled layer dissolution. Without changing the qualitative predictions, mechanical stresses in coherent layers due to lattice mismatch can easily be included in the present one-dimensional analysis by a modification of the free energy $`f(c)`$ (see e. g. ref. ).

The evaluation of the long-time stability of multilayers requires an additional analysis of the evolution of lateral composition perturbations including the effect of stresses in individual layers. The roughening instability of interfaces between strained layers has been investigated for example in ref.. A comprehensive analysis of the composition evolution under the influence of stresses owing to lattice mismatch or to the presence of dislocations has been given recently in a series of papers by Léonard and Desai .

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft, Sonderforschungsbereich 422.     

## VIII Appendix A

To derive (13), equation (4) is linearised with respect to the concentration difference $`\delta c_b(x)=bc(x)`$. With the approximation $`f^{\prime \prime }(b)f^{\prime \prime }(\beta )`$, one obtains

$$\left(\frac{d^2}{dx^2}\frac{1}{\xi _b^2}\right)\delta c_b(x)=\frac{\mu f^{}(b)}{2\kappa }$$
(22)   
($`\xi _b^22\kappa /f^{\prime \prime }(\beta )`$). Choosing the middle of the ’b’-layer as origin of the $`x`$-coordinate and requiring $`\delta c_b(0)=\delta c_b^{}(0)=0`$, the solution of (22) is obtained as

$$\delta c_b(x)=\frac{\mu f^{}(b)}{f^{\prime \prime }(\beta )}\left(\mathrm{cosh}\frac{x}{\xi _b}1\right).$$
(23)   
Although this result was derived for small values of $`\delta c_b`$, it is extrapolated to larger values in order to get an estimate for the layer thickness $`d_b`$ defined by $`c(\pm d_b/2)=(a+b)/2`$. Assuming further $`d_b\xi _b`$ and approximating $`ba`$ by $`\beta \alpha `$, one finds

$$\delta c_b(\pm d_b/2)=\frac{ba}{2}\frac{(\beta \alpha )}{2}=\frac{1}{\rho _b}\frac{\mu f^{}(b)}{2f^{\prime \prime }(\beta )}\mathrm{exp}(d_b/2\xi _b).$$
(24)   
The correction factor $`\rho _b`$ is introduced in (24) in order to account for the error caused by the linearisation of (22) within the interface. The last equation can be rewritten using the expansion $`f^{}(b)f^{}(\beta )+f^{\prime \prime }(\beta )(b\beta )`$ and the fact that $`\mu f^{}(\beta )`$ is of the order of $`𝒪((\beta b)^2,(a\alpha )^2)`$. Neglecting these higher order terms and inserting $`\mu f^{}(b)f^{\prime \prime }(\beta )(\beta b)`$ into (24), one obtains equation (13).

## IX Appendix B

In the following, the stability of two stationary solutions for the same $`d`$ and $`\overline{c}`$, $`c_1(x)`$ and $`c_2(x)`$ differing by $`\delta c(x)=c_2(x)c_1(x)`$, is analysed for $`\delta c0`$. The two solutions fulfil equation (4) with the corresponding Lagrange multipliers $`\mu _1`$ and $`\mu _2`$. Expansion of (4) for $`c_2=c_1+\delta c`$ with respect to $`\delta c`$ and $`\delta \mu =\mu _2\mu _1`$ yields the following equation for $`\delta c`$: $`f^{\prime \prime }(c_1)\delta c2\kappa d^2\delta c/dx^2=\delta \mu `$. Using the estimate $`\delta \mu 𝒪((\delta c)^2)`$, one obtains to first order in $`\delta c`$

$$\left[f^{\prime \prime }(c_1)2\kappa \frac{d^2}{dx^2}\right]\delta c(x)=0.$$
(25)   
This equation is equivalent to marginal stability of $`c_1(x)`$. Indeed, the second variation (19) of $`F[c]`$ can be transformed by partial integration into

$$\delta ^2F=𝑑x\delta c\left[f^{\prime \prime }(c)2\kappa \frac{d^2}{dx^2}\right]\delta c.$$
(26)   
In deriving (26), the periodicity of solutions $`c_1(x)`$ and $`c_2(x)`$, and consequently of $`\delta c(x)`$, was used. Comparison of (26) and (25) leads to $`\delta ^2F=0`$. In summary, if two solutions merge, marginal stability results. This happens at the boundaries $`\overline{c}^S`$ as well as $`\overline{c}^M`$ (Fig. 7).

. Present address: Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Str. 38, D-01187 Dresden, Germany

Fig. 1 

Fig. 2 

Fig. 3 

Fig. 4 

Fig. 5 

Fig. 6 

Fig. 7