Numerical approximation of the Chan-Hillard equation with memory effects in the dynamics of phase separation

We consider the modifi ed Cahn-Hilliard equation for phase separation 
suggested to account for spinodal decomposition in deeply supercooled binary alloy 
systems or glasses. This equation contains, as additional term, the second-order 
time derivative of the concentration multiplied by a positive coefficient $\Tau_d$ (time 
for relaxation). We consider a numerical approximation scheme based on Fourier 
spectral method and perform numerical analysis of the scheme. We present results 
of numerical simulations for three spatial dimensions, and examine the stability and 
convergence of the scheme.


1.
Introduction. The aim of the present paper is to investigate the modified Cahn-Hilliard (CH) equation of a conserved parameter c as proposed by Galenko and co-workers [11,12]. Here τ D > 0 is a relaxation time, and the term τ D ∂ 2 c/∂t 2 is added in order to model non-equilibrium decompositions. Using such a model, one can predict the very first stages of decomposition when the atomic diffusion flux rapidly changes and relaxes to its local equilibrium steady state. Such description allows for prediction of the earliest stages of decomposition, intermediate regimes and the latest stages of decomposition.

NICOLAS LECOQ, HELENA ZAPOLSKY AND PETER GALENKO
to thermodynamic entropy consideration, also appeared in papers [17,18] as well as in the 1958 paper of Cahn and Hilliard [2,3,4]. The present study is concerned with developping splitting methods with such a potential and getting insured of the decrease of energy and of the conservation of the parameter c. We state the stability result along with proof and point out some related results.
2. The governing equations. We consider a binary system free of imperfections and neglect the effect of coherency strain (for crystalline or amorphous solids [2,3,20]), anisotropy (as in a blend of liquid crystals or a polymer solution [8]), and shear flow (in liquids [1,19]) which may occur in a phase separation. In this case, the process of phase separation is controlled by atomic diffusion. Whole stages of the dynamics of transition from unstable to metastable or even stable phase state can be described through introduction of the prehistory of the system evolution [11]. For the conserved dynamics of phase transition with variable c(t, r) having meaning of concentration of a second chemical component in the binary system, one may write the flux with M R (t − t * ) the memory function which adopts the correlation between the system's dynamics at the past moment of time t * and the current dynamics at the time t.
Together with the balance law ∂c/∂t = −∇ · j, Eq. (3) gives the following equation for the arbitrary, slow or fast, dynamics: Taking the free energy E in the following functional form and the memory function of the exponential form M R (t−t * ) = (D/τ D ) exp[−(t−t * )/τ D ], Eq. (4) can be rewritten as [12] τ D ∂ 2 c ∂t 2 + where v 0 is the subvolume of the system, τ D the time for relaxation of the diffusion flux j to its steady-state, M the atomic mobility, 2 c the parameter characterizing the length between the phase interfaces (which is proportional to the surface energy), f c is the derivative of the free energy density f (T, c) with respect to concentration c, and T is the temperature. In the case τ D → 0, Eq. (6) reduces to the well-known Cahn-Hilliard equation (cf. [2]) which has been widely analyzed and used to describe phase separations. This kind of phenomenon is rather complex and one of its typical feature is the so-called spinodal decomposition.
In the present form, Eq. (6) can be considered as a modified Cahn-Hilliard equation which is a partial differential equation of hyperbolic type with the term τ D ∂ 2 c/∂t 2 representing the delay in the separation of phases.
Let us assume the mobility M to be independent on concentration in Eq. (6). To obtain normalized constitutive equation, the following dimensionless variables and parameters are introduced:t = tM ∆f /l 2 D dimensionless time,∇ = ∇/l D normalized operator "nabla",f c ≡ ∂f /∂c = f c /∆f dimensionless characteristic for the free energy density,˜ 2 c = 2 c l 2 D /∆f dimensionless length scale of the interface, andτ D = τ D M ∆f /l 2 D dimensionless relaxation time. Then, Eq. (6) to be numerically solved takes the form : A natural boundary condition, following from the variation of functional (5), is given by where sub-index n indicates a vector projection on the normal vector to the outer boundary of the system. Obviously, simple boundary conditions which satisfy Eq. (8) are described by ∇.j n = 0 and ∇c = 0.
These conditions must hold at any point on the outer surface around the volume v 0 . These boundary conditions lead to the two very important properties which will be our guide to know whether our causal-CH numerical solutions are convergent, namely the conservation of the total mass M t of the system and the dissipation or decrease of the total energy : 3. Unconditionally stable algorithm for the Cahn-Hilliard equation. For the Cahn-Hilliard (CH) equation, with a free energy functional described with a doublewell structure, and a gradient-squared term associated with an energy cost to an interface between two phases, several unconditionally stable algorithms were developped [9,10,25].Considering the simple case where the double-well potential is represented in a polynomial form, the CH equation takes the form : A mix of explicit and implicit terms in the update is introduced. For this CH equation, we have : Implicit term (c n ) are on the left and the explicit term (c n−1 ) are on the right, a 1 and a 2 are fixed real numbers. Unconditionally algorithms are obtained for a 1 > 2 and a 2 < 1/2 [9,10,25]. Chen and Shen [5] proposed the use of a semi-implicit Fourier-spectral method for solving the time-dependent Ginzburg Landau and CH equations. The advantage of the method proposed in Chen's paper is that no assumption is made on the bulk free energy (double-well potential). With the semi-implicit scheme proposed, the time step can be increased by about 100 times as compared to the one used in an explicit scheme. However, it should be emphasized that this scheme is still only first-order accurate in time, and hence is not sufficient if one is interested in time-dependent solutions. We note that Zhu et al. [26] used a similar approach to treat a variable mobility case. They left however the nonlinear, potential term explicit in their time discretization which carries a consequent time step restriction.
If we applied the Fourier-spectral method to Eq. (12), this equation can be directly solved as in [6,7] where the k -dependent effective time step is : and ∂ĉ/∂t is the Fourier transform of ∂c/∂t in Eq. (11) and is a function ofĉ k (t). As Eq. (14) reveals, the unconditionally stable algorithm has a mode-dependent effective time step ∆t ef f (k, ∆t).

Algorithms for the causal-Cahn-Hilliard equation.
In this section, we describe algorithm to solve the causal Cahn-Hilliard equation. Previously, semi-implicit scheme was used by Lecoq et al. [22] to characterize spinodal decomposition with a time delay process. The splitting into contractive and expansive parts of the free energy in the evolution of the modified CH system Eq. (7) and Eq. (9) leads to decreasing the free energy with arbitrary time step as it has been proved in [21]. We extend here the splitting method presented in Sec.
(3) to the causal-Cahn-Hilliard equation Eq. (7). As previously, a mix of explicit and implicit terms in the update is introduced in Eq. (7).
where f 1c (c) and f 2c (c) are part of the derivative of Gibbs free energy given in Eq. (2). More precisely, we introduce the following splitting : can either be taken at the timet n ort n−1 and the non-linear term f 2c (c) is taken att n−1 .
The fourth order derivative in position can either be taken at the timet n ort n−1 . To formulate a semi-implicit expansion for Eq. (15), it is preferable to take this last term at the new timet n . Therefore, expanding the terms δc n and δc n−1 in Eq. (15), one gets Using the Fourier-spectral method for the spatial variables, and, treating the linear fourth-order operators implicitly and the non-linear terms explicitly, the first-order semiimplicit Fourier-spectral scheme is described by Hereĉ n k is the Fourier transform of c n , f 2c (c n−1 ) k is the Fourier transform of the non-linear termf 2c (c n−1 ),k = (k 1 , k 2 , k 3 ) is a vector in the Fourier space, andk = k 2 1 + k 2 2 + k 2 3 is the magnitude ofk. Collecting terms and multiplying left and right hand-sides of Eq. (17) byτ d ∆t, substituting f 1c (c) by −2Ωc, finally one gets :  13) proposed by [6,7]. Taking a 1 = 0 and a 2 = 1, the numerical scheme is equivalent to the one proposed by Lecoq et al. [22]. 5. Stability analysis. In order to prove the stability of the time discretization of equation (7), we will use the decomposition introduced by Eyre [9] for the standard Cahn-Hilliard equation. For clarity, the symbol˜is skipped. We assume that and where f 2c : R → R is a C 1 convex function with a derivative f 2c which is Lipschitz continuous on R, i.e. there exists a positive constant L such that The scheme reads : for c n−2 and c n−1 given, find c n which solves Eq. (16). In actual 3d-computations, this scheme is solved by a fast Fourier transform method, so that we consider periodic boundary conditions on the parallelepiped V 0 = Π d k=1 ]0, L i [ (here, d = 1, 2 or 3 is the dimension and L i > 0 for 1 ≤ i ≤ d is the length in the ith direction). Note that no-flux boundary conditions or homogeneous Dirichlet boundary conditions could also be considered in a similar way. Define the total mass of the system as in Eq. (10) : If c n−2 , c n−1 and c n are regular enough, and if M t (c n−1 ) = M t (c n−2 ), then on integrating equation (16) on V 0 and integrating by parts the right side, it is easily seen that M t (c n ) = M t (c n−1 ). By induction, if M t (c 1 ) = M t (c 0 ), then M t (c n ) = M t (c 0 ) for all n; this means that the total mass of the system is conserved as in the continuous case (10).
To be more specific, we introduce now the L 2 (V 0 )-Hilbert space with L 2 -scalar prod- and we denote V = H 1 per (V 0 ) the L 2 -Sobolev space on V 0 endowed with periodic boundary conditions; V is the space of functions with finite energy. We also writė Using the inclusionV ⊂L 2 (V 0 ) ⊂V whereV =Ḣ −1 per (V 0 ) is the dual space ofV , it can be seen (see [16] for details) that the operator −∆ is a linear isomorphism fromV ontoV defined by For all v ∈ V and for allv ∈V , we denote The free energy of the system is As a shortcut, we also denote ∂c n = c n − c n−1 .
The well-posedness and stability result reads : Theorem 5.1. Assume that assumptions (19) and (20) are satisfied, and that c 0 , c 1 are chosen in V such that M t (c 1 ) = M t (c 0 ). If a 2 > 0 and a 1 ≤ 0, or if a 2 > 0, a 1 > 0 and 4a 2 1 Ω 2 < a 2 2 c /∆t, then the scheme (16) defines a unique sequence (c n ) n≥0 in V . Moreover, M t (c n ) = M t (c 0 ) for all n, and for all n ≥ 1. In particular, if the scheme (16) is unconditionnally stable.
Proof. We first recall that for any real number s, the operator −∆ can be extended as a linear operator from H s per (V 0 ) into H s−2 per (V 0 ), where the spaces H s per (V 0 ) can be defined by spectral theory (see [16] for details). Moreover, by assumption (20), for every v ∈ L 2 (V 0 ), the function f 2c (v) belongs to L 2 (V 0 ). The scheme (16) can therefore be written Bc n = G(c n−1 , c n−2 ), where and where G(c n−1 , c n−2 ) belongs to the space H −3 per (V 0 ). If a 2 > 0 and a 1 ≤ 0, then . This proves the well-posedness of the scheme when a 2 > 0 and a 1 ≤ 0. If a 2 > 0 and a 1 > 0, then by the Cauchy-Schwarz inequality, Moreover, if 4a 2 1 Ω 2 < a 2 2 c /∆t, then we write 4a 2 1 Ω 2 = a 2 2 c /∆t for some constants ∆t > ∆t and a 2 < a 2 , and we have Thus, B defines an isomorphism form H 2 per (V 0 ) onto H −2 per (V 0 ) and we conclude as previously.
For the conservation of the total mass, the argument used above (assuming enough regularity) can also be extended to the case where c n−1 and c n−2 belong to V only. For the stability, we first write the scheme as τ D ∆t 2 (∂c n − ∂c n−1 ) + Next, we apply the operator (−∆) −1 to (23), we multiply by ∂c n ∈V , and we integrate on V 0 . We find τ D ∆t 2 (∂c n − ∂c n−1 , (−∆) −1 ∂c n ) + 1 ∆t |∂c n | 2 −1 − 2Ω(c n , c n − c n−1 ) +2(1 − a 1 )Ω|∂c n | 2 0 − 2 (∆c n , c n − c n−1 ) − (1 − a 2 ) 2 c |∂c n | 2 1 = −(f 2c (c n−1 ), ∂c n ) We use three times the identity (a, a − b) = 1 2 (|a| 2 − |b| 2 + |a − b| 2 ) (one time for every norm | · | 0 , | · | 1 and | · | −1 ), and we obtain In the right side, we have used the Taylor inequality which is a consequence of (20). Inequality (24) implies the stability estimate (21), and the proof is complete. In actual computations, the concentration is seen to stay in the interval [β, 1 − β] with β > 0 small. As in the polynomial case, we can modify f 2c outside of [β, 1 − β] so as to satisfy assumption (20), with Theorem 5.1 can be applied with these values (under the additionnal assumption that c n remains in [β, 1 − β] for all n). Example 3 : if we choose a 1 = 0 and a 2 = 1, we recover the scheme used in [22]. Under the assumptions of Theorem 5.1, if L = L − 2Ω ≥ 0 (as in Example 1), the scheme is not unconditionnally stable. However, using the interpolation inequality we deduce from (21) that The scheme is therefore stable if ∆t < 8 2 c /L 2 . Remark 1 The stability results states that the sum of the potential energy, E(c n ), and of the kinetic energy, τ D 2 ∂c n ∆t 2 −1 , is nonincreasing. (20) is essential in the proof of stability. It was introduced in [24] for the proof of stability of linearly implicit schemes for Allen-Cahn and Cahn-Hilliard equations. If one wants to work with the polynomial nonlinearity, then nonlinear implicit schemes should be considered (see for instance [15]). Remark 3 The analysis of the error estimate is more involved, but it can be said roughly that the scheme (16) has order O(∆t). We refer the reader to [24] for an error estimate of fully discretized linearly implicit schemes for the Cahn-Hilliard equation, and to [16] for error estimates of a space discretization of the hyperbolic Cahn-Hilliard equation.

Remark 2 Assumption
6. Numerical Stability Tests. Tests of the numerical scheme given by Eq. (16) were performed with several couples of (a 1 , a 2 ). For clarity, only results with (a 1 = 0.25, a 2 = 0.80) are presented below, but the findings are available for all the ones tested (a 1 ∈ [−1, 1] and a 2 ∈ [0, +1]). We compare the maximum size of time step that one can use in the scheme with various couples (a 1 , a 2 ). This maximum was estimated from the numerical computations as the maximum value that one can use still produces the same patterns and structure factor as those obtain using smaller size of time steps. The largest value ∆t = 0.25 was obtained with (a 1 = 0.25, a 2 = 0.80), which is about 200-300 times larger than the time step used in the no-splitting scheme. It should be noted that to avoid any problems in the implementation of the logarithmic nonlinearity Eq.
The splitting semi-implicit scheme with (a 1 = 0.0, a 2 = 1.0) was used previously [22] to numerically predict evolution of patterns, structure factor S(k, t) in delay spinodal decomposition. It has been found that the inclusion of memory into the model leads to a delay of the phase separation kinetics due to the appearance of inertia in the causal CH model. The features obtained with several pairs (a 1 , a 2 ) are similar and not presented here for clarity (see Fig. (1-3) of [22] for more details). We found that the mean of c is preserved within 3-4 digits (Eq. (10)) and this for all the values of τ D tested. Fig. 1 shows the time behaviour of the potential energy. This energy decreases monotically and smoothly for the Cahn-Hilliard case throughout the entire computation as required by Eq. (10). For the causal CH case, the potential energy decreases with oscillations, and as the time increases, energy tends towards the one obtained with the CH case. These oscillations depend only on the value of τ D and are linked to the ones observed on the structure factor S(k) where a wave behaviour was evaluated (see Fig. (3) of reference [22]) and not on numerical instabilities.