# Numerical approximation of the Cahn–Hilliard equation with memory effects in the dynamics of phase separation

N. Lecoq, Helena Zapolsky, Peter Galenko

► To cite this version:

N. Lecoq, Helena Zapolsky, Peter Galenko. Numerical approximation of the Cahn–Hilliard equation with memory effects in the dynamics of phase separation. AIMS Proceedings, 2011, 31, pp.953-962. 10.3934/proc.2011.2011.953 . hal-01777456

HAL Id: hal-01777456

https://normandie-univ.hal.science/hal-01777456v1

Submitted on 20 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
NUMERICAL APPROXIMATION OF THE CAHN-HILLIARD
EQUATION WITH MEMORY EFFECTS IN THE DYNAMICS OF
PHASE SEPARATION

NICOLAS LECOQ AND HELENA ZAPOLSKY

Groupe de Physique des Matériaux, UMR CNRS 6634
Université de Rouen
76801 Saint-Étienne du Rouvray, France

PETER GALENKO

German Aerospace, Institut für Materialphysik im Weltraum
Deutsches Zentrum für Luft- und Raumfahrt (DLR)
51170 Köln, Germany

ABSTRACT. We consider the modified 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 coefficient $\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$

$$ \tau_D \frac{\partial^2 c}{\partial t^2} + \frac{\partial c}{\partial t} = \nabla \cdot [M \nabla (f'_c - \epsilon_c^2 \nabla^2 c)], \quad (1) $$

as proposed by Galenko and co-workers [11, 12]. Here $\tau_D > 0$ is a relaxation time, and the term $\tau_D \partial^2 c / \partial 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.

Recently, the causal CH model Eq. (1) has been analysed mathematically [13, 16] and more recently, the analysis was extended to the three-dimension Cahn-Hilliard equation including the memory term [14]. Most of these study introduce a nonconvex potential which has a double-well structure where the two mimina correspond to the two equilibrium phases. In the present study, the potential is assumed to be a regular solution [23] the Gibbs free energy has the form :

$$ f(T, c) = \Omega(T)c(1-c) + RT \{c \log(c) + (1-c) \log(1-c)\} \quad (2) $$

where the interaction parameter $\Omega$ is a function of temperature $T$ and $R$ is the gas constant. The use of logarithmic terms in the free energy, which arises naturally due

<footnote n="2000">2000 Mathematics Subject Classification. Primary: 65M12, 65M22, 68U20; Secondary: 74D10.</footnote>

<footnote>Key words and phrases. Cahn-Hilliard equation, Spinodal decomposition, Memory effects, Numerical analysis.</footnote>
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 developing 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, \vec{r})$ having meaning of concentration of a second chemical component in the binary system, one may write the flux

$$ \vec{j} = - \int_{-\infty}^{t} M_R(t-t^*) \nabla \left( \frac{\delta E(t^*, \vec{r})}{\delta c} \right) dt^*, \quad (3) $$

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 $\partial c/\partial t = -\nabla \cdot \vec{j}$, Eq. (3) gives the following equation for the arbitrary, slow or fast, dynamics:

$$ \frac{\partial c}{\partial t} = \nabla \cdot \int_{-\infty}^{t} M_R(t-t^*) \nabla \left( \frac{\delta E(t^*, \vec{r})}{\delta c} \right) dt^*. \quad (4) $$

Taking the free energy $E$ in the following functional form

$$ E[c(t^*, \vec{r})] = \int_{V_0} \left( f(T, c) + \frac{\epsilon_c^2}{2} |\nabla c|^2 \right) dv_0, \quad (5) $$

and the memory function of the exponential form $M_R(t-t^*) = (D/\tau_D) \exp[-(t-t^*)/\tau_D]$, Eq. (4) can be rewritten as [12]

$$ \tau_D \frac{\partial^2 c}{\partial t^2} + \frac{\partial c}{\partial t} = \nabla \cdot [M \nabla (f'_c - \epsilon_c^2 \nabla^2 c)], \quad (6) $$

where $v_0$ is the subvolume of the system, $\tau_D$ the time for relaxation of the diffusion flux $\vec{j}$ to its steady-state, $M$ the atomic mobility, $\epsilon_c^2$ 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 $\tau_D \to 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 $\tau_D \partial^2 c/\partial 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: $\tilde{t} = tM\Delta f/l_D^2$ dimensionless time, $\tilde{\nabla} = \nabla/l_D$ normalized operator “nabla”, $\tilde{f}'_c \equiv \partial\tilde{f}/\partial c = f'_c/\Delta f$ dimensionless characteristic for the free energy
density, $\tilde{\epsilon}_c^2 = \epsilon_c^2 l_D^2 / \Delta f$ dimensionless length scale of the interface, and $\tilde{\tau}_D = \tau_D M \Delta f / l_D^2$ dimensionless relaxation time. Then, Eq. (6) to be numerically solved takes the form :

$$ \frac{\partial^2 c}{\partial \tilde{t}^2} + \frac{1}{\tilde{\tau}_D} \frac{\partial c}{\partial \tilde{t}} = \frac{1}{\tilde{\tau}_D} \left[ \tilde{\nabla}^2 \tilde{f}'_c - \tilde{\epsilon}_c^2 \tilde{\nabla}^4 c \right]. \quad (7) $$

A natural boundary condition, following from the variation of functional (5), is given by

$$ \epsilon_c^2 \frac{\partial c}{\partial t} \nabla_n c = (f'_c - \epsilon_c^2 \nabla^2 c) j_n, \quad (8) $$

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

$$ \nabla \cdot j_n = 0 \quad \text{and} \quad \nabla c = 0. \quad (9) $$

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 :

$$ \frac{d}{dt} M_t = \frac{d}{dt} \int_{V_0} c(r, t) dv_0 = 0 \quad \frac{d}{dt} E[c(t^*, \vec{r})] \le 0 \quad (10) $$

### 3. Unconditionally stable algorithm for the Cahn-Hilliard equation.
For the Cahn-Hilliard (CH) equation, with a free energy functional described with a double-well structure, and a gradient-squared term associated with an energy cost to an interface between two phases, several unconditionally stable algorithms were developed [9, 10, 25]. Considering the simple case where the double-well potential is represented in a polynomial form, the CH equation takes the form :

$$ \frac{\partial c}{\partial t} = -\nabla^2(\nabla^2 c + c - c^3) \quad (11) $$

A mix of explicit and implicit terms in the update is introduced. For this CH equation, we have :

$$ \begin{aligned} c^n - (a_1 - 1)\Delta t \nabla^2 c^n - (a_2 - 1)\Delta t \nabla^4 c^n \\ = c^{n-1} - \Delta t(a_1 \nabla^2 c^{n-1} + a_2 \nabla^4 c^{n-1} - \nabla^2(c^n)^3) \end{aligned} \quad (12) $$

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]:

$$ \hat{c}_k^n = \hat{c}_k^{n-1} + \Delta t_{eff}(k, \Delta t) \frac{\partial \hat{c}_k^n}{\partial t} \quad (13) $$

where the k-dependent effective time step is:

$$ \Delta t_{eff}(k, \Delta t) \equiv \frac{\Delta t}{1 + \Delta t k^2 [(a_1 - 1) + (1 - a_2)k^2]} \quad (14) $$

and $\partial \hat{c}/\partial t$ is the Fourier transform of $\partial c/\partial t$ in Eq. (11) and is a function of $\hat{c}_k(t)$. As Eq. (14) reveals, the unconditionally stable algorithm has a mode-dependent effective time step $\Delta t_{eff}(k, \Delta t)$.

### 4. 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).

$$ \frac{\tilde{\tau}_D}{\Delta t} (\delta c^n - \delta c^{n-1}) + \delta c^{n-1} = \tilde{\nabla}^2 (a_1 \tilde{f}'_{1c}(c^n) + (1-a_1) \tilde{f}'_{1c}(c^{n-1}) + \tilde{f}'_{2c}(c^{n-1}) \\ -a_2 \tilde{\epsilon}_c^2 \tilde{\nabla}^2 c^n - (1-a_2) \tilde{\epsilon}_c^2 \tilde{\nabla}^2 c^{n-1}) \quad (15) $$

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:

$$ f'_{1c}(c) = -2\tilde{\Omega}c \quad f'_{2c}(c) = \log(c) - \log(1-c) + \tilde{\Omega} $$

with $\tilde{\Omega} = \Omega/RT$ so that $f'_{1c}(c)$ can either be taken at the time $\tilde{t}_n$ or $\tilde{t}_{n-1}$ and the non-linear term $f'_{2c}(c)$ is taken at $\tilde{t}_{n-1}$.

The fourth order derivative in position can either be taken at the time $\tilde{t}_n$ or $\tilde{t}_{n-1}$. To formulate a semi-implicit expansion for Eq. (15), it is preferable to take this last term at the new time $\tilde{t}_n$. Therefore, expanding the terms $\delta c^n$ and $\delta c^{n-1}$ in Eq. (15), one gets

$$ \frac{\tilde{\tau}_D}{\Delta t} \left( \frac{c^n - c^{n-1}}{\Delta t} - \frac{c^{n-1} - c^{n-2}}{\Delta t} \right) + \frac{c^n - c^{n-1}}{\Delta t} = \tilde{\nabla}^2 \left( a_1 \tilde{f}'_{1c}(c^n) + (1-a_1) \tilde{f}'_{1c}(c^{n-1}) \right. \\ \left. + \tilde{f}'_{2c}(c^{n-1}) - a_2 \tilde{\epsilon}_c^2 \tilde{\nabla}^2 c^n - (1-a_2) \tilde{\epsilon}_c^2 \tilde{\nabla}^2 c^{n-1} \right) \quad (16) $$
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 semi- implicit Fourier-spectral scheme is described by

$$
\begin{equation}
\begin{split}
& \frac{1}{\Delta t} \left( \frac{\hat{c}_k^n - \hat{c}_k^{n-1}}{\Delta t} - \frac{\hat{c}_k^{n-1} - \hat{c}_k^{n-2}}{\Delta t} \right) + \tilde{\tau}_D^{-1} \frac{\hat{c}_k^n - \hat{c}_k^{n-1}}{\Delta t} = \\
& \qquad - \hat{k}^2 \tilde{\tau}_D^{-1} \left( a_1 \tilde{f}_{1c}^{\prime}(\hat{c}_k^n) + (1-a_1) \tilde{f}_{1c}^{\prime}(\hat{c}_k^{n-1}) + \widehat{(\tilde{f}_{2c}^{\prime}(c^{n-1}))}_k \right) \\
& \qquad + a_2 \hat{k}^2 \tilde{\epsilon}_c^2 \hat{c}^n + (1-a_2) \hat{k}^2 \tilde{\epsilon}_c^2 \hat{c}^{n-1}
\end{split}
\tag{17}
\end{equation}
$$

Here $\hat{c}_k^n$ is the Fourier transform of $c^n$, $\widehat{(\tilde{f}_{2c}^{\prime}(c^{n-1}))}_k$ is the Fourier transform of the non-linear term $\tilde{f}_{2c}^{\prime}(c^{n-1})$, $\hat{k} = (k_1, k_2, k_3)$ is a vector in the Fourier space, and $\hat{k} = \sqrt{k_1^2 + k_2^2 + k_3^2}$ is the magnitude of $\hat{k}$. Collecting terms and multiplying left and right hand-sides of Eq. (17) by $\tilde{\tau}_d \Delta t$, substituting $f_{1c}^{\prime}(c)$ by $-2\tilde{\Omega}c$, finally one gets :

$$
\begin{equation}
\begin{split}
\hat{c}_k^n [\tilde{\tau}_D \Delta t^{-1} + 1 + 2\hat{k}^2 a_1 \tilde{\Omega} + a_2 \hat{k}^4 \Delta t \hat{\epsilon}_c^2] = \\
\hat{c}_k^{n-1} [2\tilde{\tau}_D \Delta t^{-1} - 1 + \hat{k}^2 \Delta t ((1-a_1)\tilde{\Omega} + (1-a_2)\hat{k}^2 \hat{\epsilon}_c^2)] \\
- \tilde{\tau}_D \Delta t^{-1} \hat{c}_k^{n-2} - \hat{k}^2 \Delta t (\widehat{\tilde{f}_{2c}^{\prime}(c^{n-1})})_k
\end{split}
\tag{18}
\end{equation}
$$

Eq. (18) represents the splitting semi-implicit numeric scheme for solution of the hyperbolic equation. It allows to obtain the concentration $\hat{c}_k^n$ at the current time step $n$ by previously established values $\hat{c}_k^{n-1}$ and $\hat{c}_k^{n-2}$ at the time steps $n-1$ and $n-2$, respectively. Taking $\tau_D = 0$, Eq. (18) is equivalent to the formulation of Eq. (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 $\tilde{\cdot}$ is skipped. We assume that

$$
f_c = f_{1c}(c) + f_{2c}(c), \quad \text{with } f_{1c}(c) = -\Omega c^2 \quad (\Omega > 0), \tag{19}
$$

and where $f_{2c}: \mathbb{R} \to \mathbb{R}$ is a $C^1$ convex function with a derivative $f_{2c}'$ which is Lipschitz continuous on $\mathbb{R}$, i.e. there exists a positive constant $L$ such that

$$
|f'_{2c}(c) - f'_{2c}(\hat{c})| \le L |c - \hat{c}|, \quad \forall c, \hat{c} \in \mathbb{R}. \tag{20}
$$

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 = \Pi_{k=1}^d]0, L_i[$ (here, $d=1, 2$ or $3$ is the dimension and $L_i > 0$ for $1 \le i \le d$ is the length in the $i$th 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) :

$$
M_t(c) = \int_{V_0} c \, dv_0
$$
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 product

$$ (u, v) = \int_{V_0} uv \, dv_0 \quad \forall u, v \in L^2(V_0), $$

and we denote $V = H_{per}^1(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 write

$$ \dot{L}^2(V_0) = \{v \in L^2(V_0) : M(v) = 0\} \quad \text{and} \quad \dot{V} = \{v \in V : M(v) = 0\}. $$

Using the inclusion $\dot{V} \subset \dot{L}^2(V_0) \subset \dot{V}'$ where $\dot{V}' = \dot{H}_{per}^{-1}(V_0)$ is the dual space of $\dot{V}$, it can be seen (see [16] for details) that the operator $-\Delta$ is a linear isomorphism from $\dot{V}$ onto $\dot{V}'$ defined by

$$ (-\Delta u, v) = (\nabla u, \nabla v) \quad \forall u, v \in \dot{V}. $$

For all $v \in V$ and for all $\dot{v} \in \dot{V}$, we denote

$$ |v|_0 = (v, v)^{1/2}, \quad |v|_1 = (\nabla v, \nabla v)^{1/2} \quad \text{and} \quad |\dot{v}|_{-1} = |(-\Delta)^{-1} \dot{v}|_1. $$

The free energy of the system is

$$ E(c) = \int_{V_0} f_c(c) + \frac{\epsilon_c^2}{2} |\nabla c|^2 dv = (f_{1c}(c), 1) + (f_{2c}(c), 1) + \frac{\epsilon_c^2}{2} |c|^2_1 $$

As a shortcut, we also denote $\partial 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 \le 0$, or if $a_2 > 0$, $a_1 > 0$ and $4a_1^2\Omega^2 < a_2\epsilon_c^2/\Delta t$, then the scheme (16) defines a unique sequence $(c^n)_{n \ge 0}$ in $V$. Moreover, $M_t(c^n) = M_t(c^0)$ for all $n$, and

$$ \begin{aligned} E(c^n) + \frac{\tau_D}{2\Delta t^2} |\partial c^n|_{-1}^2 + \epsilon_c^2(a_2 - \frac{1}{2}) |\partial c^n|_1^2 &+ \left[ 2\Omega(\frac{1}{2} - a_1) - \frac{L}{2} \right] \\ |\partial c^n|_0^2 + \frac{1}{\Delta t} |\partial c^n|_{-1}^2 &\le E(c^{n-1}) + \frac{\tau_D}{2\Delta t^2} |\partial c^{n-1}|^2. \end{aligned} \quad (21) $$

for all $n \ge 1$. In particular, if

$$ a_2 \ge \frac{1}{2} \quad \text{and} \quad 2\Omega(\frac{1}{2} - a_1) - \frac{L}{2} \ge 0, \quad (22) $$

the scheme (16) is unconditionally stable.

*Proof.* We first recall that for any real number $s$, the operator $-\Delta$ can be extended as a linear operator from $H_{per}^s(V_0)$ into $H_{per}^{s-2}(V_0)$, where the spaces $H_{per}^s(V_0)$ can be defined by spectral theory (see [16] for details). Moreover, by assumption (20), for every $v \in L^2(V_0)$, the function $f'_{2c}(v)$ belongs to $L^2(V_0)$. The scheme (16) can therefore be written

$$ Bc^n = \mathcal{G}(c^{n-1}, c^{n-2}), $$

where

$$ Bc^n = \left( \frac{\tau_D}{\Delta t^2} + \frac{1}{\Delta t} \right) c^n + 2a_1 \Omega \Delta c^n + a_2 \epsilon_c^2 \Delta^2 c^n, $$
and where $\mathcal{G}(c^{n-1}, c^{n-2})$ belongs to the space $H_{per}^{-3}(V_0)$.
If $a_2 > 0$ and $a_1 \le 0$, then

$$ (Bc, c) \geq \min(\Delta t^{-1}, a_2 \epsilon_c^2)(|c_0^2| + |\Delta c_0^2|) \quad \text{for all } c \in H_{per}^2(V_0), $$

so that $B$ is a linear isomorphism from $H_{per}^2(V_0)$ onto $H_{per}^{-2}(V_0)$. By extension, $B$ maps isomorphically $H_{per}^s(V_0)$ onto $H_{per}^{s-4}(V_0)$ for all $s \in \mathbb{R}$, so that $c^n = B^{-1}\mathcal{G}(c^{n-1}, c^{n-2})$ is well-defined in $H_{per}^1(V_0)$. This proves the well-posedness of the scheme when $a_2 > 0$ and $a_1 \le 0$.

If $a_2 > 0$ and $a_1 > 0$, then by the Cauchy-Schwarz inequality,

$$ (Bc, c) \geq \frac{1}{\Delta t} |c_0^2 - 2|a_1|\Omega|\Delta c_0|c_0 + a_2 \epsilon_c^2 |\Delta c_0|^2. $$

Moreover, if $4a_1^2\Omega^2 < a_2\epsilon_c^2/\Delta t$, then we write $4a_1^2\Omega^2 = a_2'\epsilon_c^2/\Delta t'$ for some constants $\Delta t' > \Delta t$ and $a_2' < a_2$, and we have

$$ (Bc, c) \geq \left( \frac{1}{\Delta t} - \frac{1}{\Delta t'} \right) |c_0^2 + (a_2 - a_2')\epsilon_c^2 |\Delta c_0|^2. $$

Thus, $B$ defines an isomorphism from $H_{per}^2(V_0)$ onto $H_{per}^{-2}(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

$$ \frac{\tau_D}{\Delta t^2}(\partial c^n - \partial c^{n-1}) + \frac{\partial c^n}{\Delta t} - \Delta f_{1c}'(c^n) + (1-a_1)\Delta[f_{1c}'(c^n) - f_{1c}'(c^{n-1})] \\ -(1-a_2)\epsilon_c^2\Delta^2\partial c^n + \epsilon_c^2\Delta^2 c^n = \Delta f_{2c}'(c^{n-1}), \quad (23) $$

Next, we apply the operator $(-\Delta)^{-1}$ to (23), we multiply by $\partial c^n \in \dot{V}$, and we integrate on $V_0$. We find

$$ \frac{\tau_D}{\Delta t^2} (\partial c^n - \partial c^{n-1}, (-\Delta)^{-1} \partial c^n) + \frac{1}{\Delta t} |\partial c^n|_{-1}^2 - 2\Omega(c^n, c^n - c^{n-1}) \\ +2(1-a_1)\Omega|\partial c^n|_0^2 - \epsilon^2(\Delta c^n, c^n - c^{n-1}) - (1-a_2)\epsilon_c^2|\partial c^n|_1^2 = -(f_{2c}'(c^{n-1}), \partial c^n) $$

We use three times the identity

$$ (a, a-b) = \frac{1}{2}(|a|^2 - |b|^2 + |a-b|^2) $$

(one time for every norm $|\cdot|_0, |\cdot|_1$ and $|\cdot|_{-1}$), and we obtain

$$ \frac{\tau_D}{2\Delta t^2} \left[ |\partial c^n|_{-1}^2 - |\partial c^{n-1}|_{-1}^2 + |\partial c^n - \partial c^{n-1}|_{-1}^2 \right] + \frac{1}{\Delta t} |\partial c^n|_{-1}^2 \\ - \Omega \left( |c^n|_0^2 - |c^{n-1}|_0^2 + |c^n - c^{n-1}|_0^2 \right) + 2(1 - a_1)\Omega |\partial c^n|_0^2 \\ + \frac{\epsilon_c^2}{2} \left[ |c^n|_1^2 - |c^{n-1}|_1^2 + |c^n - c^{n-1}|_1^2 \right] - (1 - a_2)\epsilon_c^2 |\partial c^n|_1^2 \\ \le (f_{2c}(c^{n-1}), 1) - (f_{2c}(c^n), 1) + \frac{L}{2} |c^n - c^{n-1}|_0^2. \quad (24) $$

In the right side, we have used the Taylor inequality

$$ \left| f_{2c}(c^n) - f_{2c}(c^{n-1}) - f_{2c}'(c^{n-1})(c^n - c^{n-1}) \right| \le \frac{L}{2} |c^n - c^{n-1}|^2, $$

which is a consequence of (20). Inequality (24) implies the stability estimate (21), and the proof is complete. $\square$
**Example 1:** a standard choice for $f_c$ is $f_c(c) = (c^2 - 1)^2/4$. The renormalized concentration $c$ is expected to belong to $[-1, 1]$ [10]. We have $f'_c(c) = c^3 - c$ and we can choose $f'_{1c}(c) = -c$ and $f'_{2c}(c) = c^3$ on $[-1, 1]$ and $f''_{2c}(c)$ constant if $|c| \ge 1$. With these choices, we have $\Omega = 1/2$ and $\text{decoq}$, $L = \sup_{c \in [-1,1]} f''_{2c}(c) = 3$; by Theorem 5.1, the scheme (16) is unconditionally stable if

$$a_2 \ge 1/2 \quad \text{and} \quad a_1 \le -1.$$

This result generalizes the values found by Chen and Shen in [5] for the limit case $\tau_D = 0$.

**Example 2:** we consider the standard logarithmic potential $f_c = f_{1c} + f_{2c}$ with

$$f'_{1c}(c) = -2\Omega c \quad \text{and} \quad f'_{2c}(c) = \log c - \log(1-c) + \Omega.$$

In actual computations, the concentration is seen to stay in the interval $[\beta, 1 - \beta]$ with $\beta > 0$ small. As in the polynomial case, we can modify $f'_{2c}$ outside of $[\beta, 1 - \beta]$ so as to satisfy assumption (20), with

$$L = \sup_{c \in [\beta, 1-\beta]} f''_{2c}(c) = \frac{1}{\beta(1-\beta)}.$$

Theorem 5.1 can be applied with these values (under the additional assumption that $c^n$ remains in $[\beta, 1 - \beta]$ 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\Omega \ge 0$ (as in Example 1), the scheme is not unconditionally stable. However, using the interpolation inequality

$$\frac{L'}{2} |\partial c^n|_0^2 \le \frac{L'}{2} |\partial c^n|_1 |\partial c^n|_{-1} \le \frac{\epsilon_c^2}{2} |\partial c^n|_1^2 + \frac{L'^2}{8\epsilon_c^2} |\partial c^n|_{-1}^2,$$

we deduce from (21) that

$$E(c^n) + \frac{\tau_D}{2\Delta t^2} |\partial c^n|_{-1}^2 + \left( \frac{1}{\Delta t} - \frac{L'^2}{8\epsilon_c^2} \right) |\partial c^n|_{-1}^2 \le E(c^{n-1}) + \frac{\tau_D}{2\Delta t^2} |\partial c^{n-1}|^2.$$

The scheme is therefore stable if $\Delta t < 8\epsilon_c^2/L'^2$.

**Remark 1** The stability results states that the sum of the potential energy, $E(c^n)$, and of the kinetic energy, $\frac{\tau_D}{2} \left| \frac{\partial c^n}{\partial t} \right|_{-1}^2$, is nonincreasing.

**Remark 2** Assumption (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(\Delta 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.

**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 \in [-1, 1]$ and $a_2 \in [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 $\Delta 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. (2), when the concentration $c$ reaches a value close to 0 (respectively 1), we truncate it to $10^{-6}$ (resp. $(1 - 10^{-6})$). This constraint is in agreement with the assumptions of Theorem 5.1 and the Example 2.

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 $\tau_D$ tested. Fig. 1 shows the time behaviour of the potential energy. This energy decreases monotonically 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 $\tau_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.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Behaviour of the potential energy $E[c(t, r)]$ in time for the CH and Causal CH equations. Time corresponds to time steps in units of $\tau_D$; ($a_1, a_2$) = (0.25, 0.80), $\Delta t = 0.25$.</figcaption></figure>

**Acknowledgments.** The authors acknowledge M. Pierre (UMR CNRS 6068, University of Poitiers, France) for valuable discussions. N. Lecoq and H. Zapolsky acknowledge support by the computer-time Grants 2006007 of the Centre de Ressources Informatiques de Haute-Normandie (CRIHAN). They acknowledge too the PEPS-CNRS Maths-ST2I CPCRISTAL "Champ de Phase CRISTalins. Transformation de phase dans les alliages ALuminium" for financial support. P. Galenko acknowledges financial support from the German Research Foundation (DFG) under the Project No. HE 160/19 and DLR Agency under contract 50WM1140.
REFERENCES

[1] V.E. Badalassi, H.D. Ceniceros and S. Banerjee *Computation of multiphase systems with phase field models*, J. Comput. Phys. **190**, (2003) 371–397.

[2] J.W. Cahn, J.E. Hilliard, *Free energy of a nonuniform system. I. Interfacial free energy* J. Chem. Phys. **28**, (1958) 258–267.

[3] J.W. Cahn, *On Spinodal Decomposition* Acta Metall. **9**, (1961) 795–801.

[4] J.W. Cahn, Trans. Metall. Soc. AIME, **242**: 166 (1968).

[5] L.Q. Chen and J. Shen, *Applications of semi-implicit Fourier-spectral method to phase field equations*, Computer Phys. Com. **108**, (1998) 147–158.

[6] M. Cheng and A. Rutemberg, *Maximally fast coarsening algorithms*, Phys. Rev. E **72**, 055701(R) (2005).

[7] M. Cheng and J. Warren, *Controlling the accuracy of unconditionally stable algorithms in the Cahn-Hilliard equation*, Phys. Rev. E **75**, 017702 (2007).

[8] R.L.H. Essery, R.C. Ball, Europhys. Lett. **16**, (1991) 379.

[9] D.J. Eyre in *Computational and Mathematical Models of Microstructural Evolution* edited by J.W Bullard et al (The Material Research Society, Warrendale, PA, 1998) 39–46.

[10] D.J. Eyre, http://www.math.utah.edu/~eyre/research/methods/stables.ps

[11] P. Galenko, D. Jou, *Diffuse-interface Model for Rapid Phase Transformations in Non-equilibrium Systems* Phys. Rev. E **71**, (2005) 046125–1–13.

[12] P. Galenko, *Phase-Field Model with Relaxation of Diffusion Flux in Rapid Solidification*, Phys. Lett. A **287**, (2001) 190–198.

[13] S. Gatti, M. Grasselli, A. Miranville and V. Pata, *On the huperbolic relaxation of the one-dimensional Cahn-Hilliard equation*, J. Math. Anal. Appli. **312** (2005), 230–247.

[14] M. Grasselli, G. Schimperna, A. Segatti and S. Zelik, *On the 3D Cahn-Hilliard equation with inertial term*, J. Evol. Equ. **9** (2009), 371–404.

[15] M. Grasselli, N. Lecoq and M. Pierre, *A long-time stable fully discrete approximation of the Cahn-Hilliard equation with inertial term*, Discrete Contin. Dyn. Syst., to appear.

[16] M. Grasselli and M. Pierre, *A splitting method for the Cahn-Hilliard equation with inertial term*, Math. Models Methods Appl. Sci. **20** (2010), 1–28

[17] M. Hillert, *A theory of nucleation for solid metallic solutions*, D.Sc. thesis, MIT, 1956.

[18] M. Hillert, *A solid-solution model for inhomogeneous system*, Acta Metall. **9**, (1961)525–539.

[19] D. Jou, J. Casas-Vazquez, M. Criado-Sancho, *Thermodynamics of fluids under flow* (Springer, Berlin 2000).

[20] A.G. Khachaturyan, *Theory of structural transformations in solids* (Wiley, New York 1983).

[21] V. Lebedev, A. Sysoeva and P. Galenko, *Unconditionally gradient-stable computational schemes in problems of fast phase transitions* Phys Rev. E **83**, 026705 (2011)

[22] N. Lecoq, H. Zapolsky and P. Galenko, *Evolution of the structure factor in a hyperbolic model of spinodal decomposition* Eur. Phys J. Special Topics **177**, (2009) 165–175.

[23] D.A. Porter and K.E. Easterling, *Phase Transformations in Metals and Alloys*, Van Nostrand Reinhold (International), London, 1981.

[24] J. Shen and X. Yang, *Numerical approximations of Allen-Cahn and Cahn-Hilliard equations*, Discrete Contin. Dyn. Syst., **28** (2010), n. 4, 1669–1691.

[25] B.P. Vollmayr-Lee and A.D. Rutenberg, *Fast and accurate coarsening simulation with an unconditionally stable time step*. Phys. Rev. E **75**, 66703 (2007)

[26] J. Zhu, L.Q. Chen, J. Shen and V. Tikare, *Coarsening kinetics from a variable-mobility Cahn-Hilliard equation: Application of a semi-implicit Fourier spectral method*, Phys. Rev E, **60**, (1999), 3564–3572.

Received July 2010; revised May 2011.

E-mail address: nicolas.lecoq@univ-rouen.fr

E-mail address: helena.zapolsky@univ-rouen.fr

E-mail address: peter.galenko@dlr.de