2.1. The Cahn–Hilliard model

Let $\Omega$ be a bounded domain in $\mathbb{R}^d$ , $d\in \{1,2,3\}$ , containing two chemical species. We introduce a normalised difference ${\varphi }$ of the concentrations of two chemical components that is governed by the following Cahn–Hilliard equation with chemical reactions [Reference Zwicker, Seyboldt, Weber, Hyman and Jülicher38]:

(2.1a) \begin{align} \partial _t {{\varphi }} &= \textrm {div}(m({{\varphi }})\nabla \mu ) + S_\varepsilon ({{\varphi }}) \qquad && \text{in $Q \;:\!=\; (0,T) \times \Omega $,} \end{align}

(2.1b) \begin{align} \mu &= \beta \bigg(\!-\!\varepsilon \Delta {{\varphi }} + \frac 1 \varepsilon \psi '({{\varphi }})\bigg) \qquad && \text{in $Q$,} \end{align}

(2.1c) \begin{align} \partial _{{\boldsymbol{n}}} \mu &= \partial _{{\boldsymbol{n}}} {{\varphi }} =0 \qquad && \text{on $\Gamma \;:\!=\; (0,T) \times \partial \Omega $,} \end{align}

(2.1d) \begin{align} {{\varphi }}(0) &={{\varphi }}_0 \qquad && \text{in $\Omega $.} \end{align}

Here, $\mu$ is the associated chemical potential, $m{\;:\; {\mathbb{R}} \to {\mathbb{R}}_{\gt 0}}$ is the concentration dependent, strictly positive mobility function, $\beta \gt 0$ is a parameter related to a surface energy density, $\varepsilon \gt 0$ is a small length scale proportional to the thickness of the diffuse interface function, $\psi {\;:\; {\mathbb{R}} \to {\mathbb{R}}_{\geq 0}}$ is a double-well potential, $ \partial _{{\boldsymbol{n}}}$ is the derivative in the direction of the unit outer normal $\boldsymbol{n}$ to $\partial \Omega$ and ${{\varphi }}_0$ serves as initial data for ${\varphi }$ . We consider the above equations on the space-time cylinder $Q$ with a fixed but arbitrary time $T\gt 0$ .

The source term $S_\varepsilon \;:\; {\mathbb{R}} \to {\mathbb{R}}$ is given as

(2.2) \begin{equation} S_\varepsilon (r)= S_1(r)+\frac 1\varepsilon S_2(r){,\quad r \in {\mathbb{R}}}. \end{equation}

Here, the term $S_2$ will later lead to a fast reaction in the interfacial region. On choosing $r_c \in (0,1]$ , we set for constants $S_+$ , $S_-$ , $K_+$ , $K_-$ , and $L$

(2.3) \begin{align} S_1(r)= \begin{cases} {S_+} & \quad \text{if $r \geq {r_c}$,} \\ {S_-} + G_1(r)({S_+}-{S_-}) & \quad \text{if $r\in (\!-\!{r_c},{r_c})$,} \\ {S_-} & \quad \text{if $r \leq -{r_c}$,} \end{cases} \end{align}

and

(2.4) \begin{align} S_2(r)= \begin{cases} - {K_+} (r-1) & \quad \text{if $r \geq {r_c}$,} \\ \widehat {S}_2(r) & \quad \text{if $r\in (\!-\!{r_c},{r_c})$,} \\ - {K_-} (r+1)& \quad \text{if $r \leq -{r_c}$,} \end{cases} \end{align}

where we define for $r\in (\!-\!{r_c},{r_c})$

\begin{equation*} \widehat {S}_2(r)= - {{K_-}} G_2(r) - {{K_+}} G_3(r)+ L G_4(r) {- K_+(r_c-1)G_1(r) - K_-(1 - r_c)(1-G_1(r))}. \end{equation*}

Here, $G_1,G_2,G_3, G_4 \;:\; [-1,1] \to {\mathbb{R}}$ are suitable differentiable interpolation functions, to be introduced below, satisfying

(2.5) \begin{align} & G_1({r_c})=1, \quad G_1(\!-\!{r_c}) =0, \quad G_2(\pm {r_c}) =G_3(\pm {r_c}) =G_4(\pm {r_c}) =0, \end{align}

(2.6) \begin{align} & G_1'(\pm {r_c})=G_4'(\pm {r_c})=0, \quad G_2'({r_c})=G_3'(\!-\!{r_c})=0, \quad G_2'(\!-\!{r_c})=G_3'({r_c})=1, \end{align}

so that the source term $S_\varepsilon$ is differentiable on $\mathbb{R}$ . We often use $r_c=1$ , which considerably simplifies the expression for $\widehat {S}_2$ and the matched asymptotic expansions, which we use later to derive a sharp interface limit. However, other choices can also be considered, such as $r_c = \frac 12$ as in [Reference Zwicker, Seyboldt, Weber, Hyman and Jülicher38].

The system (2.1) is related to the following free energy

\begin{align*} \mathcal{E}({{\varphi }}) = \beta \Big ( \frac \varepsilon 2 {\int _\Omega } |\nabla {{\varphi }}|^2 + \frac 1\varepsilon {\int _\Omega } \psi ({{\varphi }}) \Big ). \end{align*}

In what follows, we assume that $\psi$ is even, that is $\psi (r)=\psi (\!-\!r)$ for $r \in {\mathbb{R}}$ , and satisfies $\psi (\pm 1) = 0$ and $\psi ''(\pm 1) \neq 0$ . Typically, we will choose the quartic potential

(2.7) \begin{align} \psi (r) = \frac 14 (1-r^2)^2, \quad r \in {\mathbb{R}}. \end{align}

2.2. Possible choice of the interpolation functions

As mentioned above, in the theoretical analysis to follow, we consider any interpolation functions $G_1$ , $G_2$ , $G_3$ and $G_4$ such that (2.5)–(2.6) are fulfilled. Here, we present a possible choice related to the double-well potential $\psi$ . We set, for every $r \in [-1,1]$ ,

(2.8) \begin{align} G_1(r) { = \widehat {G}_1\left(\frac {r}{r_c}\right), \quad \widehat {G}_1(r) } & = \frac 34 (r+1)^2- \frac 14(r+1)^3, \end{align}

(2.9) \begin{align} G_2(r) { = r_c \widehat {G}_2\left(\frac {r}{r_c}\right), \quad \widehat {G}_2(r) } & = -{\frac 12}\frac {1}{ \sqrt {\psi ''(\!-\!1)}} (r-1) \sqrt {2 \psi (r)}, \quad G_3 (r) = - G_2(\!-\!r), \end{align}

and observe that

\begin{align*} \widehat {G}_1(1)=1, \quad \widehat {G}_1(\!-\!1) =G_1'(\pm 1)=0, \quad \widehat {G}_2(\pm 1) =0, \quad \widehat {G}_2'(1)=0, \quad \widehat {G}_2'(\!-\!1) = 1. \end{align*}

We just provide the details for verifying $\widehat {G}_2'(\!-\!1) = {1}$ as the others are straightforward. For convenience, let us set $c^*\;:\!=\;-{\frac 12}\frac {1}{ \sqrt { \psi ''(\!-\!1)}}$ so that

\begin{align*} {\widehat {G}_2'(\!-\!1) = \frac d {dr} \widehat {G}_2(r)\Big |_{r=-1} }= c^* \Big (\sqrt {2\psi (r)} + (r-1) \frac {\text{d}}{\text{d}r} \sqrt {2\psi (r)}\Big )\Big |_{r=-1} = -2 {c^*}\frac {\text{d}}{\text{d}r}\sqrt {2\psi (r)}\Big |_{r=-1} . \end{align*}

For the argument in the latter expression, using Taylor’s expansion, it holds that

\begin{align*} \psi (r) = \underbrace {\psi (\!-\!1)}_{=0} + \underbrace {\psi '(\!-\!1)(r+1)}_{=0} + \frac 12\psi ''(\!-\!1)(r+1)^2 + o \big ((r+1)^2\big ), \quad r \in {\mathbb{R}}, \end{align*}

whence we infer that

\begin{align*} \frac {\text{d}}{\text{d}r}\sqrt {2\psi (r)}\Big |_{r=-1} = \lim _{r \searrow -1} \frac {\sqrt {\psi ''(\!-\!1)} \,|r+1|}{r+1} = \sqrt { \psi ''(\!-\!1)}. \end{align*}

Thus, using the definition of $\widehat {G}_2$ as stated above, we infer that $\widehat {G}_2'(\!-\!1)=1$ , as claimed. Note that due to the relation $\widehat {G}_3(r) = - \widehat {G}_2(\!-\!r)$ the condition $\widehat {G}_3'(1) = 1$ can be inferred in the same way. In particular, for the quartic potential (2.7), we have $\psi ''(\!-\!1) = 2$ , and so

\begin{equation*} \widehat {G}_2(r) = - \frac {1}{4} (1-r^2)(r-1), \quad \widehat {G}_3(r) = - \widehat {G}_2(\!-\!r) = - \frac {1}{4}(1-r^2)(r+1) \end{equation*}

which fulfil (2.5) and (2.6), respectively. In addition, we set

(2.10) \begin{equation} G_4(r) { = \widehat {G}_4 \Big (\frac {r}{r_c} \Big ), \quad \widehat {G}_4(r) } = 2\psi (r) {,}\quad r\in \mathbb{R}{,} \end{equation}

which clearly satisfies (2.5) and (2.6).

2.3. The sharp interface model

We now present the corresponding sharp interface system related to the Cahn–Hilliard model above. We will later use the method of formally matched asymptotic expansions, see, e.g., [Reference Bänsch, Deckelnick, Garcke and Pozzi4], to derive this free boundary problem. We denote by $\Omega ^\pm =\Omega ^\pm (t)$ the two regions occupied by the droplet and background phases and by $\Sigma =\Sigma (t)$ the evolving interface separating the two phases. In addition, let $m_\pm \;:\!=\;m(\pm 1)$ be the mobilities in the two phases. The free boundary problem corresponding to (2.1) reads as follows:

(2.11a) \begin{align} - m_+ \Delta \mu &= {S_+} - \rho _+ \mu && \! \text{in $\Omega ^+$,} \end{align}

(2.11b) \begin{align} - m_- \Delta \mu & ={S_-} - \rho _- \mu && \! \text{in $\Omega ^-$,} \end{align}

(2.11c) \begin{align} [\mu ]^+_- & = 0 && \! \text{on $\Sigma $,} \end{align}

(2.11d) \begin{align} \mu &= \frac {\gamma \beta \kappa }2 && \! \text{on $\Sigma $,} \end{align}

(2.11e) \begin{align} - 2 \mathcal{V} & = [m\nabla \mu ]^+_- \cdot {{\boldsymbol \nu }} + S_I \qquad && \! \text{on $\Sigma $,} \end{align}

(2.11f) \begin{align} \partial _{{\boldsymbol{n}}} \mu & =0 \qquad\qquad &&\! \text{on $\partial \Omega $,} \end{align}

where $ \rho _\pm = \frac {K_\pm }{\beta \psi ''(\pm 1)}$ , $\gamma$ is a constant depending on the choice of the double-well potential $\psi$ and $S_I$ is the interface reaction constant defined in (4.15), which depends on the choice of the interpolation functions $G_2$ , $G_3$ and $G_4$ . In the case of the quartic potential and $r_c=1$ , we will obtain

\begin{equation*} \rho _\pm = \frac {K_\pm }{2\beta }, \quad \gamma = \frac {2 \sqrt {2}}{3}, \quad S_I = \frac {1}{\sqrt {2}} \Big (K_{+} - K_{-} + \frac {4}{3} L \Big ). \end{equation*}

In the above, we also use $\kappa$ to denote the mean curvature of $\Sigma$ that is given as the sum of the principal curvatures of $\Sigma$ , $\boldsymbol \nu$ is the unit normal to the interface, and $\mathcal{V}$ is the normal velocity of the interface in the direction of the normal $\boldsymbol \nu$ . In addition, for $x\in \Sigma (t)$ and a function $u$ , we define its jump across the interface at $(t,x)$ as

\begin{equation*} [u]^+_-(t,x)\;:\!=\;\lim _{\substack {y\to x\\y\in \Omega ^+(t)}} u(t,y)- \lim _{\substack {y\to x\\y\in \Omega ^-(t)}} u(t,y)\,. \end{equation*}

Further information about the notation used can be found in [Reference Bänsch, Deckelnick, Garcke and Pozzi4].

Notice that the above system is connected to the well-known Mullins–Sekerka free boundary problem [Reference Bänsch, Deckelnick, Garcke and Pozzi4], with the difference that here we have a constant source term $S_I$ on the right-hand side of (2.11e) as well as an affine linear term in the quasi-static diffusion (2.11a) and (2.11b).

2.4. Nondimensionalization for the sharp interface problem

We now perform a nondimensionalization argument in order to identify important dimensionless parameters, but for the subsequent analysis, we mainly work with the original model (2.11). Choosing units $\widetilde x$ , $\widetilde t$ , $\widetilde \mu$ and $\widetilde{\mathcal{V}}$ for length, time, chemical potential and normal velocity we introduce the nondimensional variables

\begin{equation*} \widehat{x} =\frac x{\widetilde{x}},\quad \widehat{t} =\frac t{\widetilde{t}}, \quad \widehat {\mu } =\frac {\mu}{\widetilde{\mu}}, \quad \widehat{\mathcal{V}} = \frac{\mathcal{V}}{\widetilde{\mathcal{V}}}. \end{equation*}

Setting $\widetilde{\mathcal{V}} = \widetilde {x}/\widetilde {t}$ , we now consider the rescaled variant of the system (2.1) on the rescaled domains $\widehat {\Omega }^+ = \widehat {\Omega }^+(t)$ , $\widehat {\Omega }^-= \widehat {\Omega }^-(t)$ , and the corresponding interface $\widehat {\Sigma } = \widehat {\Sigma }(t)$ . Denoting by $\widehat {\nabla }$ , $\widehat {\Delta }$ and $\widehat {\kappa }$ the gradient, Laplacian and mean curvature with respect to $\widehat {x}$ , for the new nondimensional variables, we obtain from (2.11) the following system

\begin{alignat*} {3} - \widehat {\Delta } \widehat {\mu } &= \frac {(\widetilde x)^2}{m_+ \widetilde \mu }{S_+} - \frac {(\widetilde x)^2\rho _+}{m_+} \widehat {\mu } \quad && \text{in $\widehat {\Omega }^+$}, \\ - \widehat {\Delta } \widehat {\mu } &= \frac {(\widetilde x)^2}{m_- \widetilde \mu } {S_-} - \frac {(\widetilde x)^2\rho _-}{m_-} \widehat {\mu } \quad && \text{in $ \widehat {\Omega }^-$,} \\ [\widehat {\mu }]^+_- &= 0 \qquad && \text{on $\widehat {\Sigma }$,} \\ \widehat {\mu } &= \frac {\gamma \beta {\widehat {\kappa }}}{2\widetilde \mu \widetilde x} \qquad && \text{on $\widehat {\Sigma }$,} \\ - 2 { \widehat{\mathcal{V}}} & = \frac {\widetilde t\widetilde \mu } {(\widetilde x)^2}[m\widehat {\nabla } \widehat {\mu }]^+_- \cdot {\widehat {\boldsymbol \nu }} +\frac {\widetilde t}{\widetilde x} S_I \qquad && \text{on ${\widehat {\Sigma }}$,} \\ \partial _{{\boldsymbol{n}}} \widehat {\mu } &= 0 \qquad && \text{on ${\partial \widehat {\Omega }}$.} \end{alignat*}

To obtain simple nondimensional equations, and on assuming that $\rho _-\gt 0$ and $S_- \gt 0$ , we set

\begin{equation*} \widetilde x =\sqrt {\frac {{m_-}}{\rho _-}}, \quad \widetilde \mu = \frac {(\widetilde x)^2 S_- }{m_-}= \frac {S_-}{\rho _-}, \quad \widetilde t = \frac {(\widetilde x)^2}{\widetilde \mu m_-} =\frac 1{S_-} . \end{equation*}

We now define the nondimensional parameter

\begin{equation*} \beta ^* =\frac {\beta }{\widetilde x} \frac 1{\widetilde \mu }= \frac {\beta \rho _-^{\frac 32}}{m_-^{\frac 12} S_-}= \frac {c_{l}}{\widetilde x}, \end{equation*}

where we call

(2.13) \begin{equation} c_{l}= \beta \frac {\rho _-}{S_-}{,} \end{equation}

in analogy to solidification problems, the modified capillary length. In addition, we introduce the relative mobility $m^*$ , the relative reaction coefficients $S^*$ and $\rho ^*$ , and the nondimensional interface reaction term $S_I^*$ as follows

\begin{equation*} m^*=\frac {m_+} {m_-},\quad S^*=\frac {S_+} {S_-},\quad \rho ^*=\frac {\rho _+} {\rho _-}, \quad S_I^* =\frac {\widetilde t}{\widetilde x}S_I= \frac {1}{S_-} \sqrt {\frac {\rho _-}{m_-}} S_I. \end{equation*}

We then obtain, dropping the hat notation for convenience,

\begin{alignat*} {3} - m^*\Delta \mu &= S^* - \rho ^* \mu \quad && \text{in $ \Omega ^+$,} \\ - \Delta \mu &= 1 - \mu \quad && \text{in $ \Omega ^-$}, \\ [ \mu ]^+_- &= 0 \qquad && \text{on $ \Sigma $,} \\ \mu & = \frac {\gamma \beta ^*{\kappa }}2 \qquad & &\text{on $ \Sigma $,} \\ - 2 { \mathcal{V}} &= m^*\nabla \mu _+ \cdot { {\boldsymbol \nu }} - \nabla \mu _- \cdot { {\boldsymbol \nu }} + S_I^* \qquad && \text{on $\Sigma $,} \\ \partial _{{\boldsymbol{n}}} \mu &=0 \qquad && \text{on $\partial \Omega $,} \end{alignat*}

and observe that the evolution critically depends on the nondimensional number $\beta ^*$ , which relates the influence of surface tension $\beta$ to a generalised supersaturation $\frac {S_-}{\rho _-} \sqrt {\frac {m_-}{\rho _-}}$ stemming from chemical reactions.

3.1. Notation

Let $\Omega$ be a bounded domain in $\mathbb{R}^d$ , $d\in \{1,2,3\}$ . The Lebesgue measure of $\Omega$ is denoted by $|\Omega |$ , while the Hausdorff measure of the boundary $\partial \Omega$ is denoted by $|\partial \Omega |$ .

For any Banach space $X$ , its norm is represented as $\| \cdot \|_X$ , its dual space as $X^*$ , and the duality pairing between $X^*$ and $X$ is denoted by $\langle \cdot , \cdot \rangle _X$ . In the case where $X$ is a Hilbert space, the inner product is denoted by $(\cdot ,\cdot )_X$ .

For each $1 \leq p \leq \infty$ and $k \geq 0$ the standard Lebesgue and Sobolev spaces defined on $\Omega$ are denoted as $L^p(\Omega )$ and $W^{k,p}(\Omega )$ , with their respective norms $\| \cdot \|_{L^p(\Omega )}$ and $\| \cdot \|_{W^{k,p}(\Omega )}$ . For simplicity we may often use $\| \cdot \|_{L^p}$ instead of $\| \cdot \|_{L^p(\Omega )}$ , and employ similar shorthand notation for other norms. We adopt the convention $H^k(\Omega ) \;:\!=\; W^{k,2}(\Omega )$ for all $k \in \mathbb{N}$ , and denote the mean value of a functional $h \in {(H^{1}(\Omega ))^*}$ as

\begin{equation*} h_\Omega \;:\!=\; \frac {1}{|\Omega |} \langle h,1 \rangle _{H^1}. \end{equation*}

We now introduce a tool commonly employed in the investigation of problems associated with equations of Cahn–Hilliard type. Given $\phi \in {(H^{1}(\Omega ))^*}$ , we seek $u \in H^{1}(\Omega )$ such that

(3.1) \begin{align} \int _{\Omega } \nabla u \cdot \nabla v = \langle \phi , v \rangle _{H^1}, \quad v \in H^{1}(\Omega ). \end{align}

This corresponds to the standard weak formulation of the homogeneous Neumann problem for the Poisson equation $-\Delta u = \phi$ for $\phi \in L^{2}(\Omega )$ . The solvability of (3.1) for $\phi \in {(H^{1}(\Omega ))^*}$ relies on the condition that $\phi$ possesses a zero mean value, that is, $\phi _\Omega =0$ . If this condition is satisfied, a unique solution with a zero mean value exists, and the operator $\mathcal{N}\;:\; \text{dom}(\mathcal{N})= \{\phi \in {(H^{1}(\Omega ))^*} \;:\; \phi _\Omega = 0\} \to \{u \in H^{1}(\Omega ) \;:\; u_\Omega = 0\}$ defined by mapping $\phi$ to the unique solution $u$ to (3.1) with $u_\Omega = 0$ is well-defined. This operator yields an isomorphism between the mentioned spaces. Additionally, the norm

(3.2) \begin{align} \phi \mapsto \|\phi \|_{*}^2 \;:\!=\; \|\nabla \mathcal{N} (\phi - {\phi _\Omega })\|^2_{L^2} + |\phi _\Omega |^2, \quad \phi \in {(H^{1}(\Omega ))^*}, \end{align}

is proven to define a Hilbert norm in $(H^{1}(\Omega ))^*$ that is equivalent to the standard dual norm. From these definitions, it directly follows that

\begin{equation*} \begin{alignedat}{3} &\int _{\Omega } \nabla \mathcal{N}\phi \cdot \nabla v = \langle \phi , v \rangle _{H^1} \quad &&\text{for every $\phi \in \text{dom}(\mathcal{N})$ and $v \in H^{1}(\Omega )$},\\ &\langle \phi , \mathcal{N}\zeta \rangle _{H^1} = \langle \zeta , \mathcal{N}\phi \rangle _{H^1}\quad && \text{for every $\phi , \zeta \in \text{dom}(\mathcal{N})$}, \\ &\langle \phi , \mathcal{N}\phi \rangle _{H^1}= \mathopen \| \nabla \mathcal{N}\phi \mathclose \|_{L^2}^2 = \|\phi \|_{*}^2 \quad && \text{for every $\phi \in \text{dom}(\mathcal{N})$}. \end{alignedat} \end{equation*}

Moreover, it is established that, see e.g. [Reference Garcke and Lam22],

\begin{align*} \int _{0}^{t} \langle \partial _t v(s) , \mathcal{N} v(s) \rangle _{H^1} \,\text{d}s = \int _{0}^{t} \langle v(s) , \mathcal{N}(\partial _t v(s)) \rangle _{H^1} \,\text{d}s = \frac {1}{2} \|v(t)\|_{*}^2 - \frac {1}{2} \|v(0)\|_{*}^2 \end{align*}

for every $t \in [0, T]$ and $v \in H^{1}(0,T;{{(H^{1}(\Omega )^*))}}$ such that $v_\Omega (t) = 0$ for every $t \in [0,T]$ .

3.2. Assumptions

For the well-posedness, we require the following assumptions:

1. (A1) The symbols $K_\pm$ and $S_\pm$ denote real-valued constants, whereas $\beta$ and $\varepsilon$ denote positive constants.

2. (A2) The potential $\psi \;:\; {\mathbb{R}} \to [0,\infty )$ is twice differentiable and can be decomposed as $\psi =\psi _1 + \psi _2$ , with $\psi _1$ convex and $\psi _2$ a quadratic perturbation. Namely, we require that there exist a positive constant $C_1$ such that it holds

   \begin{align*} |\psi _2(r)| \leq C_1 (|r|^2 +1){,} \quad {r \in {\mathbb{R}},} \end{align*}
   and in addition, we require
   \begin{align*} \forall \eta \gt 0 \, \exists \,C_\eta \;:\; {\forall r \in {\mathbb{R}}} \quad |\psi '(r)| & \leq \eta \psi (r) + C_\eta . \end{align*}

3. (A3) We require the source $S_\varepsilon$ to be Lipschitz continuous. Consequently, there exists a positive constant $C_S$ such that

   \begin{align*} |S_\varepsilon (r)| \leq C_S (|r|+1), \quad r \in \mathbb{R}. \end{align*}

4. (A4) The mobility function $m\;:\;{\mathbb{R}} \to {\mathbb{R}}$ is continuous and there exist positive constants $m_*$ and $M^*$ such that

   \begin{align*} 0 \lt m_* \leq m(r) \leq M^*, \quad r \in {\mathbb{R}}. \end{align*}

Let us notice that assumption (3.3) restricts the class of admissible double-well potentials, but still includes the quartic potential in (2.7). For this latter, referring to (A2), we employ the splitting

\begin{align*} \psi _1(r) = \frac 14 r^4, \quad \psi _2(r) = \frac 14 (1- 2 r^2), \quad r \in {\mathbb{R}}. \end{align*}

Here is our main result.

Theorem 3.1. Suppose that (A1)–(A4) are fulfilled. Then, for every given ${{\varphi }}_0\in {H^{1}(\Omega )}$ there exists a weak solution $({{\varphi }},\mu )$ to (2.1) such that

\begin{align*} {{\varphi }} & \in H^{1}(0,T;\; {(H^{1}(\Omega )^*))} \cap L^{\infty }(0,T;\; {H^{1}(\Omega )} ) \cap L^{2}(0,T;\; {H^{2}(\Omega )}), \\ \mu & \in L^{2}(0,T;\; {H^{1}(\Omega )}), \end{align*}

and satisfying

\begin{align*} & \langle \partial _t {{\varphi }}, v\rangle _{{H^1}} + {\int _\Omega } m({{\varphi }}) \nabla \mu \cdot \nabla v = {\int _\Omega } S_\varepsilon ({{\varphi }}) v, \\ & {\int _\Omega } \mu v = \beta \varepsilon {\int _\Omega } \nabla {{\varphi }} \cdot \nabla v + \frac \beta \varepsilon {\int _\Omega } \psi '({{\varphi }})v, \end{align*}

for every $v \in H^{1}(\Omega )$ and almost every $t \in (0,T)$ , along with attainment of the initial condition $\varphi (0) = \varphi _0$ holding for almost every $x \in \Omega$ .

Moreover, let $\{({{\varphi }}_i,\mu _i)\}_i$ , $i=1,2$ , denote two arbitrary solutions to (2.1) associated with initial data ${{\varphi }}_{0,i}\in H^{1}(\Omega )$ , $i=1,2$ and to a constant mobility $m$ . In addition, let us assume that there exists an exponent $p \in [1,7)$ such that $\psi _1'$ satisfies the following pointwise growth condition

(3.3) \begin{align} |\psi _1'(r)-\psi _1'(s)| \leq C (1+ |r|^p+|s|^p) |r-s|, \quad r,s \in {\mathbb{R}}{.} \end{align}

Then, it holds that

\begin{align*} & \mathopen \| ({{\varphi }}_1-{{\varphi }}_2)-({{\varphi }}_1-{{\varphi }}_2)_\Omega \mathclose \|_{L^{\infty }(0,T;\; {(H^{1}(\Omega ))^*}) \cap L^{2}(0,T;\; {H^{1}(\Omega )}) } + \mathopen \| ({{\varphi }}_1)_\Omega -({{\varphi }}_2)_\Omega \mathclose \|_{L^\infty (0,T)} \\ & \quad \leq {C{^*}} \big ( {\mathopen \| ({{\varphi }}_{0,1} -{{\varphi }}_{0,2}) - (({{\varphi }}_{0,1})_\Omega - ({{\varphi }}_{0,2})_\Omega )\mathclose \|_*} +|({{\varphi }}_{0,1})_\Omega -({{\varphi }}_{0,2})_\Omega | \big ), \end{align*}

for a positive constant ${C}{^*}$ just depending on $\Omega$ , $T$ and the nonlinearity $\psi$ . Consequently, under these conditions, the weak solution to (2.1) is unique.

Proof of Theorem 3.1. We begin with the existence part of the theorem. In the subsequent discussion, we adopt a formal approach, leveraging standard procedures to derive estimates for the solution. While these computations are formal, they suggest that the same estimates can be rigorously applied to the $k$ -dimensional system obtained via a Faedo–Galerkin scheme, constructed using the first $k$ eigenfunctions of the Laplace operator with homogeneous Neumann boundary conditions. These bounds can then be used to pass to the limit as $k\to \infty$ , thereby constructing a solution to the problem that satisfies (2.1). A rigorous proof can be easily adapted within this approximation framework, see for instance [Reference Garcke and Lam21] for details of a Faedo–Galerkin approximation applied to a similar Cahn–Hilliard model with source terms. Furthermore, in what follows, since we will need to make numerous estimates, we will use the symbol $C$ to denote any nonnegative constant depending on the system’s data, which may change its value from line to line, or even within the same line.

First estimate: Before starting with the proper estimates, let us introduce a new function related to the phase-dependent mobility function $m$ . We set

\begin{align*} M'(z) = \int _0^z \frac 1{m(s)} ds, \quad z \in {\mathbb{R}}, \end{align*}

and set $M(r)$ to be the antiderivative of the above function with $M(0) = 0$ . We note that $M \in C^2({\mathbb{R}})$ and $M''(r)=\frac 1 {m(r)}$ , $r \in {\mathbb{R}}$ . Due to the bounds on $m$ in (A4), we also have that there exist two positive constants $c_1$ and $c_2$ such that

\begin{align*} c_1(1+|r|) \leq |M'(r)| \leq c_2(1+|r|), \quad r \in {\mathbb{R}}. \end{align*}

We then test (2.1a) with $M'({{\varphi }})$ , (2.1b) with $-\Delta {{\varphi }}$ and add the resulting identities to obtain

\begin{align*} \frac 12 \frac d {dt} \mathopen \| M({{\varphi }})\mathclose \|_{L^1} + \beta \varepsilon \mathopen \| \Delta {{\varphi }}\mathclose \|^2_{L^2} + \frac \beta \varepsilon {\int _\Omega } \psi _1''({{\varphi }})|\nabla {{\varphi }}|^2 ={\int _\Omega } S_\varepsilon ({{\varphi }}) M'({{\varphi }}) - \frac \beta \varepsilon {\int _\Omega } \psi _2''({{\varphi }})|\nabla {{\varphi }}|^2{,} \end{align*}

where we notice that the third term on the left-hand side is nonnegative due to (A2). For the second term on the right-hand side, we use that $\psi _2''$ is bounded and the Neumann boundary condition to see

\begin{equation*} - \frac {\beta }{\varepsilon } {\int _\Omega } \psi _2''({{\varphi }}) |\nabla {{\varphi }}|^2 \leq C \|\nabla {{\varphi }} \|_{L^2}^2 \leq C \| \Delta {{\varphi }} \|_{L^2} \| {{\varphi }} \|_{L^2} \leq \frac {\beta \varepsilon }{2} \| \Delta {{\varphi }} \|_{L^2}^2 + C \| {{\varphi }} \|_{L^2}^2. \end{equation*}

Now, since $M'$ is growing linearly, using (A3), we readily infer that

\begin{equation*} {\int _\Omega } S_\varepsilon ({{\varphi }}) M'({{\varphi }}) \leq C \Big ( 1 + \| {{\varphi }} \|_{L^2}^2 \Big ). \end{equation*}

Hence, we obtain

\begin{equation*} \frac {d}{dt} \| M({{\varphi }}) \|_{L^1} + \beta \varepsilon \| \Delta {{\varphi }} \|_{L^2}^2 \leq C \Big ( 1 + \| {{\varphi }} \|_{L^2}^2 \Big ) \leq C \Big ( 1 + \| M ({{\varphi }}) \|_{L^1} \Big ). \end{equation*}

Using Grönwall’s inequality then produces

\begin{align*} \mathopen \| M({{\varphi }})\mathclose \|_{L^{\infty }(0,T;\; {L^{1}(\Omega )}) } + \mathopen \| \Delta {{\varphi }}\mathclose \|_{L^{2}(0,T;\; {L^{2}(\Omega )}) } \leq C. \end{align*}

Besides, since $M'$ is growing linearly, $M$ grows quadratically, so that from the above inequality we obtain, using also elliptic regularity theory, that

\begin{align*} \mathopen \| {{\varphi }}\mathclose \|_{L^{\infty }(0,T;\; {L^{2}(\Omega ))} \cap L^{2}(0,T;\; {H^{2}(\Omega )})} \leq C. \end{align*}

Second estimate: Next, recalling (A2), upon testing (2.1b) with $\frac 1{|\Omega |}$ yields

\begin{align*} \mathopen \| \mu _\Omega \mathclose \|_{L^2(0,T)} \leq C, \end{align*}

where we also used the Sobolev embedding $H^2(\Omega ) {\hookrightarrow } L^{\infty }(\Omega )$ .

Third estimate: We now perform the usual energy estimate in the context of the Cahn–Hilliard equation by testing (2.1a) with $\mu$ , (2.1b) with $-\partial _t {{\varphi }}$ and adding the resulting identities to obtain

\begin{align*} \frac {\text{d}}{\text{d}t} \mathcal{E}({{\varphi }}) + {\int _\Omega } m({{\varphi }}) |\nabla \mu |^2 = {\int _\Omega } S_\varepsilon ({{\varphi }}) \mu = {\int _\Omega } S_\varepsilon ({{\varphi }}) (\mu -\mu _\Omega ) + {\int _\Omega } S_\varepsilon ({{\varphi }})\mu _\Omega = I_1 + I_2. \end{align*}

Now, employing the Poincaré–Wirtinger inequality, the Young inequality and (A3), we infer that

\begin{align*} |I_1| &\leq C ( \mathopen \| {{\varphi }}\mathclose \|_{L^2}+1)\mathopen \| \nabla \mu \mathclose \|_{L^2} \leq \frac 12 \mathopen \| \nabla \mu \mathclose \|^2_{L^2} + C {(}\mathopen \| {{\varphi }}\mathclose \|^2_{L^2}{+1)}, \\ |I_2| &\leq C ( \mathopen \| {{\varphi }}\mathclose \|_{L^2}+1)|\mu _\Omega | \leq C {(}\mathopen \| {{\varphi }}\mathclose \|^2_{L^2}{+1)} + |\mu _\Omega |^2. \end{align*}

Thus, using the previous estimates readily entails that

\begin{align*} \mathopen \| {{\varphi }}\mathclose \|_{L^{\infty }(0,T;\; {H^{1}(\Omega )})} + \mathopen \| \mu \mathclose \|_{L^{2}(0,T;\; {H^{1}(\Omega )})} \leq C. \end{align*}

Fourth estimate: Finally, from the boundedness of the mobility (A4) and the linear growth (A3) of $S_\varepsilon$ , it is standard to infer from the previous estimates and the weak formulation (2.1a) that

\begin{equation*} \| \partial _t {{\varphi }} \|_{{(H^{1}(\Omega ))^*}} \leq C \Big (\| \nabla \mu \|_{L^2} + \| {{\varphi }} \|_{L^2} + 1 \Big ), \end{equation*}

so that

\begin{align*} \mathopen \| \partial _t {{\varphi }}\mathclose \|_{L^{2}(0,T;\; {(H^{1}(\Omega ))^*})} \leq C. \end{align*}

This concludes the existence part of the proof.

Moving to the uniqueness part, recalling that the mobility is assumed to be constant, and with loss of generality, we set $m \equiv 1$ . Now we consider two solutions $\{({{\varphi }}_i,\mu _i)\}_i$ , $i=1,2$ , associated with initial data ${{\varphi }}_{0,i}$ , $i=1,2$ . Then, we introduce the notation

\begin{align*} {{\varphi }}= {{\varphi }}_1-{{\varphi }}_2, \quad \mu = \mu _1- \mu _2, \quad {{\varphi }}_0={{\varphi }}_{0,1}-{{\varphi }}_{0,2}, \end{align*}

and consider the system (2.1) written for the differences. Considering the difference of (2.1a), recalling the Lipschitz continuity of $S_\varepsilon$ , and testing it with ${{\varphi }}_\Omega {=({{\varphi }}_{1})_\Omega -({{\varphi }}_{2})_\Omega }$ produces

(3.4) \begin{align} \frac {\text{d}}{\text{d}t}\bigg ( \frac 12 |{{\varphi }}_\Omega |^2 \bigg ) \leq C |{{\varphi }}_\Omega |^2 + C \mathopen \| {{\varphi }}\mathclose \|^2_{L^2}. \end{align}

Then, we consider the difference of (2.1a) minus its mean value and test it with $\mathcal{N}({{\varphi }}-{{\varphi }}_\Omega )$ , the difference of (2.1b) and test it with $-({{\varphi }}-{{\varphi }}_\Omega )$ and, upon adding the resulting equalities, we obtain that

(3.5) \begin{equation} \begin{aligned} & \frac 12 \frac {\text{d}}{\text{d}t}\mathopen \| {{\varphi }}-{{\varphi }}_\Omega \mathclose \|^2_* + {\int _\Omega } {(\mu -\mu _\Omega )({{\varphi }}-{{\varphi }}_\Omega )} \\ & \quad = {\int _\Omega } \big ( S_\varepsilon ({{\varphi }}_1) - S_\varepsilon ({{\varphi }}_2) - (S_\varepsilon ({{\varphi }}_1) - S_\varepsilon ({{\varphi }}_2) )_\Omega \big ) \mathcal{N}({{\varphi }}-{{\varphi }}_\Omega ) \\ & \quad \leq C \mathopen \| {{\varphi }}\mathclose \|^2_{L^2} + C \mathopen \| {{\varphi }}-{{\varphi }}_\Omega \mathclose \|^2_* . \end{aligned} \end{equation}

Besides, it holds that

\begin{align*} {\int _\Omega } {(\mu -\mu _\Omega )({{\varphi }}-{{\varphi }}_\Omega )} = \beta \varepsilon \mathopen \| \nabla {{\varphi }}\mathclose \|^2_{L^2} + \frac \beta \varepsilon {\int _\Omega } (\psi '({{\varphi }}_1)-\psi '({{\varphi }}_2))({{\varphi }}-{{\varphi }}_\Omega ). \end{align*}

For the last term, using (A2), we have

\begin{align*} & \frac \beta \varepsilon {\int _\Omega } (\psi '({{\varphi }}_1)-\psi '({{\varphi }}_2))({{\varphi }}-{{\varphi }}_\Omega ) \\ & \quad = \underbrace {\frac \beta \varepsilon {\int _\Omega } (\psi _1'({{\varphi }}_1)-\psi _1'({{\varphi }}_2)){{\varphi }}}_{\geq 0} + \underbrace {\frac \beta \varepsilon {\int _\Omega } (\psi '_2({{\varphi }}_1)-\psi '_2({{\varphi }}_2)){{\varphi }}}_{\geq -C \mathopen \| {{\varphi }}\mathclose \|^2_{L^2}} \\ & \qquad -\frac \beta \varepsilon {\int _\Omega } (\psi _1'({{\varphi }}_1)-\psi _1'({{\varphi }}_2)){{\varphi }}_\Omega -\underbrace {\frac \beta \varepsilon {\int _\Omega } (\psi _2'({{\varphi }}_1)-\psi _2'({{\varphi }}_2)){{\varphi }}_\Omega }_{\geq - C (\mathopen \| {{\varphi }}\mathclose \|^2_{L^2} + |{{\varphi }}_\Omega |^2)}. \end{align*}

We then move the third term on the right-hand side of this latter identity to the right-hand side of the estimate (3.5) and continue with the estimation. Recalling the growth condition in (3.3) and the continuous embedding $H^{1}(\Omega ) {\hookrightarrow } L^{6}(\Omega )$ holding in three dimensions, we apply Hölder’s inequality to bound

\begin{align*} & \Big |\frac \beta \varepsilon {\int _\Omega } (\psi _1'({{\varphi }}_1)-\psi _1'({{\varphi }}_2)){{\varphi }}_\Omega \Big | \leq C \left (1+ \mathopen \| |{{\varphi }}_1|^p\mathclose \|_{L^{\frac 65}} + \mathopen \| |{{\varphi }}_2|^p\mathclose \|_{L^{\frac 65}} \right ) \mathopen \| {{\varphi }}\mathclose \|_{L^6} |{{\varphi }}_\Omega | \\ & \quad \leq \frac {\beta \varepsilon }{4} \mathopen \| {{\varphi }}\mathclose \|^2_{H^1} + C \Big (1+ \mathopen \| {{\varphi }}_1\mathclose \|_{L^{\frac {6p}5}}^{2p} + \mathopen \| {{\varphi }}_2\mathclose \|_{L^{\frac {6p}5}}^{2p} \Big ) |{{\varphi }}_\Omega |^2. \end{align*}

Given that

\begin{align*} {{\varphi }}_i \in L^{\infty }(0,T;\; {H^{1}(\Omega )}) \cap L^{2}(0,T;\; {H^{2}(\Omega )}) {\hookrightarrow } L^{\frac {4q}{q-6}}(0,T;\; {L^{q}(\Omega )}) \end{align*}

for $i=1,2$ and $q \in [6,\infty )$ in three spatial dimensions, with the convention that $\frac {4q}{q-6}\;:\!=\; +\infty$ when $q=6$ , we deduce that any exponent up to $p = 5$ is admissible as we notice that

\begin{equation*} t \mapsto \Lambda (t)\;:\!=\;(1+ \mathopen \| {{\varphi }}_1(t)\mathclose \|_{L^6}^{10} + \mathopen \| {{\varphi }}_2(t)\mathclose \|_{L^6}^{10} ) \in L^{\infty }(0,T). \end{equation*}

Besides, selecting $q = \frac {6p}5$ in the above interpolation embedding, we infer that the resulting time exponent $\frac {4p}{p/5}$ is strictly bigger than $2p$ for any $p \in (5,7)$ , entailing that $t \mapsto \Lambda (t) \in L^1(0,T)$ for any $p \in (5,7)$ . Regarding the term $C \mathopen \| {{\varphi }}\mathclose \|^2$ on the right-hand side of (3.5), we recall the norm $\| \cdot \|_*$ defined in (3.2) that is equivalent to the norm in $(H^1(\Omega ))^*$ , and apply the Poincaré inequality to deduce that

\begin{align*} C \mathopen \| {{\varphi }}\mathclose \|^2_{L^2} & \leq C \| {{\varphi }} - {{\varphi }}_\Omega \|_{L^2}^2 + C |{{\varphi }}_\Omega \|_{L^2}^2 \\ & \leq C \| {{\varphi }} - {{\varphi }}_{\Omega } \|_{H^1} \| {{\varphi }} - {{\varphi }}_\Omega \|_* + C |{{\varphi }}_\Omega |^2 \\ & \leq C \| \nabla {{\varphi }} \|_{L^2} \| {{\varphi }} - {{\varphi }}_\Omega \|_* + C |{{\varphi }}_\Omega |^2 \\ & \leq \frac {\beta \varepsilon }{4} \mathopen \| \nabla {{\varphi }}\mathclose \|^2_{L^2} + C \mathopen \| {{\varphi }}-{{\varphi }}_\Omega \mathclose \|^2_* + C |{{\varphi }}_\Omega |^2. \end{align*}

Thus, rearranging the terms and adding the above estimates, we end up with

\begin{align*} & \frac 12 \frac {\text{d}}{\text{d}t} \Big (\mathopen \| {{\varphi }}-{{\varphi }}_\Omega \mathclose \|^2_* + |{{\varphi }}_\Omega |^2 \Big ) {+\frac {\beta \varepsilon }{2}} \mathopen \| \nabla {{\varphi }}\mathclose \|^2_{L^2} \leq C \mathopen \| {{\varphi }}-{{\varphi }}_\Omega \mathclose \|^2_* {+ C (\Lambda +1 )|{{\varphi }}_\Omega |^2.} \end{align*}

Finally, we integrate over time and employ Grönwall’s inequality to conclude the proof.

4.1. Outer expansion

In what follows, we suppose that the solution variables ${{\varphi }}_\varepsilon$ and $\mu _\varepsilon$ , far away from the interface, can be expressed as

\begin{align*} {{\varphi }}_\varepsilon = \sum _{i=0}^{\infty } \varepsilon ^i {{\varphi }}_i, \quad \mu _\varepsilon = \sum _{i=0}^{\infty } \varepsilon ^i \mu _i. \end{align*}

Then, we consider (2.1b) to leading order $\varepsilon ^{-1}$ to infer that

\begin{align*} - \beta \psi '({{\varphi }}_0)=0. \end{align*}

Here, we accounted for the conditions $\psi (\pm 1)=\psi '(\pm 1)=0$ in (4.1). Since ${{\varphi }}_0=0$ is an unstable solution, this leads to ${{\varphi }}_0=\pm 1$ . We then set

\begin{align*} \Omega ^\pm \;:\!=\;\{x \in \Omega \;:\; {{\varphi }}_0(x) = \pm 1\}. \end{align*}

To the next order $\varepsilon ^0$ , (2.1b) yields

\begin{align*} \mu _0 = \beta \psi ''({{\varphi }}_0){{\varphi }}_1. \end{align*}

Besides, to the same order, since $ {{\varphi }}_0$ equals $\pm 1$ in $\Omega ^\pm$ , we infer from (2.1a) that we have

\begin{align*} - \textrm {div}(m({{\varphi }}_0)\nabla \mu _0) = S_\pm - K_\pm {{\varphi }}_1. \end{align*}

Combining these two equations and recalling $\rho _{\pm } = \frac {K_{\pm }}{\beta \psi ''(\pm 1)}$ yields

\begin{align*} - m_+ \Delta \mu _{0} = S_{+} - \rho _+ \mu _0 & \quad \text{ in } \Omega ^+, \\ - m_- \Delta \mu _{0} = S_{-} - \rho _- \mu _0 & \quad \text{ in } \Omega ^-{,} \end{align*}

where we recall $m_\pm = m(\pm 1)$ .

4.2. Inner expansion and matching conditions

To explore the behaviour of $\Sigma _\varepsilon =\{{{\varphi }}_\varepsilon =0\}$ as $\varepsilon \to 0$ , we introduce a new coordinate system. Let $d$ denote the signed distance function to limiting hypersurface $\Sigma _0$ , and define $z= \frac d\varepsilon$ the rescaled distance variable. Additionally, we select $d$ in such a way that $d({x})\gt 0$ in $\Omega ^+$ and $d({x})\lt 0$ in $\Omega ^-$ . Thus, it follows that $\nabla d = {\boldsymbol \nu }$ is the unit normal of $\Sigma _0$ and points from $\Omega ^-$ towards $\Omega ^+$ . Next, let us consider a parametrisation of $\Sigma _0$ by arc-length, denoted as $g(t, s)$ . Then, within a tubular neighbourhood of $\Sigma _0$ , for a sufficiently smooth function $f(t,x)$ , we obtain the reparametrisation rule

\begin{align*} f(t,x) = f(t, g(t,s)+\varepsilon z {\boldsymbol \nu } (g(t,s))) \;=\!:\; F(t,s,z). \end{align*}

Under the smoothness assumption for $\Sigma _0$ , we can express $s$ and $z$ as functions of $(t,x)$ , so that the following identities relating partial derivatives of $f$ and $F$ can be derived, see e.g. [Reference Eck, Garcke and Knabner17, Section 7.9, p. 468] or [Reference Garcke and Stinner26, Appendix B] for a derivation:

(4.2) \begin{equation} \begin{aligned} \partial _t f & = - \frac 1 \varepsilon \mathcal{V} \partial _z F + \text{ h.o.t.}, \\ \nabla _x f & = \frac 1\varepsilon \partial _z F {\boldsymbol \nu } + \nabla _{\Sigma _0} F + \text{ h.o.t.}, \quad \Delta _x f = \frac 1{\varepsilon ^2} \partial _{zz} F - \frac 1\varepsilon \kappa \partial _z F + \text{ h.o.t.}, \end{aligned} \end{equation}

where h.o.t. denotes higher-order terms in $\varepsilon$ . Here, $\nabla _{\Sigma _0}$ stands for the surface gradient on $\Sigma _0$ , $\kappa =-\textrm {div}_{\Sigma _0} {\boldsymbol \nu }$ for the corresponding mean curvature. The solution in the inner region is assumed to possess the following expansion

\begin{align*} \Phi _\varepsilon = \sum _{i=0}^{\infty } \varepsilon ^i \Phi _i, \quad \Lambda _\varepsilon = \sum _{i=0}^{\infty } \varepsilon ^i \Lambda _i. \end{align*}

The postulated convergence of the level sets $\Sigma _\varepsilon = \{{{\varphi }}_\varepsilon =0\}$ to the hypersurface $\Sigma _0$ translates to the condition

(4.3) \begin{align} \Phi _0(t,s,z=0)=0. \end{align}

Furthermore, we assume

\begin{align*} \lim _{z \to + \infty } \Phi _\varepsilon (t,s,z) =1, \quad \lim _{z \to -\infty } \Phi _\varepsilon (t,s,z) =-1. \end{align*}

With reference to [Reference Eck, Garcke and Knabner17, Reference Garcke, Lam, Sitka and Styles23, Reference Garcke and Stinner26], we employ the following matching conditions

(4.4) \begin{align} & \lim _{z \to \pm \infty } \Phi _0 (t,s,z) = {{\varphi }}_0^{\pm }(t,{x}), \quad && \lim _{z \to \pm \infty } \Lambda _0 (t,s,z) = \mu _0^{\pm }(t,{x}), \end{align}

(4.5) \begin{align} & \lim _{z \to \pm \infty } \partial _z \Phi _0 (t,s,z) = 0, \quad && \lim _{z \to \pm \infty } \partial _z \Lambda _0 (t,s,z) = 0, \end{align}

(4.6) \begin{align} & \lim _{z \to \pm \infty } \partial _z \Phi _1 (t,s,z) = \nabla {{\varphi }}_0^{\pm }(t,{x})\cdot {\boldsymbol \nu }, \quad && \lim _{z \to \pm \infty } \partial _z \Lambda _1 (t,s,z) = \nabla \mu _0^{\pm }(t,{x})\cdot {\boldsymbol \nu }, \end{align}

where ${{\varphi }}_0^\pm (t,{x})\;:\!=\; \lim _{\delta \to 0} {{\varphi }}_0 (t,{x}\pm \delta {\boldsymbol \nu })$ , ${x} \in \Sigma _0$ , and similarly for $\mu _0$ . This allows us to introduce the corresponding jump across $\Sigma _0$ by using the notation

\begin{align*} [{{\varphi }}_0]^+_-=[{{\varphi }}_0(t,{x}))]^+_-\;:\!=\;{{\varphi }}_0^+ (t,{x})-{{\varphi }}_0^- (t,{x}), \end{align*}

and similarly for $\mu _0$ .

4.3. Leading order expansions in the interfacial region

From (2.1b), to leading order $\varepsilon ^{-1}$ , we find

\begin{align*} \partial _{zz} \Phi _0 - \psi '(\Phi _0)=0. \end{align*}

Using (4.3), we observe that this entails $\Phi _0$ is just a function of $z$ , resulting in the ODE relation

(4.7) \begin{align} \partial _{zz} \Phi _0(z)-\psi '(\Phi _0(z))=0, \quad \Phi _0(0)=0, \quad \Phi _0(\pm \infty ) = \pm 1. \end{align}

Upon testing (4.7) with $\partial _{z}\Phi _0$ , we obtain the so-called equipartition of energy:

(4.8) \begin{align} \frac 12 |\partial _{z}\Phi _0(z)|^2 = \psi (\Phi _0(z)), \quad z \in {\mathbb{R}}. \end{align}

Hence, we find the identity

(4.9) \begin{align} \int _{-\infty }^{+ \infty } |\partial _{z}\Phi _0(z)|^2\,\text{d}z = \int _{-\infty }^{+ \infty } 2 \psi (\Phi _0(z))\,\text{d}z = \int _{-1}^1 \sqrt {2 \psi (s)} \,\text{d}s \;=\!:\; \gamma . \end{align}

From $\psi '(\!-\!z)=-\psi '(z)$ , we see that $-\Phi _0(z)=\Phi _0(\!-\!z)$ . Let us point out that, in the special case of $\psi$ being the quartic potential (2.7), i.e., $\psi (r) = \frac 14 (1-r^2)^2,$ we obtain

\begin{equation*} \gamma = {\frac 1{\sqrt {2}}} \int _{-1}^1 (1-s^2 ) \,\text{d}s = \frac {2\sqrt {2}}3 . \end{equation*}

Next, considering (2.1a) to order $\varepsilon ^{-2}$ produces

\begin{align*} \partial _{z}(m(\Phi _0)\partial _z \Lambda _0 ) = 0. \end{align*}

Upon integration and using the matching condition (4.5) to $\Lambda _0$ , we infer that

(4.10) \begin{align} m(\Phi _0(z))\partial _{z} \Lambda _0(t,s,z) =0, \quad {z \in {\mathbb{R}}}. \end{align}

Since $m(\Phi _0)\gt 0$ , (4.10) implies $\partial _{z} \Lambda _0(t,s,z) = 0$ . We integrate once more and use the matching condition (4.4) to obtain the jump condition

\begin{equation*} [\mu _0]^+_-=0. \end{equation*}

4.4. Higher-order expansions in the interfacial region

Moving to $\varepsilon ^0$ order in (2.1b), we derive that

(4.11) \begin{align} \Lambda _0 = \beta \psi ''(\Phi _0)\Phi _1 - \beta \partial _{zz}\Phi _1 + \beta \kappa \partial _{z}\Phi _0. \end{align}

Testing the above by ${\partial }_z \Phi _0$ and integrating from $-\infty$ to $+\infty$ produces

\begin{align*} \int _{-\infty }^{+ \infty } \Lambda _0(t,s) \partial _{z}\Phi _0(z)\,\text{d}z = \beta \int _{-\infty }^{+ \infty } \big ( \partial _{z} (\psi '(\Phi _0(z))) \Phi _1 - \partial _{zz}\Phi _1 \partial _{z}\Phi _0(z) + \kappa |\partial _{z}\Phi _0(z)|^2 \big ) \,\text{d}z. \end{align*}

Using (4.4) and (4.5) for $\Phi _0$ , integrating by parts and using $\psi '(\pm 1)=0$ , we infer

\begin{align*} & \int _{-\infty }^{+ \infty } \partial _{z} (\psi '(\Phi _0(z)))\Phi _1 - \partial _{zz}\Phi _1 \partial _{z}\Phi _0(z) \,\text{d}z \\ & \quad = \big [\psi '(\Phi _0)\Phi _1 - \partial _{z}\Phi _1 \partial _{z}\Phi _0\big ]^{+\infty }_{-\infty } -\int _{-\infty }^{+ \infty } \partial _z \Phi _1 (\psi '(\Phi _0(z))-\partial _{zz}\Phi _0(z)) \,\text{d}z=0, \end{align*}

as both terms on the right-hand side vanish, whence, recalling (4.8) and (4.9), this entails that

\begin{equation*} 2 \mu _0 = \gamma \beta \kappa . \end{equation*}

In order to utilise the equations to order $\varepsilon ^{-1}$ , we now exploit the preliminary assumption $r_c=1$ . From (2.1a), we obtain

\begin{align*} -\mathcal{V} \partial _{z}\Phi _0 = \partial _{z}(m(\Phi _0)\partial _z \Lambda _1) -[ {K_-} G_2(\Phi _0)-{K_+} G_2(\!-\!\Phi _0)]+L G_4(\Phi _0) . \end{align*}

Integrating from $-\infty$ to $+\infty$ and using (4.6) for $\Lambda _1$ leads to

\begin{align*} -2 \mathcal{V} & = [m({{\varphi }}_0)\nabla \mu _0]^+_- \cdot {\boldsymbol \nu } + \int _{- \infty }^{{+\infty }} {K_+} G_2(\!-\!\Phi _0(z)) - {K_-} G_2(\Phi _0(z)) + L G_4(\Phi _0(z))\,\text{d}z. \end{align*}

As $\Phi _0$ is a fixed function, the integral on the right-hand side yields a constant $S_I$ depending on $K_\pm$ , $L$ and the double-well potential $\psi$ . By $-\Phi _0(z)=\Phi _0(\!-\!z)$ , it holds that

(4.12) \begin{align} S_I & \;:\!=\; \int _{- \infty }^{{+\infty }} {K_+} G_2(\!-\!\Phi _0(z)) - {K_-} G_2(\Phi _0(z)) \,\text{d}z + L \int _{- \infty }^{{+\infty }} G_4(\Phi _0(z))\,\text{d}z \\ \notag &= {({K_+}- {K_-})} \int _{- \infty }^{{+\infty }} G_2(\Phi _0(z)) \,\text{d}z + L \int _{- \infty }^{{+\infty }} G_4(\Phi _0(z))\,\text{d}z, \end{align}

where, upon setting $\widetilde c \;:\!=\; \frac {1}{2 \sqrt {\psi ''(\!-\!1)}}$ , we realise with the help of the equipartition of energy (4.8) that

(4.13) \begin{equation} \begin{aligned} & \int _{- \infty }^{{+\infty }} G_2(\Phi _0(z)) \,\text{d}z = -\widetilde c \int _{-\infty }^{{+\infty }} (\Phi _0(z)-1)\sqrt {2 \psi (\Phi _0(z))}\,\text{d}z \\ & \quad =-\widetilde c \int _{-\infty }^{{+\infty }} (\Phi _0(z)-1)\partial _{z}\Phi _0(z)\,\text{d}z = -\widetilde c \int _{-1}^1 (s-1)\,\text{d}s = 2 \widetilde c. \end{aligned} \end{equation}

In addition, we notice that

\begin{equation*} G_4(\Phi _0(z))=2 \psi (\Phi _0(z))=\sqrt {2 \psi (\Phi _0(z))} \sqrt {2 \psi (\Phi _0(z))} =\sqrt {2 \psi (\Phi _0(z))}\partial _{z} \Phi _0(z) \end{equation*}

and compute

(4.14) \begin{align} \int _{- \infty }^{{+\infty }} G_4(\Phi _0(z))\,\text{d}z= \int _{- \infty }^{{+\infty }} \sqrt {2 \psi (\Phi _0(z))} \partial _{z}\Phi _0(z) {\,\text{d}z}= \int _{-1}^1 \sqrt {2 \psi (s)} \,\text{d}s=\gamma . \end{align}

Thus, recalling the definition of $\widetilde c$ , upon substituting (4.13) and (4.14) into (4.12), we find that

(4.15) \begin{align} S_I= \frac {{{K_+}- {K_-}}}{\sqrt {\psi ''(\!-\!1)}}+\gamma L. \end{align}

This concludes our analysis concerning the sharp interface system (2.11). Let us point out that, for the quartic potential (2.7), it holds that $\widetilde c = \frac {\sqrt 2}4$ , and so

(4.16) \begin{align} S_I= \frac {\sqrt 2}2 ({{K_+}- {K_-}})+\frac {2\sqrt {2}}3 L. \end{align}

We remark that in the case that $r_c\in (0,1)$ the term $S_I$ can be computed if we insert $\Phi _0$ in the definition of the source term, see (2.2), and integrate from $-\infty$ to $\infty$ . This would then replace the integrals in (4.12).

Combining the calculations above, we have demonstrated that as $\varepsilon \to 0$ , the phase-field system (2.1) formally converges to the limit described by the free-boundary problem (2.11).

5.1. Setting

For $d \in \{2,3\}$ , a simple solution can be computed in a planar geometry. In particular, we will later use planar solutions to validate the asymptotic analysis from the previous section. More precisely, we will compare numerical computations for the phase field model with planar solutions of the sharp interface limit. We now consider the free boundary problem (2.11) in the special domain

\begin{equation*} \Omega =(0,\mathcal{L})\times (0,\widetilde {\mathcal L})^{d-1}, \end{equation*}

for $\mathcal{L}, \widetilde {\mathcal L}\gt 0$ and look for a planar solution under the geometry

\begin{align*} \Omega _+(t)=(0,q(t)) \times (0, \widetilde {\mathcal L})^{d-1}, \quad \Omega _-(t)=(q(t),\mathcal{L}) \times (0, \widetilde {\mathcal L})^{d-1}, \end{align*}

where $q(t)$ encodes the location of the moving interface $\Sigma _0$ . In addition, we require a $90^\circ$ degree boundary condition at points where the interface meets the external boundary. For $x\in {\mathbb{R}}^d$ , we write $x=(z,\widehat {x})$ with $z\in {\mathbb{R}}$ and $\widehat {x}\in {\mathbb{R}}^{d-1}$ . As $\boldsymbol{\nu }=(\!-\!1,{\boldsymbol{0}})^\top$ , we get

\begin{equation*} \mathcal{V} = -\frac {dq(t)}{dt} {=-\dot q}{,} \end{equation*}

with the dot denoting the time derivative. For this planar setting, we make the ansatz

\begin{equation*} { \mu _\pm (t, x )= } \mu _\pm (t, (z,\widehat {x}) ) =\widehat {\mu }_\pm (t,z) , \quad x=(z,\widehat {x}) \in (0, \mathcal{L})\times (0, \widetilde {\mathcal L})^{d-1}. \end{equation*}

On the interface, we obtain

\begin{equation*} \nabla \mu _\pm {(t,x)}\cdot \boldsymbol{\nu } =- \widehat {\mu }_\pm ^\prime (t,z), \end{equation*}

where prime denotes the partial derivative with respect to $z$ . In what follows, we drop the hat-notation for convenience. Due to the planar setting, we have $\kappa =0$ and hence (2.11d) and (2.11c) can be replaced by

(5.1) \begin{equation} \mu _+(t,q(t)) =\mu _-(t,q(t)) =0, \end{equation}

while (2.11e) can be reformulated as

(5.2) \begin{equation} 2 {\dot q} = - [m \mu ']^+_- +S_I. \end{equation}

Equation (2.11f) can be replaced by

(5.3) \begin{equation} \mu _+^\prime (t,0)=0, \quad \mu _-^\prime (t,\mathcal{L})=0. \end{equation}

In addition, (2.11a) and (2.11b) become

(5.4) \begin{align} -m_+ \mu ''_+ &=S_+ -\rho _+ \mu _+ \qquad \text{for } z\in (0,q(t)) , \end{align}

(5.5) \begin{align} -m_- \mu ''_- &=S_- -\rho _-\mu _- \qquad \text{for } z\in (q(t),\mathcal{L} ) . \end{align}

We remark that in the setting where $S_{\pm } = 0$ and $\rho _{\pm } = 0$ , similar stability analysis for planar solutions can be found in [Reference Gurtin27, Reference Mullins and Sekerka31].

5.2. Solution formula and evolution equation for the interface position

From now on, we restrict ourselves to the case $\rho _\pm \gt 0$ which is relevant for applications, as shown in [Reference Zwicker, Seyboldt, Weber, Hyman and Jülicher38]. We obtain the following solutions

(5.6) \begin{align} \mu _+(t,z)&= d_+ \bigg (1-\frac {\cosh {({\lambda }_+ z)}}{\cosh {({\lambda }_+ q{(t)})}} \bigg ), \end{align}

(5.7) \begin{align} \mu _-(t,z)&= d_- \bigg (1-\frac {\cosh {({\lambda }_- (\mathcal{L}-z))}}{\cosh {({\lambda }_- (\mathcal{L}-q{(t)}))}} \bigg ){,} \end{align}

where the constants $d_{\pm }$ and $\lambda _\pm$ are defined as

\begin{equation*} d_{\pm } \;:\!=\; \frac {S_{\pm }}{\rho _{\pm }} = \beta \frac {S_{\pm } \psi ''(\pm 1)}{K_{\pm }}, \quad {\lambda _\pm } \;:\!=\; \sqrt {\frac {\rho _{\pm }}{m_{\pm }}}. \end{equation*}

It can be readily verified that (5.6) and (5.7) solve the ordinary differential (5.4) and (5.5), and fulfil the boundary conditions (5.1) and (5.3).

Meanwhile, the evolution (5.2) for the interface position ${q=}q(t)$ becomes

(5.8) \begin{equation} \begin{aligned} {\dot q} &= \frac {1}{2} \Big (\!-\!m_+ \mu '_+ (q) + m_- \mu '_- (q) +S_I \Big ) \\ &= \frac {1}{2} \Big ( d_+m_+ {\lambda }_+ \tanh {( {\lambda }_+ q)}+ d_-m_- {\lambda }_- \tanh {( {\lambda }_- (\mathcal{L}-q))} +S_I \Big ) \\ & \;=\!:\; \mathcal{H}(q). \end{aligned} \end{equation}

We notice that

\begin{align*} \mathcal{H} (0) = \frac 12[d_-m_- {\lambda }_- \tanh ({\lambda }_- \mathcal{L}) +S_I], \quad \mathcal{H} (\mathcal{L}) = \frac 12[d_+m_+ {\lambda }_+ \tanh ({\lambda }_+ \mathcal{L} ) +S_I]. \end{align*}

In the case where $S_I = 0$ , by having $d_+ \lt 0$ and since the other parameters $m_{\pm }$ , $\lambda _\pm$ and $d_-$ are positive, it holds that $\mathcal{H}(0) \gt 0$ , $\mathcal{H}(\mathcal{L}) \lt 0$ and $\mathcal{H}'(q) \lt 0$ . Due to the continuity of $\mathcal{H}(q)$ , we are guaranteed the existence of exactly one root $q^*$ where $\mathcal{H}(q^*) = 0$ . For example, setting $S_+ = -1$ and all other parameters equal to 1, we find that $q^* = \frac {\mathcal{L}}{2}$ is a root of $\mathcal{H}$ . Note that when $S_I \neq 0$ it is possible that roots of $\mathcal{H}$ may not exist at all. However, after fixing parameters $d_{\pm }$ and $\lambda _\pm$ (hence fixing also $K_{\pm })$ , one can adjust the parameter $L$ in $S_I$ to help ensure the existence of a root for $\mathcal{H}(q)$ .

5.3. Linear stability of planar solutions

We now aim to analyse the stability of the planar solutions, with their position denoted by $q^*$ and the corresponding chemical potentials denoted by $\mu _\pm ^*$ . Specifically, we consider a perturbed interface of the form $w \;:\!=\; q^* + \epsilon Y$ , where $0 \lt \epsilon \lt 1$ and $Y = Y(t, \widehat {x})$ . The idea is to start from the stationary front characterised by $q^*$ , as identified above, and proceed with a stability analysis around this equilibrium. Thus, we define perturbed domains as

\begin{align*} \Omega _{Y,\epsilon }^+ & \;:\!=\; \{ x \in \Omega \,\, : \,\, {z} \lt w(t, \widehat {x}) \}, {\quad \text{and} \quad } \Omega _{Y,\epsilon }^- \;:\!=\; \{ x \in \Omega \,\, : \,\, {z} \gt w(t, \widehat {x}) \}. \end{align*}

We make the ansatz

\begin{equation*} \mu _\pm (t,x) = \mu _\pm ^* ({z}) + \epsilon u_\pm (t,x), \end{equation*}

and demand that those solve the free boundary problem on the perturbed domains, given by

(5.9a) \begin{align} &- m_+ \Delta (\mu _{+}^* + \epsilon u_+) = S_+ - \rho _+ (\mu _+^* + \epsilon u_+) \quad && \text{ in } \Omega _{Y,\epsilon }^+, \end{align}

(5.9b) \begin{align} &- m_- \Delta (\mu _{-}^* + \epsilon u_-) = S_- - \rho _- (\mu _-^* + \epsilon u_-) \quad & & \text{ in } \Omega _{Y,\epsilon }^-, \end{align}

(5.9c) \begin{align} &\mu _{+}^* + \epsilon u_+ = \mu _{-}^* + \epsilon u_- \quad && \text{ on } \{{z} = w\}, \end{align}

(5.9d) \begin{align} & 2 (\mu _{\pm }^* + \epsilon u_{\pm }) = \gamma \beta \kappa \quad && \text{ on } \{{z} = w\}, \end{align}

(5.9e) \begin{align} & - 2 \mathcal{V} = [m \nabla (\mu ^* + \epsilon u)]_{-}^+ \cdot \boldsymbol{\nu } + S_I \quad && \text{ on } \{{z} = w\}, \end{align}

(5.9f) \begin{align} & (\mu _{+}^* + \epsilon u_+)'(t, 0) = 0, \quad (\mu _-^* + \epsilon u_-)'(t, \mathcal{L}) = 0. && \end{align}

In the literature, see e.g. [Reference Eck, Garcke and Knabner17, Reference Garcke, Lam, Sitka and Styles23, Reference Gurtin27, Reference Mullins and Sekerka31], $Y$ is directly chosen with a specific ansatz, which we specify later. We then linearise the above equations about the original interface $\{{z} = q^*\}$ , using

(5.10) \begin{align} \mu _{\pm }^*(w) & = \mu _{\pm }^*(q^*) + (\mu _{\pm }^*)'\vert _{{z} = q^*} (w - q^*) + \text{ h.o.t.} = 0 + \epsilon (\mu _{\pm }^*)' \vert _{{z} = q^*} Y + \text{ h.o.t.}, \\ \notag (\mu _{\pm }^*)'(w) & = (\mu _{\pm }^*)'(q^*) + (\mu _{\pm }^*)''\vert _{{z} = q^*} (w - q^*) + \text{ h.o.t.} = (\mu _{\pm }^*)'(q^*) + \epsilon (\mu _{\pm }^*)'' \vert _{{z} = q^*} Y + \text{ h.o.t.}, \end{align}

as well as (5.1)–(5.5) to simplify. For instance, consider (5.9a) where upon linearisation about $\{z = q^*\}$ we have

\begin{equation*} - m_+ \Delta (\mu _+^* + \varepsilon u_+) = {\mathcal{S}}_+ - \rho _+(\mu _+^* + \varepsilon u_+) \quad \text{ in } \{z \lt q^*\}. \end{equation*}

Then, by (5.4), we obtain

\begin{equation*} - m_+ \Delta u_+ = - \rho _+ u_+ \quad \text{ in } \{ z \lt q^*\}. \end{equation*}

Meanwhile for the interfacial condition (5.9c), upon linearisation about $\{z = q^*\}$ we have by (5.10),

\begin{equation*} \mu _+^*(q^*) + \varepsilon (\mu _+^*)'(q^*)Y + \varepsilon u_+(q^*) + \text{ h.o.t.} = \mu _-^*(q^*) + \varepsilon (\mu _-^*)'(q^*)Y + \varepsilon u_-(q^*) + \text{ h.o.t.} \end{equation*}

leading to the relation

\begin{equation*} (\mu _+^*)'(q^*) Y + u_+ = (\mu _-^*)' (q^*) Y + u_- \quad \text{ on } \{ z = q^*\}. \end{equation*}

Hence, upon linearising (5.9) about the original interface $\{z = q^*\}$ , we obtain the following system for $Y$ and $u_{\pm }$ :

(5.11a) \begin{align} &- m_+ \Delta u_+ = - \rho _+ u_+ \quad && \text{ in } \{{z} \lt q^*\}, \end{align}

(5.11b) \begin{align} & - m_- \Delta u_- = - \rho _- u_- \quad && \text{ in } \{{z} \gt q^*\}, \end{align}

(5.11c) \begin{align} & (\mu _{+}^*)' \vert _{{z} = q^*} Y + u_+ = (\mu _{-}^*)' \vert _{{z} = q^*} Y + u_- \quad && \text{ on } \{{z} = q^*\}, \end{align}

(5.11d) \begin{align} &2( (\mu _{\pm }^*)' \vert _{{z} = q^*} Y + u_{\pm }) =- \gamma \beta \Delta _{\widehat {x}} Y \quad && \text{ on } \{ {z} = q^*\}, \end{align}

(5.11e) \begin{align} &2 {\partial }_tY = -m_+ ((\mu _+^*)'' \vert _{{z} = q^*} Y + u_+') + m_- ((\mu _-^*)'' \vert _{{z} = q^*} Y +u_-') \quad && \text{ on } \{ {z} = q^* \}, \end{align}

(5.11f) \begin{align} & (u_+)'(t, 0) = 0 , \quad (u_-)'(t, \mathcal{L}) = 0, && \end{align}

where the linearisation of the mean curvature operator gives

\begin{equation*} \kappa (t, \widehat {x}) \approx -\Delta _{\widehat {x}} (\epsilon Y(t, {\widehat {x}})){,} \end{equation*}

with $\Delta _{\widehat {x}}$ denoting the Laplace operator with respect to the $(d-1)$ -dimensional coordinates $\widehat {x}$ . We make the ansatz

(5.12) \begin{align} Y(t,\widehat {x}) = \delta (t) W(\widehat {x}) \end{align}

where $W$ as an eigenfunction of the $\Delta _{\widehat {x}}$ -operator with Neumann boundary conditions such that for $\widehat {\boldsymbol{\ell }} = (\ell _2, \ldots , \ell _d) \in {\mathbb{N}}_0^{d-1}$ ,

\begin{equation*} \begin{cases} \displaystyle \Delta _{\widehat {x}} W = \frac {\zeta _{\widehat{\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2} W & \text{ in } (0, \widetilde{\mathcal{L}})^{d-1},\\ \partial _{{\boldsymbol{n}}} W = 0 & \text{ on } \partial (0, \widetilde{\mathcal{L}})^{d-1}, \end{cases} \end{equation*}

and $\frac {\zeta _{\widehat {\boldsymbol{\ell}}, d}}{(\widetilde{\mathcal{L}})^2}$ serves as a corresponding eigenvalue. This leads to perturbed interfaces that fulfil a $90^\circ$ angle condition at points where the interface intersects the outer boundary. A possible eigenfunction is

\begin{equation*} W(x_2, \ldots , x_d) = \cos \Big ( \frac {\pi }{\widetilde{\mathcal{L}}} \ell _2 x_2 \Big ) \times \cdots \times \cos \Big ( \frac {\pi }{\widetilde{\mathcal{L}}} \ell _d x_d \Big ) \end{equation*}

with

\begin{equation*} \zeta _{\widehat {\boldsymbol \ell }, d} = - \pi ^2 (\ell _2^2 + \cdots + \ell _d^2). \end{equation*}

Furthermore, it holds that

\begin{equation*} \displaystyle \Delta _{\widehat {x}} Y = y(t) \Delta _{\widehat {x}} W(\widehat {x}) = y(t) \frac {\zeta _{\widehat {\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2} W(\widehat {x}) = \frac {\zeta _{\widehat {\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2} Y. \end{equation*}

Similarly, we consider the ansatz

\begin{equation*} u_{\pm }(t,x) = \delta (t) v_{\pm }(z) W(\widehat {x}), \end{equation*}

so that (5.11a)–(5.11b) reduces to (where we recall prime denotes the partial derivative with respect to $z$ )

\begin{equation*} - m_{\pm } \delta (t) (v_{\pm }'' W + v_{\pm } \Delta _{\widehat {x}} W) = - \rho _{\pm } \delta (t) v_{\pm } W \end{equation*}

and this yields the equation for $v_{\pm }$ to be solved

\begin{equation*} v_{\pm }'' = \Big ({\lambda _\pm ^2} - \frac {\zeta _{\widehat {\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2} \Big ) v_{\pm } = (\Gamma _{\pm }^{\widehat {\boldsymbol \ell }})^2 v_{\pm }, \end{equation*}

where we recall ${\lambda _\pm ^2} = \frac {\rho _{\pm }}{m_{\pm }}$ and set $\Gamma _{\pm }^{\widehat {\boldsymbol \ell }}= \sqrt {{\lambda _\pm ^2} - \frac {\zeta _{\widehat {\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2}}$ . We then consider

\begin{align*} v_{+}(z) = a_+ \cosh (\Gamma _{+}^{\widehat {\boldsymbol \ell }} z), \quad v_{-}(z) = a_- \cosh (\Gamma _{-}^{\widehat {\boldsymbol \ell }}(\mathcal{L} - z)), \end{align*}

for unknown coefficients $a_+$ and $a_-$ to be determined. We now need to ensure that the boundary conditions (5.11c)–(5.11f) are fulfilled. On recalling (5.6) and (5.7), we see that

\begin{equation*} \begin{alignedat}{3} (\mu _+^*)'(q^*) & = - d_+ {\lambda }_+ \tanh ({\lambda }_+ q^*), \quad && (\mu _+^*)''(q^*)&& = - d_+ {\lambda }_+ ^2, \\ (\mu _-^*)'(q^*) & = d_- {\lambda }_- \tanh ({\lambda }_- ( \mathcal{L} - q^*)), \quad && (\mu _-^*)''(q^*) &&= -d_- {\lambda }_- ^2, \end{alignedat} \end{equation*}

and so (5.11c) and (5.11d) become

\begin{align*} & - d_+ {\lambda }_+ \tanh ({\lambda }_+ q^*) + a_+ \cosh (\Gamma _{+}^{\widehat {\boldsymbol \ell }} q^*) \\ & \quad = d_- {\lambda }_- \tanh ({\lambda }_- (\mathcal{L} - q^*)) + a_- \cosh (\Gamma _{-}^{\widehat {\boldsymbol \ell }} (\mathcal{L} - q^*)) = - \frac {\gamma \beta }{2} \frac {\zeta _{\widehat {\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2}. \end{align*}

Hence, we derive the following:

\begin{align*} a_+ (q^*)& = \frac {1}{\cosh (\Gamma _{+}^{\widehat {\boldsymbol \ell }} q^*)} \Big ( d_+ {\lambda }_+ \tanh ({\lambda }_+ q^*) - \frac {\gamma \beta }{2} \frac {\zeta _{\widehat {\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2} \Big ), \\ a_-(q^*) & = -\frac {1}{\cosh (\Gamma _{-}^{\widehat {\boldsymbol \ell }} (\mathcal{L} - q^*))} \Big ( d_- {\lambda }_- \tanh ({\lambda }_- (\mathcal{L} - q^*)) + \frac {\gamma \beta }{2} \frac {\zeta _{\widehat {\boldsymbol \ell }, d}}{(\widetilde{\mathcal{L}})^2} \Big ). \end{align*}

Meanwhile, for the evolution (5.11e), we have

\begin{align*} 2 {\partial }_t Y & = m_+ d_+ {\lambda }_+ ^2 Y - m_- d_- {\lambda }_- ^2 Y \\ & \quad - m_+ a_+(q^*) \Gamma _+^{\widehat {\boldsymbol \ell }} \sinh (\Gamma _+^{\widehat {\boldsymbol \ell }} q^*) Y - m_- a_-(q^*) \Gamma _{-}^{\widehat {\boldsymbol \ell }} \sinh (\Gamma _-^{\widehat {\boldsymbol \ell }} (\mathcal{L} - q^*)) Y. \end{align*}

Using the explicit form of $Y$ in (5.12), the above equality can be reduced to an ordinary differential equation for the perturbation magnitude $\delta (t)$ :

\begin{align*} 2\dot {\delta } & = (S_+ - S_- -m_+ a_+(q^*) \Gamma _+^{\widehat {\boldsymbol \ell }} \sinh (\Gamma _+^{\widehat {\boldsymbol \ell }} q^*) - m_- a_-(q^*) \Gamma _{-}^{\widehat {\boldsymbol \ell }} \sinh (\Gamma _-^{\widehat {\boldsymbol \ell }} (\mathcal{L} - q^*)) ) \delta , \end{align*}

where we recall the definition $d_{\pm } = \frac {S_{\pm }}{\rho _{\pm }} = \frac {S_{\pm }}{\lambda _\pm ^2 m_{\pm }}$ . The amplification factor is the prefactor coefficient on the right-hand side of the above identity

(5.13) \begin{align} S_+ - S_- -m_+ a_+(q^*) \Gamma _+^{\widehat {\boldsymbol \ell }} \sinh (\Gamma _+^{\widehat {\boldsymbol \ell }} q^*) - m_- a_-(q^*) \Gamma _{-}^{\widehat {\boldsymbol \ell }} \sinh (\Gamma _-^{\widehat {\boldsymbol \ell }} (\mathcal{L} - q^*)). \end{align}

For the rest of this section, we consider specific parameters ensuring the amplification factor (5.13) is positive so as to induce the development of instabilities from interface perturbations. We fix $\widehat {\boldsymbol \ell } = (\ell _2,\ldots ,\ell _d)$ , and

\begin{equation*} S_+ = -1, \quad S_{-} = m_{\pm } = \rho _{\pm } = 1{,} \end{equation*}

so that the stationary state is $q^* = \frac {\mathcal{L}}{2}$ , along with $d_{+} = -1$ , $d_{-} = 1$ , ${\lambda _\pm } = 1$ ,

\begin{align*} & \Gamma _{\pm }^{\widehat {\boldsymbol \ell }} = \sqrt {1 + \frac {|\widehat {\boldsymbol \ell }|^2 \pi ^2}{\widetilde{\mathcal{L}}^2}}, \quad a_+(q^*) = \frac {-\tanh ( \frac {\mathcal{L}}{2}) + \frac {\gamma \beta }{2} |\widehat {\boldsymbol \ell }|^2 \pi ^2/\widetilde {\mathcal{L}}^2 } {\cosh {\Big (} \frac {\mathcal{L}}{2} \sqrt {1 + |\widehat {\boldsymbol \ell }|^2 \pi ^2/\widetilde{\mathcal{L}}^2}{\Big )}}, \quad \\ & a_-(q^*) = \frac {-\tanh ( \frac {\mathcal{L}}{2}) +\frac {\gamma \beta }{2} |\widehat {\boldsymbol \ell }|^2 \pi ^2/\widetilde {\mathcal{L}}^2}{\cosh {\Big (} \frac {\mathcal{L}}{2} \sqrt {1 + |\widehat {\boldsymbol \ell }|^2 \pi ^2/\widetilde{\mathcal{L}}^2}{\Big )}}{,} \end{align*}

and the amplification factor (5.13) reduces in this setting to

(5.14) \begin{align} -2 + \sqrt {1 + \frac {|\widehat {\boldsymbol \ell }|^2 \pi ^2}{\widetilde {\mathcal{L}}^2}} \tanh \left ( \frac {\mathcal{L}}{2}\sqrt {1 + \frac {|\widehat {\boldsymbol \ell }|^2 \pi ^2}{\widetilde {\mathcal{L}}^2}} \right ) \left (2 \tanh \left ( \frac {\mathcal{L}}{2} \right ) - {\gamma \beta } \frac {|\widehat {\boldsymbol \ell }|^2 \pi ^2}{\widetilde{\mathcal{L}}^2} \right ). \end{align}

In the case $\widehat {\boldsymbol \ell } ={(0,\ldots ,0)}$ , this corresponds to interface perturbations of the form $w = q^* + \epsilon$ , i.e., translational perturbations. We immediately see that the amplification factor becomes

\begin{equation*} -2 + 2\tanh ^2 \Big ( \frac {\mathcal{L}}{2} \Big ) = -2 \mathrm{sech}^2 \Big ( \frac {\mathcal{L}}{2} \Big ) \lt 0. \end{equation*}

Thus, we obtain stability with respect to translational perturbations. For $|\widehat {\boldsymbol \ell }| \gt 0$ , we obtain interface perturbations of the form $w = q^* + \epsilon \cos (\pi \ell _2 x_2 / \widetilde {\mathcal{L}})\times {\cdots }\times \cos (\pi \ell _d x_d / \widetilde {\mathcal{L}})$ . We see that (5.14), for a given perturbation whose wave length is related to $\widehat {\boldsymbol \ell }$ , can be positive if $\beta \lt \beta _{\mathrm{crit}}(\widehat {\boldsymbol \ell })$ , where

(5.15) \begin{align} \beta _{\mathrm{crit}} (\widehat {\boldsymbol \ell })= \frac {2}{\gamma } \frac {\widetilde {\mathcal{L}}^2}{|\widehat {\boldsymbol \ell }|^2 \pi ^2} \left ( \tanh \Big ( \frac {\mathcal{L}}{2} \Big ) - \frac {1}{\sqrt {1 + |\widehat {\boldsymbol \ell }|^2 \pi ^2/\widetilde{\mathcal{L}}^2} \tanh \big ( \frac {\mathcal{L}}{2} \sqrt {1 + |\widehat {\boldsymbol \ell }|^2 \pi ^2 / \widetilde{\mathcal{L}}^2} \big )} \right ). \end{align}

We remark that in this setting $\beta$ is the modified capillary length $c_l$ introduced in (2.13) in section 2.4, i.e., we obtain instability in cases where the modified capillary length is small enough. As $|\widehat {\boldsymbol \ell }|^2 = (\ell _2)^2+\cdots +(\ell _d)^2$ with ${\widehat {\boldsymbol \ell }=(} \ell _2,\ldots ,\ell _d) \in \mathbb{N}_0^{d-1}$ , we obtain that the possible perturbations lead to amplification factors via sum of squares, namely:

\begin{equation*} |\widehat {\boldsymbol \ell }|^2 \in \begin{cases} \{0, 1, 4, 9, 16, \ldots \} & \text{ if } d = 2, \\ \{0, 1, 2, 4, 5, 8, 9, 10, \ldots \} & \text{ if } d = 3. \end{cases} \end{equation*}

6.1. Finite element method

We assume that $\Omega$ is a polyhedral domain and let $\mathcal{T}_h$ be a regular triangulation of $\Omega$ into disjoint open simplices. Associated with $\mathcal{T}_h$ is the piecewise linear finite element space

\begin{equation*} S^h = \{ \zeta \in C^0(\overline {\Omega }) \, : \, \zeta \vert _{o} \in P_1(o) \quad \forall o \in \mathcal{T}_h\}, \end{equation*}

where we denote by $P_1(o)$ the set of all affine linear functions on $o$ , see [Reference Ciarlet13]. Let $(\cdot ,\cdot )_{L^2}^h$ be the usual mass lumped $L^2$ -inner product on $\Omega$ associated with $\mathcal{T}_h$ , and let $\pi ^h \;:\; C^0(\overline {\Omega }) \to S^h$ be the standard interpolation operator. In addition, let $\tau$ denote a chosen uniform time step size. Then our finite element approximation of (2.1) is given as follows. Let $\phi _h^0 \in S^h$ , e.g., $\phi _h^0 = \pi ^h \varphi _0$ if $\varphi _0 \in C^0(\overline \Omega )$ . Then, for $n \geq 0$ , find $(\phi _h^{n+1}, \mu _h^{n+1}) \in S^h \times S^h$ such that

(6.1a) \begin{align} \frac {1}{\tau }(\phi _h^{n+1}-\phi _h^n, \chi _h)_{L^2}^h + (m(\phi _h^n) \nabla \mu _h^{n+1}, \nabla \chi _h)_{L^2}^h = (S_\varepsilon (\phi _h^n), \chi _h)_{L^2}^h \quad & \forall \chi _h\in S^h, \end{align}

(6.1b) \begin{align} \beta \varepsilon (\nabla \phi _h^{n+1}, \nabla \eta _h)_{L^2} + \frac {\beta }{\varepsilon } (\psi '(\phi _h^{n+1}), \eta _h)_{L^2}^h - (\mu _h^{n+1}, \eta _h)_{L^2}^h = 0 \quad &\forall \eta _h\in S^h. \end{align}

We implemented the scheme (6.1) with the help of the finite element toolbox ALBERTA, see [Reference Schmidt and Siebert35]. To increase computational efficiency, we employ adaptive meshes, which have a finer mesh size $h_{f}$ within the diffuse interfacial regions and a coarser mesh size $h_{c}$ away from them. In particular, we use the strategy from [Reference Barrett, Nürnberg and Styles5, Reference Baňas and Nürnberg7], and refine an element $o$ if $\eta _o = |\max _o |\phi _h^n| - 1 | \gt 0.5$ , unless it is already of size $h_f$ , and similarly coarsen an element $o$ if $\eta _o \lt 0.1$ , unless it is already of size $h_c$ . For simplicity, we assume from now on that $\Omega = \prod _{i=1}^d (0,L_i) \subset \mathbb R^d$ , with $L_1 \geq L_2 \geq { \ldots \geq } L_d$ , and then let $h_f = \frac {L_d}{N_f}$ , $h_c = \frac {L_d}{N_c}$ for two chosen integer parameters $N_f \gt N_c$ . Here, unless otherwise specified, for the computations with a phase field parameter $\varepsilon = (2^k\pi )^{-1}$ , $k \in \mathbb N$ , we choose $N_f = 8N_c = 2^{3+k} L_d$ . This ensures that the interfacial regions are accurately resolved while using a relatively coarser mesh in the pure regions. For the time discretisation, we choose $\tau =10^{-3}$ , unless stated otherwise.

The nonlinear system of equations arising at each time level of (6.1) is solved with the help of Newton’s method. The resulting linear systems at each iteration in two spatial dimensions are solved by direct factorisation using the package UMFPACK, see [Reference Davis15], and in three spatial dimensions with a $V$ -cycle multigrid solver using a block Gauss–Seidel smoother.

For the initial data $\varphi _0$ , we in general choose a diffuse interface representation of a desired sharp interface, with signed distance function $d_0 \;:\; \overline \Omega \to \mathbb R$ . In particular, unless otherwise stated, we let

(6.2) \begin{equation} \varphi _0(x) = \tanh \bigg ( \frac {d_0(x)}{\varepsilon \sqrt {2}} \bigg) . \end{equation}

For the definition of $S_2$ in (2.2), recall (2.4), we need to specify $K_\pm$ and $L$ . For convenience, for the numerical simulations to follow, we will define the relations $\rho _\pm = \frac {K_\pm }{2\beta }$ , from which $K_\pm$ can then be inferred. The mobility function $m$ is chosen as

\begin{align*} m(s) & = \begin{cases} m_+ & \text{ for } s \geq 1, \\ m_+ \frac {1+s}{2} + m_- \frac {1-s}{2} & \text{ for } |s| \leq 1, \\ m_- & \text{ for } s \leq -1, \end{cases} \end{align*}

where $m_{\pm } \gt 0$ are chosen fixed parameters. Unless otherwise stated, we set $m_\pm =\rho _\pm =r_c=1$ and $L=0$ throughout.

Figure 1. Convergence experiment for a moving front in $(0,1)^2$ . We compare the true solution $q$ of the sharp interface problem with the discrete approximations $q_h$ of the Cahn–Hilliard equation for $\varepsilon =(2^k\pi )^{-1}$ , $k=2,4,6$ .

6.2. Numerical computations: Planar solutions and their stability

We begin with a convergence experiment for the exact solution of a moving flat vertical front in $\Omega = (0,1)^2$ , whose $x$ -position $q(t)$ satisfies the ODE (5.8). To this end, we solve the ODE numerically and compare the obtained approximation of $q(t)$ with the position of the diffuse front for solutions of our finite element scheme (6.1) for decreasing values of $\varepsilon$ . In particular, for a given $\phi ^n_h$ we will compare $q(t_n)$ , where $t_n = n\tau$ , with the unique value $q_h(t_n) \in (0,1)$ such that $\phi ^n(q_h(t_n), 0) = 0$ . In particular, for the physical parameters we set $\beta =0.1$ , $S_-=4$ , $S_+=-1$ , $\rho _-=0.1$ , $\rho _+=1$ and $L=-1$ , and we compute the evolutions over the time interval $[0,1]$ . For the initial position of the flat interface, we choose $q(0) = 0.3$ . Then we calculate the evolutions for (6.1) for $\varepsilon =(2^{k}\pi )^{-1}$ , $k=2,\ldots ,6$ . A comparison between $q$ and the various $q_h$ can be seen in Figure 1, where we note an excellent agreement between the true solution $q$ and the numerical approximations $q_h$ when $\varepsilon$ is sufficiently small. In order to more closely investigate the convergence behaviour of the phase field solution to the sharp interface solution as $\varepsilon \to 0$ , we also compute the error between $q(T)$ and $q_h(T)$ at time $T=1$ and display these values in Table 1. The experimental  order of convergence suggests a quadratic convergence in $\varepsilon$ . This is backed by the fact that the first order correction $\Phi _1$ , solving (4.11), is zero in the case that the curvature $\kappa$ is zero, which holds for a flat interface, see also [Reference Garcke and Stinner26].

For completeness, we repeat the same convergence experiment also in three dimensions on the unit cube $\Omega =(0,1)^3$ . As expected, the observed results are nearly identical to the earlier two-dimensional computations, see Figure 2 and Table 2.

For stability investigations, we let $\Omega =(0,1)^2$ and set $\beta = 0.1$ , $S_\pm = \mp 8$ . This encourages growth of the $2$ -mode, as can be seen in Figure 3 for a run with $\varepsilon =\frac 1{32\pi }$ . Here, as initial data, we use a flat front at the middle of the domain, with an added perturbation of magnitude less than $0.1$ given by a sum of modes from 1 to 20 with random coefficients. Computing the amplification factors in (5.13) shows that ${\widehat {\boldsymbol \ell }}=(\ell _2)=(2)$ is the most unstable mode. This agrees with the numerical solution plotted in Figure 3, where the $2$ -mode is the one which is most amplified.

As an alternative, we compute on $\Omega = (0,4) \times (0,2)$ with $\beta =0.1$ , $S_\pm = \mp 1.5$ . For the same type of initial perturbation as before, this once again encourages growth of the $2$ -mode, see Figure 4 for a run with $\varepsilon =\frac 1{16\pi }$ . Also in this case, the linear stability analysis predicts that the $2$ -mode is most unstable.

Moreover, on the domain $\Omega = (0,4) \times (0,1)$ we compute with $\beta =0.1$ , $S_\pm = \mp 3$ . See Figure 5 for a simulation with $\varepsilon =\frac 1{32\pi }$ , where we observe the growth of the 1-mode, which is also predicted by the linear stability analysis when computing the numbers for the amplification factors in (5.13) for different values of $\widehat {\boldsymbol \ell }$ . At later times, the long and nearly horizontal interface becomes unstable for higher modes. Here, as initial data, we use a flat front at the middle of the domain, with an added perturbation of magnitude less than $0.025$ given by a sum of modes from 1 to 20 with random coefficients.

Table 1. Convergence experiment for a moving front in $(0,1)^2$ , over the time interval $[0,1]$ . We also display the experimental order of convergence (EOC)

Table 2. Convergence experiment for a moving front in $(0,1)^3$ , over the time interval $[0,1]$ . We also display the experimental order of convergence (EOC)

Figure 2. Convergence experiment for a moving front in $(0,1)^3$ . We compare the true solution $q$ of the sharp interface problem with the discrete approximations $q_h$ of the Cahn–Hilliard equation for $\varepsilon =(2^k\pi )^{-1}$ , $k=2,3,5$ .

Moreover, on the domain $\Omega = (0,2)^2$ we compute with $\beta =0.01$ , $S_\pm = \mp 1.3$ . See Figure 6 for a simulation with $\varepsilon =\frac 1{64\pi }$ , where we observe the growth of the 5-mode. Here, as initial data, we use a flat front at position $0.5$ , with an added perturbation of magnitude less than $0.025$ given by a sum of modes from 1 to 20 with random coefficients. It is worth noticing that the 5-mode is also predicted to grow the most by the amplification factors obtained in the linear stability analysis. In this case, we computed the values in (5.13) for a $q(0)$ , which leads to a planar solution which is not stationary.

Figure 3. ( $\varepsilon =\frac 1{32\pi }$ , $\Omega =(0,1)^2$ ) evolution for $\beta =0.1$ , $S_-=8$ , $S_+ = -8$ . We show the solution at times $t=0,0.1,1,2,10$ .

Figure 4. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega = (0,4) \times (0,2)$ ) evolution for $\beta =0.1$ , $S_-=1.5$ , $S_+ = -1.5$ . We show the solution at times $t=0,1,2,5,10$ .

Figure 5. ( $\varepsilon =\frac 1{32\pi }$ , $\Omega = (0,4) \times (0,1)$ ) evolution for $\beta =0.1$ , $S_-=3$ , $S_+ = -3$ . We show the solution at times $t=0,1,2,5,10$ .

Figure 6. ( $\varepsilon =\frac {1}{64\pi }$ , $\Omega = (0,2)^2$ ) evolution for $\beta =0.01$ , $S_-=1.3$ , $S_+ = 1.3$ . We show the solution at times $t=0, 1, 2, 3, 20$ .

For a three-dimensional analogue of Figure 3, we use the parameters $\beta = 0.1$ , $S_\pm = \mp 4.5$ , $m_\pm =0.2$ on the unit cube $\Omega =(0,1)^3$ . The initial perturbation of a flat interface at position $q=0.5$ is made up of a single mode with maximal magnitude $0.2$ . The evolution is shown in Figure 7.

6.3. Numerical computations: Spinodal decomposition

In this subsection, we are interested in simulations that demonstrate spinodal decomposition. To this end, we choose for the discrete initial data $\varphi _0^h$ a random function with zero mean and values inside $[-0.1, 0.1]$ . On the domain $\Omega =(0,4)^2$ we then choose the physical parameters $\beta =0.002$ , $S_-=0.25$ , $S_+ = -4$ , and let $\varepsilon = \frac 1{16\pi }$ . The simulation is shown in Figure 8. Increasing the value of $\beta$ to $0.02$ yields the results shown in Figure 9. In both cases, the numerical solution shown at the final time appears to be a steady state. In particular, the natural discrete analogue of the energy (1.1) remains constant over a period of time. We can observe a suppression of the Ostwald ripening mechanism and the stable coexistence of multiple droplets with a characteristic size, consistent with the numerical simulations in [Reference Zwicker, Seyboldt, Weber, Hyman and Jülicher38, Reference Zwicker, Hyman and Jülicher39] as well as the experimental observations in [Reference Nakashima, van Haren, André, Robu and Spruijt32] for certain protocell models.

Figure 7. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,1)^3$ ) evolution for $\beta =0.1$ , $m_\pm = 0.2$ , $S_-=4.5$ , $S_+ = -4.5$ . We show the solution at times $t=0,0.1,0.5,1,1.5$ .

Figure 8. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,4)^2$ ) evolution for $\beta =0.002$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times $t=0.1,0.5,1,2,10,100$ . In our numerical computations, the solution at the final time appears to be a steady state.

Figure 9. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,4)^2$ ) evolution for $\beta =0.02$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times $t=0.1,0.5,1,2,10,100$ . In our numerical computations, the solution at the final time appears to be a steady state.

To understand the observed behaviour at the end of the simulations shown in Figures 8 and 9 a bit better, we consider an experiment with the same physical parameters as in Figure 9, but starting from three circular initial blobs with radii $0.29$ , $0.3$ and $0.31$ . As can be seen from the evolution shown in Figure 10, the three blobs very quickly reduce in size and move apart from each other. Soon after, they settle on an arrangement that, in our numerical computations, is a steady state.

Figure 10. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,4)^2$ ) evolution for $\beta =0.02$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times $t=0,1,10,100$ . In our numerical computations, the solution at the final time appears to be a steady state.

Next, we consider some analogous simulations in 3d. For the simulations in Figures 11, 12 and 13, we let $\Omega =(0,1)^3$ and always choose $S_+=-4$ . For the remaining physical parameters, we choose $(\beta ,S_-) = (0.02, 0.25), (0.02, 1), (0.1,1)$ , respectively. For the phase field parameter, we choose $\varepsilon =\frac 1{8\pi }$ , and set $\tau =10^{-4}$ .

Figure 11. ( $\varepsilon =\frac 1{8\pi }$ , $\Omega =(0,1)^3$ ) evolution for $\beta =0.02$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times $t=0.1,0.2,0.5,1,5$ .

Figure 12. ( $\varepsilon =\frac 1{8\pi }$ , $\Omega =(0,1)^3$ ) evolution for $\beta =0.02$ , $S_-=1$ , $S_+ = -4$ . We show the solution at times $t=0.01,0.05,0.1,0.5,5$ .

Figure 13. ( $\varepsilon =\frac 1{8\pi }$ , $\Omega =(0,1)^3$ ) evolution for $\beta =0.1$ , $S_-=1$ , $S_+ = -4$ . We show the solution at times $t=0.02,0.1,0.5,1,5$ .

For a three-dimensional analogue of Figure 10, we use the parameters $\beta = 0.1$ , $S_- = 0.8$ , $S_+= -10$ , $m_-=0.2$ , $m_+=0.5$ and $\rho _\pm = 0.6$ , on the cube $\Omega =(0,4)^3$ . We start from three spherical initial blobs with radii $0.34$ , $0.35$ and $0.36$ . As can be seen from the evolution shown in Figure 14, the three blobs hardly change their size, and soon settle on an arrangement that, in our numerical computations, is a steady state. As further coarsening is eventually suppressed, these two- and three-dimensional computations show that the active Cahn–Hilliard model can suppress Ostwald ripening. This is in agreement with [Reference Zwicker, Hyman and Jülicher39], which studied the suppression of Ostwald ripening in so-called active emulsions. In particular, multiple droplets can be stable. Due to the Neumann boundary conditions, we can extend the solution obtained by reflections in space, and we then observe also the occurrence of cylinders and toroidal shapes.

Figure 14. ( $\varepsilon =\frac 1{8\pi }$ , $\Omega =(0,4)^3$ ) evolution for $\beta =0.1$ , $S_- = 0.8$ , $S_+= -10$ , $m_-=0.2$ , $m_+=0.5$ , $\rho _\pm = 0.6$ . We show the solution at times $t=0,1,5$ . In our numerical computations, the solution at the final time appears to be a steady state.

Figure 15. ( $\varepsilon =\frac 1{32\pi }$ , $\Omega =(0,2)^2$ ) evolution for $\beta =0.002$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times $t=0,0.2,1,2,5$ .

6.4. Numerical computations: Active droplets in 2d

Our first simulation for active droplets shows the possible creation of shells in 2d. We let $\Omega =(0,2)^2$ and choose $\beta =0.002$ , $S_-=0.25$ , $S_+ = -4$ for the physical parameters. The initial droplet is a disk of radius $r_0=0.4$ . For the phase field parameter, we choose $\varepsilon =\frac 1{32\pi }$ . In Figure 15, we observe the development of a stable shell, which is consistent with the liquid spherical shells observed experimentally in Figures 2 and 3 of [Reference Bergmann, Bauermann and Bartolucci8], and numerically in Figure 5 of [Reference Bauermann, Bartolucci, Boekhoven, Weber and Jülicher6]. However, if we reduce the dimension of the domain to the unit square, $\Omega =(0,1)^2$ , then the shell shape is only transient. In fact, the temporary shell evolves to a smaller, stable disk. See Figure 16 for the results.

Figure 16. ( $\varepsilon =\frac 1{32\pi }$ , $\Omega =(0,1)^2$ ) evolution for $\beta =0.002$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times $t=0,0.2,1,2,5$ .

The next simulation shows pinch-off for a slightly perturbed initially circular droplet. In fact, for the ‘radius’ of the initial droplet we choose

(6.3) \begin{equation} r_0(\theta ) = 0.25 + 0.02 \cos \bigg(2\theta - \frac {2\pi }9\bigg), \end{equation}

with $\theta \in [-\pi ,\pi ]$ denoting the angle in polar coordinates. Hence, the initial blob has its widest dimension along an axis that is tilted by $20^\circ$ compared to the $x$ -axis, with the smallest dimension along an axis perpendicular to it. Letting this initial droplet evolve in the domain $\Omega = (0,2)^2$ leads to a thinning of the droplet at the centre of the domain, and eventually to a pinch-off into two separate droplets. See Figure 17 for our numerical results for $\varepsilon =\frac 1{16\pi }$ . Running the simulation on the larger domain $\Omega = (0,8)^2$ once again shows very delicate fingering patterns appear, but no pinch-off occurs. See Figure 18.

Figure 17. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,2)^2$ ) evolution for $\beta =0.02$ , $S_-=1$ , $S_+ = -4$ . We show the solution at times $t=0,2,4,6,8,10$ .

Figure 18. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,8)^2$ ) evolution for $\beta =0.02$ , $S_-=1$ , $S_+ = -4$ . We show the solution at times $t=0,2,4,6,8,10$ .

Figure 19. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,2)^3$ ) evolution for $\beta =0.002$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times $t=0,0.3,1,2,5$ .

6.5. Numerical computations: Active droplets in 3d

In our first simulations for active droplets in 3d, we attempt to create analogous shell structures as in Figure 15. To this end, we let $\Omega =(0,2)^3$ and choose the same physical parameters $\beta =0.002$ , $S_-=0.25$ , and $S_+ = -4$ . The initial droplet is a ball of radius $r_0=0.4$ . For the phase field parameter, we choose $\varepsilon =\frac 1{16\pi }$ . See Figure 19, where we observe the creation of a shell, which eventually changes into a single ball again. When we increase the radius of the initial ball to $r_0=0.6$ , then the ensuing evolution is more intricate, see Figure 20. At first, two concentric shells appear, with the inner shell merging into a ball after some time. Then the inner ball disappears, leaving just a thin outer shell, which starts to become thinner and thinner, and which eventually fragments into several much smaller blobs. These blobs become spherical and then continue to move away from each other, slowly increasing in size.

Figure 20. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega =(0,2)^3$ ) evolution for $\beta =0.002$ , $S_-=0.25$ , $S_+ = -4$ . We show the solution at times t = 0, 0.3, 0.5, 1, 1.5, 1.6, 2, 5, where the visualisations in the two rows differ.

For the remaining experiments in this subsection, we use $r_c = 0.5$ , recall (2.3), and consider the domain $\Omega = (0,8)^3$ . In our first experiment, we choose for the initial droplet a rounded cylinder of total dimension $0.8\times 0.4\times 0.4$ . In particular, it is made up of two half-balls of radius $0.2$ , which are connected with a cylinder of radius $0.2$ and height $0.4$ . For the physical parameters we choose $\beta =0.1$ , $S_-=0.8$ , $S_+ = -10$ , $m_-=0.2$ , $m_+=0.5$ , $\rho _-=0.2$ , $\rho _+=0.1$ , $L=-1$ , and we let $\varepsilon = \frac 1{16\pi }$ . The results of our numerical simulation can be seen in Figure 21, where we notice a primary pinch-off that then leads to several secondary pinch-offs.

Figure 21. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega = (0,8)^3$ ) evolution for $\beta =0.1$ , $S_-=0.8$ , $S_+ = -10$ , $m_-=0.2$ , $m_+=0.5$ , $\rho _-=0.2$ , $\rho _+=0.1$ , $L=-1$ . We show the solution at times $t=0,1,1.3,1.4,2,3.5$ .

For the next experiment, we only change the aspect ratio of the initial droplet. In particular, the initial droplet now is a rounded cylinder of total dimension $0.6\times 0.4\times 0.4$ . The results of this numerical simulation can be seen in Figure 22. In Figure 23, we show the interfaces at time $t=4.4$ within $\Omega$ from three different points of view. We observe that the evolution can be very complex with many splittings into daughter droplets and other topological changes, consistent with the numerical simulation of Figure 3A of [Reference Zwicker, Seyboldt, Weber, Hyman and Jülicher38], as well as the experimental observations in Figure 9 of [Reference Slootbeek, van Haren, Smokers and Spruijt36]. Figure 23 shows that droplets have been formed in the centre, and away from the centre the evolution is still very complex. We expect that eventually more and more round droplets will form.

Figure 22. ( $\varepsilon =\frac 1{16\pi }$ , $\Omega = (0,8)^3$ ) evolution for $\beta =0.1$ , $S_-=0.8$ , $S_+ = -10$ , $m_-=0.2$ , $m_+=0.5$ , $\rho _-=0.2$ , $\rho _+=0.1$ , $L=-1$ . We show the solution at times $t=0,1.4,1.8,2,3,4$ .

Figure 23. The evolution from Figure 22 at time $t=4.4$ viewed from the front, from the side and from above.