Numerical approximation of the Stochastic Cahn-Hilliard Equation near the Sharp Interface Limit

We consider the stochastic Cahn-Hilliard equation with additive noise term $\varepsilon^\gamma g\, \dot{W}$ ($\gamma>0$) that scales with the interfacial width parameter $\varepsilon$. We verify strong error estimates for a gradient flow structure-inheriting time-implicit discretization, where $\varepsilon^{-1}$ only enters polynomially; the proof is based on higher-moment estimates for iterates, and a (discrete) spectral estimate for its deterministic counterpart. For $\gamma$ sufficiently large, convergence in probability of iterates towards the deterministic Hele-Shaw/Mullins-Sekerka problem in the sharp-interface limit $\varepsilon \rightarrow 0$ is shown. These convergence results are partly generalized to a fully discrete finite element based discretization. We complement the theoretical results by computational studies to provide practical evidence concerning the effect of noise (depending on its 'strength' $\gamma$) on the geometric evolution in the sharp-interface limit. For this purpose we compare the simulations with those from a fully discrete finite element numerical scheme for the (stochastic) Mullins-Sekerka problem. The computational results indicate that the limit for $\gamma\geq 1$ is the deterministic problem, and for $\gamma=0$ we obtain agreement with a (new) stochastic version of the Mullins-Sekerka problem.


Introduction
We consider the stochastic Cahn-Hilliard equation with additive noise du = ∆ − ε∆u + 1 ε f (u) dt + ε γ g dW in D T := (0, T ) × D (1.1a) We fix T > 0, γ > 0, and ε > 0 is a (small) interfacial width parameter. For simplicity, we assume D ⊂ R 2 to be a convex, bounded polygonal domain, with n ∈ S 2 the outer unit normal along ∂D, and W ≡ {W t ; 0 ≤ t ≤ T } to be an R-valued Wiener process on a filtered probability space (Ω, F, {F t } t , P). The function g ∈ C ∞ (D) is such that D g dx = 0 to enable conservation of mass in (1.1), and ∂ n g = 0 on ∂D. Furthermore, we assume u ε 0 ∈ H 1 , and impose D u ε 0 dx = 0, for simplicity; generalization for arbitrary mean values is straightforward.
The nonlinear drift part f in (1.1) is the derivative of the double-well potential F (u) := The particular case g ≡ 0 in (1.1) leads to the deterministic Cahn-Hilliard equation which can be interpreted as the H −1 -gradient flow of the Ginzburg-Landau free energy. It is convenient to reformulate (1.1) as where w denotes the chemical potential.
The Cahn-Hilliard equation has been derived as a phenomenological model for phase separation of binary alloys. The stochastic version of the Cahn-Hilliard equation, also known as the Cahn-Hilliard-Cook equation, has been proposed in [12,21,22]: here, the noise term is used to model effects of external fields, impurities in the alloy, or may describe thermal fluctuations or external mass supply. We also mention [18], where computational studies for (1.1) show a better agreement with experimental data in the presence of noise. For a theoretical analysis of various versions of the stochastic Cahn-Hilliard equation we refer to [8,9,13,14]. Next to its relevancy in materials sciences, (1.1) is used as an approximation to the Mullins-Sekerka/Hele-Shaw problem; by the classical result [1], the solution of the deterministic Cahn-Hilliard equation is known to converge to the solution of the Mullins-Sekerka/Hele-Shaw problem in the sharp interface limit ε ↓ 0. A partial convergence result for the stochastic Cahn-Hilliard equation (1.1) has been obtained recently in [3] for a sufficiently large exponent γ. We extend this work to eventually validate uniform convergence of iterates of the time discretization Scheme 3.1 to the sharp-interface limit of (1.1) for vanishing numerical (time-step k), and regularization (width ε) parameters: hence, the zero level set of the solution to the geometric interface of the Mullins-Sekerka problem is accurately resolved via Scheme 3.1 in the asymptotic limit.
It is well-known that an energy-preserving discretization, along with a proper balancing of numerical parameters and the interface width parameter ε, is required for accurate simulation of the deterministic Cahn-Hilliard equation; see e.g. [16]: analytically, this balancing of scales allows to circumvent a straight-forward application of Gronwall's lemma in the error analysis, which would otherwise cause a factor in a corresponding error estimate that grows exponentially in ε −1 . The present paper pursues a corresponding goal for a structurepreserving discretization of the stochastic Cahn-Hilliard equation (1.1); we identify proper discretization scales which allow a resolution of interface-driven evolutions, and thus avoid a Gronwall-type argument in the corresponding strong error analysis. This allows for practically relevant scaling scenarios of involved numerical parameters to accurately approximate solutions of (1.1) even in the asymptotic regime where ε 1.
The proof of a strong error estimate for a space-time discretization of (1.1) which causes only polynomial dependence on ε −1 in involved stability constants uses the following ideas: (a) We use the time-implicit Scheme 3.1, whose iterates inherit the basic energy bound (see Lemma 3.1, i)) from (1.1). We benefit from a weak monotonicity property of the drift operator in the proof of Lemma 3.4 to effectively handle the cubic nonlinearity in the drift part. (b) For γ > 0 sufficiently large, we view (1.1) as a stochastic perturbation of the deterministic Cahn-Hilliard equation (i.e., (1.1) with g ≡ 0), and proceed analogically also in the discrete setting. We then benefit in the proof of Lemma 3.4 from (the discrete version of) the spectral estimate (2.1) from [11,2] for the deterministic Cahn-Hilliard equation (see Lemma 3.1,v)). (c) For the deterministic setting [16], an induction argument is used on the discrete level, which addresses the cubic error term (scaled by ε −1 ) in Lemma 3.4. This argument may not be generalized in a straightforward way to the current stochastic setting where the discrete solution is a sequence of random variables allowing for (relatively) large temporal variations. For this reason we consider the propagation of errors on two complementary subsets of Ω: on the large subset Ω 2 we verify the error estimate (Lemma 3.5), while we benefit from the higher-moment estimates for iterates of Scheme 3.1 from (a) to derive a corresponding estimate on the small set Ω \ Ω 2 (see Corollary 3.7). A combination of both results then establishes our first main result: a strong error estimate for the numerical approximation of the stochastic Cahn-Hilliard equation (see Theorem 3.8), avoiding Gronwall's lemma. (d) Building on the results from (c), and using an L ∞ -bound for the solution of Scheme 3.1 (Lemma 5.1), along with error estimates in stronger norms (Lemma 5.2), we show uniform convergence of iterates on large subsets of Ω (Theorem 5.5). This intermediate result then implies the second main result of the paper: the convergence in probability of iterates of Scheme 3.1 to the sharp interface limit in Theorem 5.7 for sufficiently large γ. In particular, we show that the numerical solution of (1.1) uniformly converges in probability to 1, −1 in the interior and exterior of the geometric interface of the deterministic Mullins-Sekerka problem (5.1), respectively. As a consequence we obtain uniform convergence of the zero level set of the numerical solution to the geometric interface of the Mullins-Sekerka problem in probability; cf. Corollary 5.8.
The error analysis below in particular identifies proper balancing strategies of numerical parameters with the interface width that allow to approximate the limiting sharp interface model for realistic problem setups, and motivates the use of space-time adaptive meshes for numerical simulations; see e.g. [24]. In Section 6, we present computational studies which evidence asymptotic properties of the solution for different scalings of the noise term. Our studies suggest the deterministic Mullins-Sekerka problem as sharp-interface limit already for γ ≥ 1; we observe this in simulations for spatially colored, as well as for the space-time white noise. In contrast, corresponding simulations for γ = 0 indicate that the sharp-interface limit is a stochastic version of the Mullins-Sekerka problem; see Section 6.4.
To sum up, the convergence analysis presented in this paper is a combination of a perturbation and discretization error analysis. The latter depends on stability properties of the proposed numerical scheme: higher-moment energy estimates for the Scheme 3.1, a discrete spectral estimate for the related deterministic variant, and a local error analysis on the sample set Ω are crucial ingredients of our approach. The techniques developed in this paper constitute a general framework which can be used to treat different and/or more general phase-field models including the stochastic Allen-Cahn equation, and apply to settings which involve multiplicative noise, driving trace-class Hilbert-space-valued Wiener processes, and bounded polyhedral domains D ⊂ R 3 , as well.
The paper is organized as follows. Section 2 is dedicated to the analysis of the continuous problem. The time discretization Scheme 3.1 is proposed in Section 3 and rates of convergence are shown, while Section 4 extends this convergence analysis to its finite-element discretization. The convergence of the numerical discretization to the sharp-interface limit is studied in Section 5. Section 6 contains the details of the implementation of the numerical schemes for the stochastic Cahn-Hilliard and the stochastic Mullins-Sekerka problem, respectively, as well as computational experiments which complement the analytical results.
2. The stochastic Cahn-Hilliard equation 2.1. Notation. For 1 ≤ p ≤ ∞, we denote by L p , · L p the standard spaces of p-th order integrable functions on D. By (·, ·) we denote the L 2 -inner product, and let · = · L 2 . For k ∈ N we write H k , · H k for usual Sobolev spaces on D, and H −1 = (H 1 ) . We define L 2 0 := {φ ∈ L 2 ; D φ dx = 0}, and for v ∈ L 2 we denote its zero mean counterpart as v ∈ L 2 0 , i.e., v := v − 1 |D| D v dx. We frequently use the isomorphism (−∆) −1 : In particular, 0 . Throughout the paper, C denotes a generic positive constant that may depend on D, T , but is independent of ε. Definition 2.1. Let u ε 0 ∈ L 2 (Ω, F 0 , P; H 1 )∩L 4 (Ω, F 0 , P; L 4 ) and denote H 2 = {ϕ ∈ H 2 , ∂ n ϕ = 0 on ∂D}. Then, the process is called a strong solution of (1.1) if it satisfies P-a.s. and for all 0 ≤ t ≤ T The following lemma establishes existence and bounds for the strong solution u of (1.1) and for the chemical potential w from (1.2b); cf. [13, Section 2.3] for a proof of i), while ii) follows similarly as part i) by the Itô formula and the Burkholder-Davis-Gundy inequality.
There exists a unique strong solution u of (1.1), and there hold

2.4.
Error bound between u of (1.1) and u CH of (1.1) with g ≡ 0. In [3] the authors study the convergence of the solution of the stochastic Cahn-Hilliard equation (1.1) to the deterministic sharp-interface limit. In particular, they show the convergence in probability of the solution u of (1.1) to the approximation u A of u CH for sufficiently large γ > 0. Apart from the spectral estimate (2.1), a central ingredient of their analysis is the use of a stopping time argument to control the drift nonlinearity. The stopping time which, in our setting, is defined as ds > ε σ 0 for some constant σ 0 > 0, enables the derivation of the estimates in Lemma 2.2 below up to the stopping time T ε on a large sample subset that satisfies P[Ω 1 ] → 1 for ε ↓ 0, for some constant κ 0 . On specifying the condition (A) below it can be shown that T ε ≡ T , which yields Lemma 2.2. In this section we extend the work [3] by showing a strong error estimate for u − u CH in Lemma 2.3.
In Section 3 we perform an analogous analysis on the discrete level by using a stopping index J ε , and a set Ω 2 which are discrete counterparts of T ε and Ω 1 , respectively. Both 5 approaches require a lower bound for the noise strength γ to ensure, in particular, positive probability of the sets Ω 1 and Ω 2 , respectively.
For the analysis in this section we require the following assumptions to hold.
Assumption (A) ensures positivity of all exponents in the estimates in the lemmas of this section. The following lemma relies on the spectral estimate (2.1) and is a consequence of [3, Theorem 3.10] for p = 3, d = 2, where a slightly different notational setup is used.
A closer inspection of the proofs in [3] (cf. [3,Lemma 4.3] in particular) reveals that the parameter l can be chosen arbitrarily large in the above theorem.
We now use Lemma 2.2 to show bounds for the difference u − u CH in different norms.

A time discretization Scheme for (1.1)
For fixed J ∈ N, let 0 = t 0 < t 1 < · · · < t J = T be an equidistant partition of [0, T ] with step size k = T J , and ∆ j W := W (t j ) − W (t j−1 ), j = 1, . . . , J. We approximate (1.1) by the following scheme: The solvability and uniqueness of {(X j , w j )} j≥1 , as well as the P-a.s. conservation of mass of {X j } j≥1 are immediate.
For the error analysis of Scheme 3.1, we use the iterates (X j CH , w j CH ) J j=0 ⊂ H 1 2 which solve Scheme 3.1 for g ≡ 0. The following lemma collects the properties of these iterates from [16,17]. We remark that, compared to [16,17], the results are stated in a simplified (but equivalent) form, which is more suitable for the subsequent analysis.
Assume in addition u ε 0 H 3 ≤ Cε −pCH . Then for k ≤ Cε lCH , and C 0 > 0 from (2.1) it holds iv) max Proof. The proof of i), ii), iv), v) is a direct consequence of [16, Lemma 3, Corollary 1, To show iii), we use the Gagliardo-Nirenberg inequality and [16, inequality (76)], ii), iv) to get the following L ∞ -error estimate for k ≤ Cε lCH , and some l CH > 0, The numerical solution of Scheme 3.1 satisfies the discrete counterpart of the energy estimate in Lemma 2.1, i). The time-step constraint in the lemma below is a consequence of the implicit treatment of the nonlinearity; see the last term in (3.2), its estimate (3.3), and (3.4); the lower bound for admissible γ has the same origin.
Adding both equations then leads to P-a.s.
Note that the third term on the left-hand side reflects the numerical dissipativity in the scheme. We can estimate the nonlinear term as (cf. [15, Section 3.1]), where we employ the notation f(u) := |u| 2 − 1, i.e., f (X j ) = f(X j )X j . The third term on the right-hand side again reflects numerical dissipativity.
By ω ∈ Ω fixed, and ϕ = (−∆) −1 [X j − X j−1 ](ω) in Scheme 3.1, we eventually have P-a.s., which together with ∆ −1/2 g ≤ C yields the estimate Hence, using this estimate, and exploiting again the inherent numerical dissipation of the scheme we can estimate 1 2ε We substitute (3.2) along with the last inequality into (3.1) and get ε 2 which motivates time-steps k < 2ε 3 . Next, by using the second equation in Scheme 3.1, we can rewrite the first term on the right-hand side as (3.6) On recalling f (X j ) = f(X j )X j , we rewrite the remaining term as (3.7) Thanks to the embeddings L s → L r (r ≤ s), and the Cauchy-Schwarz and Young's inequalities, The leading term may now be controlled by the numerical dissipation term in (3.2). Finally, by the Poincaré's inequality, we estimate By combining the above estimates for A 3,1 , A 3,2 we obtain an estimate for (3.7).
ii) The second estimate can be shown along the lines of the first part of the proof by applying max j before taking the expectation in (3.8). The additional term that arises from the terms A 2 , A 4 in (3.5) can be rewritten by using the second equation in Scheme 3.1, where the equality in the second line follows from the zero mean property of the noise.
The last sum in (3.9) is a discrete square-integrable martingale, and by the independence properties of the summands, the Poincaré inequality and the energy estimate i) we have Therefore, (3.9) can be estimated using the discrete BDG-inequality (see Lemma 3.3) and part i) by iii) We show assertion iii) for p = 2 1 . By collecting the estimates of the terms in (3.5) in part i) (cf. (3.6), 3.7)) we deduce from (3.4) that Multiply this inequality with E(X j ) and use the identity , the estimate ε 2γ+1 ≤ ε 4 0 ε 2γ−3 , Young's inequality, and the generalized Hölder's inequality to We note that to get the above estimate we employed the reformulation E( ) on the right-hand side. By Poincaré's inequality, the last term in (3.11) may be bounded as After summing-up in (3.11) and taking expectations we get for any j ≤ J that (3.12) where the third term is bounded via (3.8) in part ii), and the statement then follows from the discrete Gronwall inequality. For p = 2 r , r = 2, we may now argue correspondingly: we start with (3.11), which we now multiply with |E(X j )| 2 . Assertion iii) now follows via induction with respect to r.
iv) The last estimate follows analogously to ii) from the BDG-inequality and iii).
The error analysis of the implicit Scheme 3.1 in the subsequent Section 3.1 involves the use of a stopping index J ε , and an associated random variable 1 {j≤Jε} that is measurable w.r.t. the σ-algebra F t j , but not w.r.t. F t j−1 . This issue prohibits the use of the standard BDG-inequality since 1 {j≤Jε} is not independent of the Wiener increment ∆ j W . The following lemma contains a discrete BDG-inequality which will be used in Section 3.1. We take {F t j } J j=0 to be a discrete filtration associated with the time mesh {t j } J j=0 ⊂ [0, T ] on (Ω, F, P). Lemma 3.3. For every j = 1, . . . , J, let F j be an F t j -measurable random variable, and Proof. We start by noting that With this identity, we obtain is also a discrete square-integrable martingale. Hence, by the L 2 -maximum martingale inequality, using the independence of 1 {j≤τ } F j and ∆ W for j < it follows that The assertion of the lemma then follows from (3.13) and (3.14). 13 3.1. Error analysis. Denote Z j := X j − X j CH , use Scheme 3.1 for a fixed ω ∈ Ω, and choose ϕ = (−∆) −1 Z j (ω), ψ = Z j (ω). We obtain P-a.s.
We use Lemma 3.1, v) to obtain a first error bound.
, and let k ≤ Cε lCH with l CH ≥ 3 from Lemma 3.1 be sufficiently small. There exists C > 0, such that P-a.s. and for all 1 ≤ ≤ J, Proof. 1. Consider the last term on the left-hand side of (3.15). On recalling Z j = X j −X j CH , by a property of f , see [17, eq. (2.6)], and Lemma 3.1, iii), we get for some C > 0 (3.17) 2. In order to later keep a portion of ∇Z j 2 on the left-hand side of (3.15) we use the identity We apply Lemma 3.1, v) to get a lower bound for the first term on the right-hand side, On noting ε < 1, we estimate the remaining nonlinearities in (3.18) using Lemma 3.1, iii), 14 3. We insert the estimates from the steps 1. and 2. into (3.15), and use the bound 4. We sum the last inequality from j = 1 up to j = , and consider max j≤ . On noting Z 0 = 0, we obtain P-a.s. (3.20) Hence, the implicit version of the discrete Gronwall lemma implies for sufficiently small k ≤ k 0 (D) that P-a.s.
which concludes the proof.
In the deterministic setting (g ≡ 0), an induction argument, along with an interpolation estimate for the L 3 -norm is used to estimate the cubic error term on the right-hand side of (3.16); cf. [16]. In the stochastic setting, this induction argument is not applicable any more, which is why we separately bound errors in (3.16) on two subsets Ω 2 and Ω \ Ω 2 . In the first step, we study accumulated errors on Ω 2 locally in time, and therefore mimic a related (time-continuous) argument in [3]. We introduce the stopping index 1 ≤ J ε ≤ J where the constant σ 0 > 0 will be specified later. The purpose of the stopping index is to identify those ω ∈ Ω where the cubic error term is small enough. In the sequel, we estimate the terms on the right-hand side of (3.16), putting = J ε . Clearly, the part .20) is bounded by ε σ 0 ; the remaining part will be denoted by For 0 < κ 0 < σ 0 , we gather those ω ∈ Ω in the subset where the error terms in Lemma 3.4 which cannot be controlled by the stopping index J ε do not exceed the larger error threshold ε κ 0 . The following lemma quantifies the possible error accumulation in time on Ω 2 up to the stopping index J ε in terms of σ 0 , κ 0 > 0, and illustrates the role of k in this matter; it further provides a lower bound for the measure of Ω 2 correspondingly. ii) The proof uses the discrete BDG-inequality (Lemma 3.3), which is suitable for the implicit Scheme 3.1; we use the higher-moment estimates from Lemma 3.2, iii) to bound the last term in R Jε .
Let Ω c 2 := Ω \ Ω 2 . We use Markov's inequality to estimate P We first estimate the last term in R Jε : interpolation of L 3 between L 2 and H 1 , then of L 2 between H −1 and H 1 (D ⊂ R 2 ) and the Young's inequality yield The leading term on the right-hand side is absorbed on the left-hand side of the inequality in Lemma 3.4, which is considered on the whole of Ω; the expectation of the last term (on the whole of Ω) is bounded via Lemma 3.2, iv) by Ck 2 ε 4 |E(u ε 0 )| 2 + 1 . For the first term in R Jε we use the discrete BDG-inequality (Lemma 3.3) to bound its expectation by In order to benefit from the definition of J ε for its estimate, we split the leading summand, Putting things together leads to 3. Consider the inequality in Lemma 3.4 on Ω 2 . The estimate ii) then follows after taking expectation, using (3.23) and recalling the definition of J ε .
Compared to assumption (A), the less restrictive lower bound for γ is due to the use of the discrete spectral estimate (see Lemma 3.1, v)), which introduces a factor ε −4 that is absorbed into ε 3 2 κ 0 in the proof below. Consequently we only need to require γ ≥ 19 3 in order to ensure positive probability of Ω 2 .
Proof. 1. Assume that J ε < J on Ω 2 ; we want to verify that Use (3.22), and the estimate Lemma 3.5 i) to conclude on Ω 2 .

2.
Recall that the last part in Lemma 3.5 yields P[ We collect the requirements on the analytical and numerical parameters: For sufficiently small ε 0 ≡ (σ 0 , κ 0 ) > 0 and l CH ≥ 3 from Lemma 3.1, and arbitrary We note that, except for the higher regularity of the initial condition, the assumption (B) is less restrictive than the assumption (A) from Section 2. Furthermore, the condition E(u ε 0 ) < C can be weakened to E(u ε 0 ) < Cε −α , α > 0, cf. [17, Assumption (GA 2 )]. Lemma 3.7. Suppose (B). Then there exists C > 0 such that Proof. Recall the notation from (3.20), and split E[ Due to assumption (B) it follows directly from Lemma 3.5, ii) and Lemma 3.6 that In order to bound E[1 Ω c 2 A J ], we use the embedding L 4 ⊂ H −1 which along with the higher-moment estimate from Lemma 3.2 iv) implies that . Next, we note that by Lemma 3.5 it follows that Hence, using the Cauchy-Schwarz inequality we get After inspecting (3.24), (3.25) we note that the statement follows by assumption (B), since the latter contribution dominates the error. The dominating error contribution in Lemma 3.7 comes from the term E[1 Ω c 2 A J ]. This is in contrast to Section 2 where the error contribution from the set Ω c 1 can be made arbitrarily small, due to the additional parameter l > 0 in Lemma 2.2 which can be chosen arbitrarily large independently of the other parameters.
We are now ready to prove the first main result of this paper.
Due to condition (A) 2 it holds that σ 0 − κ 0 < 1 3 σ 0 . Consequently the contribution ε ; it is only stated explicitly to highlight the error contribution from the difference u − u CH from Section 2. Proof. We estimate the error via splitting it into three contributions, Remark 3.9. An alternative approach to Theorem 3.8 would be to follow the arguments in [23] for a related problem, which exploit a weak monotonicity property of the drift operator in (1.1), and stability of the discretization to obtain a strong error estimate for Scheme 3.1 of the form While the error tends to zero for k ↓ 0 in (3.26), this estimate is only of limited practical relevancy in the asymptotic regime where ε is small, since only prohibitively small step sizes k exp(− 1 ε ) are required in (3.26) to guarantee small approximation errors for iterates from Scheme 3.1. Moreover, the error analysis that leads to (3.26) does not provide any insight on how to numerically resolve diffuse interfaces via proper balancing of discretization parameter k and interface width ε -which is relevant in the asymptotic regime where ε 1.

Space-time Discretization of (1.1)
We generalize the convergence results in Section 3 for Scheme 3.1 to its space-time discretization. For this purpose, we introduce some further notations: let T h be a quasi-uniform triangulation of D, and V h ⊂ H 1 be the finite element space of piecewise affine, globally continuous functions, and the Riesz projection P H 1 : In what follows, we allow meshes T h for which P L 2 is H 1 -stable; see [10]. Also, we define the inverse discrete Laplacian (−∆ h ) −1 : We are ready to present the space discretization of Scheme 3.1.
, and using the definition of (−∆ h ) −1 , as well as X j h , P L 2 g ∈ L 2 0 P-a.s., such that To obtain the first identity in (3.5) for Scheme 4.1, we use ε γ (g, w j h )∆ j W = ε γ P L 2 g, w j h ∆ j W , such that the second equation in Scheme 4.1 with ψ h = P L 2 g may be applied; as a consequence, g has to be replaced by P L 2 g in the rest of equality (3.5). This modification leads to the term ∇P L 2 g in (3.6), which is again bounded by ∇g ; the bound P L 2 g L ∞ ≤ C, which is required to bound the term A 3,1 from (3.7), follows by an approximation result; cf. [7,Chapter 7]. The above steps then yield the estimate (3.8) ii'), iii'), iv') We can follow the argumentation in the proof of Lemma 3.2 without change.
for all ≤ J, provided that additionally for any β > 0, and p CH , q CH , p CH > 0. Proof. Again, we here denote by {(X j CH;h , w j CH;h )} 1≤j≤J ⊂ [V h ] 2 the solution of Scheme 4.1 for g ≡ 0, whose stability and convergence properties are studied in [16,17]. Under the assumption (4.1), [17, Theorem 3.2, (iii)] provides the bound max 0≤j≤J X j CH;h L ∞ ≤ C .
We use this bound to adapt estimate (3.17) to the present setting and get Step 2. of the proof of Lemma 3.4 involves the discrete spectral estimate (see Lemma 3.1, iv)) for {X j CH } j to handle the leading term on the right-hand side of (3.17) -which we do not have for {X j CH;h } j in the present setting. Therefore, we perturb the leading term on the right-hand side of the last inequality, and use the L ∞ -bounds for X j CH , X j CH;h , as well as the mean-value theorem to conclude The remaining steps in the proof of Lemma 3.4 now follow with only minor adjustments.

Claim 3. Additionally assume (4.1). Then Lemma 3.5 holds for {Z
Proof. The proof for Lemma 3.5 directly transfers to the present setting. Proof. We only need to adapt the interpolation argument for L 3 to the present setting, starting with the estimate Z i By the definition of the H −1 -norm, the definition and H 1 -stability of the L 2 -projection, and again the fact that (Z i h , 1) = 0, we deduce Next, we formulate a counterpart of Lemma 3.7 for the fully discrete numerical solution; as a consequence of the Claims 1 to 4 above the corollary can be proven analogically to Lemma 3.7 with the assumption (B) complemented by the additional restriction on the discretization parameters (4.1).
We are now ready to extend Theorem 3.8 to Scheme 4.1.

22
Theorem 4.2. Let u be the strong solution of (1.1), and X j h ; 1 ≤ j ≤ J the solution of Scheme 4.1. Assume (B) and (4.1). Then there exists C > 0 such that where m CH , m CH > 0.
We note that the exponents m CH , m CH > 0 in the above estimate can be determined on closer inspection of [16, Corollary 2] on assuming (4.1). Furthermore, assumption (4.1), which is a simplified reformulation of assumption 4) in [16,Corollary 2], guarantees that lim ε↓0 Proof. We split the error into three contributions, The first term is bounded by Cε 2 3 σ 0 as in Theorem 3.8. The second term is bounded by ε m CH thanks to [16, Corollary 2] (stated here in a simplified form), provided assumption (4.1) holds. The last term is bounded by C ε κ 0 max k 2 ε 4 , ε γ+ σ 0 +1 3 , ε σ 0 , ε 2γ

Sharp-interface limit
In this section, we show the convergence of iterates {X j } J j=1 of Scheme 3.1 to the solution of a sharp interface problem. Recall that in the absence of noise, the sharp interface limit of (1.1) is given by the following deterministic Hele-Shaw/Mullins-Sekerka problem: Find v MS : [0, T ] × D → R and the interface Γ MS t ; 0 ≤ t ≤ T such that for all t ∈ (0, T ] the following conditions hold: where κ is the curvature of the evolving interface Γ MS t , and V is the velocity in the direction of its normal n Γ , as well as [ ∂vMS , and F is the double-well potential; cf. [1] for a further discussion of the model. Below, we show that iterates {X j } J j=1 of Scheme 3.1 converge to the limiting Mullins-Sekerka problem (5.1); see Theorem 5.7 for a precise specification of the convergence result. For this purpose, we need sharper stability and convergence results than those available from Section 3, which also requires to tighten the assumptions (B), and so to further restrict admissible choices of γ > 0. We note that the stronger stability estimates below are derived formally using the (analytically) strong formulation of Scheme 3.1; the derivation can be justified rigorously by the respective elliptic regularity of the Laplace and the bi-Laplace operators.
In order to establish convergence to zero (for ε ↓ 0) of the right-hand side in the inequality of the lemma, we need to impose a stronger assumptions than (B); for simplicity, we assume n CH ≥ 3 2 in Lemma 3.1: (C 1 ) Assume (B), and that (σ 0 , κ 0 , γ) also satisfies For sufficiently small ε 0 ≡ (σ 0 , κ 0 ) > 0 and l CH ≥ 3 from Lemma 3.1, and arbitrary 0 < β < 1 2 the time-step satisfies Compared to assumption (B), only larger values of σ 0 , and consequently larger values of γ are admitted, as well as smaller time-steps k.
Proof. 1. We subtract Scheme 3.1 in strong form for g ≡ 0 and g ≡ 0, respectively, fix ω ∈ Ω, and multiply the resulting error equations with Z j (ω) and −∆Z j (ω), respectively. We obtain We estimate the right-hand side above as We restate the nonlinear term in (5.3) as where in the last step we used integration by parts |Z j | 2 Z j , −∆Z j = 3 Z j ∇Z j 2 . Next, we apply integration by parts to I 1 , I 2 to estimate Hence, using Poincaré, Sobolev and Young's inequalities, Lemma 3.1, ii), and assumption (B), we deduce that 2. We insert these bounds into (5.3), sum up over all time-steps, take max j≤J and expectations, We use the discrete BDG-inequality (Lemma 3.3) and the Poincaré inequality to estimate the last term as follows, We now use Lemma 3.7 to bound the right-hand side of (5.4).
p , for 2 < p < 3 via (5.5). To meet (5.8) therefore additionally requires for some p > 2 and hence A proper scenario is k = ε α for some α > 0 to meet assumption (C 1 ). We then sharpen this choice of the time-step to k = ε α for some α ≥ α > 0 to have for an arbitrary η > 0. We now choose 2 < p, s.t. p p−2 0 is sufficiently large to meet (5.10).
(3) We may proceed analogously for the second term on the right-hand side in Lemma 5.3.
Proof. We subtract Scheme 3.1 for g ≡ 0 and g ≡ 0 for a fixed ω ∈ Ω, and multiply the first error equation with −∆Z j (ω), and the second with ∆ 2 Z j (ω). We integrate by parts in the nonlinear term and obtain 1 2 We proceed as in the proof of Lemma 5.2 and rewrite the nonlinearity on the right-hand side as We estimate We estimate 3 =1 I on Ω κ,j via Lemma 3.1, ii)-iii) and the embedding H 1 → L 4 on recalling (5.6) We multiply (5.11) by 1 Ω κ,j , sum up for 1 ≤ i ≤ j, take max 1≤j≤J and expectation, employ the identity (recall, 1 Ω κ,j−1 − 1 Ω κ,j ≥ 0) use Lemmata 5.2 and 3.7 to estimate (5.12) and obtain To estimate the stochastic term we use ∂ n g = 0 on ∂D and proceed as follows, The first term on the right-hand side may be bounded by Lemma 5.2, the third term is absorbed in the left-hand side of (5.13), and for the second term we use the discrete BDGinequality (Lemma 3.3) and Lemma 3.7 to estimate Hence, the statement of the lemma follows from (5.13) and the above estimates on noting The L ∞ -estimate in the next theorem is a crucial ingredient to show convergence to the sharp-interface limit.

30
Remark 5.6. We discuss a strategy to identify admissible quadruples (σ 0 , κ 0 , γ, k) which meet assumption (C 3 ): for this purpose, we limit ourselves to a discussion of the leading term inside the maximum which defines F 2 (see Lemma 5.3), and recall Remark 5.4.
We are now ready to formulate the second main result of this paper, which is convergence in probability of the solution {X j } J j=0 of Scheme 3.1 to the solution of the deterministic Hele-Shaw/Mullins-Sekerka problem (5.1) for ε ↓ 0, provided that assumption (C 3 ) is valid, and (5.1) has a classical solution; cf. Theorem 5.7 below. The proof rests on a) the uniform bounds for {1 Ω κ,j Z j p L ∞ } J j=1 (see Theorem 5.5), and the property that lim ε↓0 max 1≤j≤J P[Ω κ,j ] = 1 (in Lemma 5.1) for the sequence {Ω κ,j } J j=1 ⊂ Ω, and b) a convergence result for {X j CH } J j=0 towards a smooth solution of the Hele-Shaw/Mullins-Sekerka problem in [17,Section 4].
For each ε ∈ (0, ε 0 ) we consider below the piecewise affine interpolant in time of the iterates {X j } J j=0 of Scheme 3.1 via Let Γ 00 ⊂ D in (5.1e) be a smooth closed curve, and (v MS , Γ MS ) be a smooth solution of (5.1) starting from Γ 00 , where Γ MS := 0≤t≤T {t} × Γ MS t . Let d(t, x) denote the signed distance function to Γ MS By [1, Theorem 5.1], assumption (D) establishes the existence of a family of smooth solutions {u ε 0 } 0≤ε≤1 which are uniformly bounded in ε and (t, x), such that if u ε CH is the corresponding solution of (1.1) with g ≡ 0, then The following theorem establishes uniform convergence of iterates {X j } J j=0 from Scheme 3.1 in probability on the sets I MS , O MS .
Remark 5.9. The numerical experiments in Section 6 suggest that the conditions on γ and k which are required for Theorem 5.7 to hold are too pessimistic; in particular, they indicate convergence to the deterministic Mullins-Sekerka/Hele-Shaw problem already for γ = 1, k = O(ε).

Computational experiments
The computational experiments are meant to support and complement the theoretical results in the earlier sections: • Convergence to the deterministic sharp-interface limit (5.1) for the space-time white noise in Section 6.3. We study pathwise convergence of the white noise-driven simulations to the deterministic sharp interface limit, which is a scenario beyond the one for regular trace-class noise where Theorem 5.7 and Corollary 5.8 establish convergence in probability. • Pathwise convergence to the stochastic sharp interface limit (6.4) (introduced in Section 6.2 below) for spatially smooth noise in Section 6.4, where we also examine the sensitivity of numerical simulations with respect to the mesh refinement.
6.1. Implementation and adaptive mesh refinement. For the computations below we employ a mass-lumped variant of Scheme 4.1 where the standard L 2 -inner product in Scheme 4.1 is replaced by the discrete (mass-lumped) In all experiments we take D = (0, 1) 2 ⊂ R 2 and g is taken to be a constant. We note that an implicit Euler finite element scheme similar to Scheme 6.1 has been used previously in [19], which also performs simulations to study long time behavior of the system for different strengths of the (space-time white) noise with fixed ε. For a given initial interface Γ 00 we construct an ε-dependent family of initial conditions where d 0 is the signed distance function to Γ 00 . Consequently, {u ε 0 } ε>0 have bounded energy and contain a diffuse layer of thickness proportional to ε along Γ 00 , and u ε 0 (x) ≈ −1, u ε 0 (x) ≈ 1 in the interior, exterior of Γ 00 , respectively. The construction ensures that D u ε 0 dx → m 0 for ε → 0, where m 0 is the difference between the respective areas of the exterior and interior of Γ 00 in D. For convenience we set u ε,h variables which approximate the increments of a Q-Wiener process on a probability space (Ω, F, P) which is given by where {e i } i∈N is an orthonormal basis in L 2 (D), {β i } i∈N are independent real-valued Brownian motion, and {λ i } i∈N are real-valued coefficients such that Qe i = λ 2 i e i , i ∈ N. In order to preserve mass the noise is required to satisfy P-a.s. D W (t, x) dx = 0, t ∈ [0, T ].
In the experiments below we consider two types of Wiener processes: a smooth (finite dimensional) noise and a L 2 0 -cylindrical Wiener process (space-time white noise). The smooth noise is given by where ∆ j β k = β k (t j ) − β k (t j−1 ) are independent scalar-valued Brownian increments. The discrete approximation of the smooth noise is then constructed as where φ (x m ) = δ m , = 1, . . . , L are the (standard) nodal basis function of V h , i.e., V h = span{φ , = 1, . . . , L}. The space-time white noise (Q = I) is approximated as (cf. [5]) In order to preserve the zero mean value property of the noise we normalize the increments as The Wiener process is simulated using a standard Monte-Carlo technique, i.e., for ω m ∈ Ω, m = 1, . . . , M , we approximate the Brownian increments in (6.2),(6.3) as ∆ j β (ω m ) ≈ √ kN j (0, 1)(ω m ), where N j (0, 1)(ω m ) is a realization of the Gaussian random number generator at time level t j . The discrete nonlinear systems related to (realizations of) the scheme (6.1) are solved using the Newton method with a multigrid linear solver.
To increase the efficiency of the computations we employ a pathwise mesh refinement algorithm. For a realization X j h,m := X j h (ω m ), ω m ∈ Ω of the V h -valued random variable X j h we define η grad (x) = max{|∇X j h,m (x)|, |∇X j−1 h,m (x)|} and refine the finite element mesh in such a way that h(x) = h min if εη grad (x) ≥ 10 −2 and h(x) ≈ h max if εη grad (x) ≤ 10 −3 ; the mesh produced at time level j is then used for the computation of X j+1 h,m . The adaptive algorithm produces meshes with mesh size h = h min along the interfacial area and h ≈ h max in the bulk where u ≈ ±1, see Figure 3 for a typical adapted mesh. In our computations we choose h max = 2 −6 and h min = π 4 ε, i.e. h min = h max for ε ≥ 1/(16π) and h min scales linearly for smaller values of ε.
In the presented simulations, mesh refinement did not appear to significantly influence the asymptotic behavior of the numerical solution. This is supported by comparison with additional numerical simulation on uniform meshes. The observed robustness of numerical simulations with respect to the mesh refinement can be explained by the fact that the asymptotics are determined by pathwise properties of the solution on a large probability 35 set. This conjecture is supported by the convergence in probability in Theorem 5.7 and Corollary 5.8. In the present setup the (possible) bias due to the pathwise adaptive-mesh refinement did not have significant impact on the results. In general, the use of adaptive algorithms with rigorous control of weak errors may be a preferable approach, cf. [24]. 6.2. Stochastic Mullins-Sekerka problem and its discretization. We consider the following stochastic modification of the Mullins-Sekerka problem (5.1) We note that the only difference between (5.1) and (6.4) is in the equations (5.1a), (6.4a), respectively. Alternatively equation (6.4a) can be stated in an integral form as For the approximation of the stochastic Mullins-Sekerka problem (6.4), we adapt the unfitted finite element approximation for the deterministic problem (5.1) from [6]. In particular, let Γ j−1 be a polygonal approximation of the interface Γ at time t j−1 , parameterized by Y j−1 h ∈ [V h (I)] 2 , where I = R/Z is the periodic unit interval, and where V h (I) is the space of continuous piecewise linear finite elements on I with uniform mesh size h. Let π h : C(I) → V h (I) be the standard nodal interpolation operator, and let ·, · denote the L 2 -inner product on I, with ·, · h the corresponding mass-lumped inner product. Then we In the above, ρ denotes the parameterization variable, so that |[Y j−1 ] ρ | is the length element on Γ j−1 and ν j−1 h ∈ [V h (I)] 2 is a nodal discrete normal vector, see [6] for the precise definitions.
6.3. Convergence to the deterministic sharp-interface limit.
6.3.1. One circle. We set γ = 1, g = 8π and consider the discrete space-time white noise (6.3). We note that the considered space-time white noise does not satisfy the smoothness assumptions required for the theoretical part of the paper (i.e., γ > 1 and tr(∆Q) < ∞), however the numerical results indicate that for ε ↓ 0 the computed evolutions still converge to the deterministic Mullins-Sekerka problem (5.1).
The numerical studies below are performed using the scheme (6.1) with adaptive mesh refinement. The time-step size for ε = 2 −i /(64π), i = 0, . . . , 4 was k i = 2 −i 10 −5 . The motivation of the different choice of the time-step is to eliminate possible effects of numerical damping and to ensure the convergence of the Newton solver for smaller values of ε.
For each ε we use the initial condition u ε,h 0 that approximates a circle with radius R = 0.2. Since circles are stationary solutions of the deterministic Mullins-Sekerka problem, the convergence of the numerical solution for the stochastic Cahn-Hilliard equation to the solution of the Mullins-Sekerka problem for ε ↓ 0 can be determined by measuring the deviations of the zero level-set of the solution X j h , j = 1, . . . , J from the circle with radius R = 0.2 for a sufficiently large computational time. We note that the zero level-set of the initial condition u ε,h 0 above, exactly approximates the corresponding stationary solution of the Mullins-Sekerka problem, but it is not a stationary solution of the corresponding (discrete) deterministic Cahn-Hilliard equation, i.e., of (6.1) with g ≡ 0. In order to obtain the optimal phasefield profile across the interfacial region, we let u ε,h 0 relax towards the discrete stationary state by computing with (6.1) for g ≡ 0 for a short time and then use that discrete solution as the actual initial condition for the subsequent simulations.
The results in Figure 1 indicate that for decreasing ε the evolution of the zero level set of the numerical solution approaches the solution of the deterministic Mullins-Sekerka model, which is represented by the stationary circle with radius 0.2. We observe that the deviations of the interface from the circle are decreasing for smaller ε.
6.3.2. Two circles. In this experiment we consider the same setup as in the previous one with an initial condition which consists of two circles with radii R 1 = 0.15 and R 2 = 0.1, respectively. The evolution of the solution is more complex than in the previous experiment as the interface undergoes a topological change. To minimize the Ginzburg-Landau energy, the left (larger) circle grows, the right (smaller) circle shrinks and the resulting steady state is a single circle with mass equal to the mass of the two initial circles; see Figure 2 for an example of a deterministic evolution with ε = 1/(512π). In Figure 3   6.4. Comparison with the stochastic Mullins-Sekerka model. We use the numerical scheme (6.1) to study the case of non-vanishing noise, i.e., γ = 0, with the discrete approximation of the smooth noise (6.2). The noise is symmetric across the center of the domain in order to facilitate an easier comparison with the Mullins-Sekerka problem. The computations below are pathwise, i.e., in the graphs below we display results computed for a single realization of the Wiener process. If not mentioned otherwise we use the time-step size k = 10 −5 .
The initial condition is taken to be the ε-dependent approximation of a circle with radius R = 0.2 as in §6. 3  to a stationary state and then use the stabilized profile X 0 h := X js h as an initial condition for the computation. The zero level-set of the stationary solution X js h is a circle with perturbed radius R = 0.2 + O(ε), where in general the perturbation O(ε) also depends on the finite element mesh. To compensate for the effect of the perturbation in the initial condition for larger values of ε we represent the interface by a level set Γ j u Γ := {x ∈ D; X j h (x) = u Γ } (i.e., Γ j 0 is the zero level set of the discrete solution at time level t j ) where the values u Γ = X js (0.2, 0), i.e., it is the "compensated" level-set for which the stationary profile Γ js u Γ coincides with the circle with radius R = 0.2. The usual value for the "compensated" level-set was u Γ ≈ 0.27 in the computations below.
We observe that in order to properly resolve the spatial variations of the noise it is necessary to use a mesh size smaller or equal to h max = 2 −7 for the discretization of the Cahn-Hilliard equation. The computations for the Mullins-Sekerka problem, using the scheme (6.5), were more sensitive to the mesh size, and an accurate resolution for the considered noise required a mesh size h max = 2 −8 , cf. Figure 4 which includes the results for h max = 2 −8 as well as h max = 2 −7 .
In Figure 4 we compare the evolution for the stochastic Cahn-Hilliard equation for ε = 1/(32π), ε = 1/(64π) on a uniform mesh with h = 2 −7 , h = 2 −8 , respectively, with the evolution of the stochastic Mullins-Sekerka problem (6.4) on uniform meshes with h = 2 −7 , h = 2 −8 , respectively, for a single realization of the noise. We also include results for ε = 1/(128π), ε = 1/(512π), where to make the computations feasible we employ the adaptive algorithm with h max = 2 −8 and h max = 2 −9 , h max = 2 −11 , respectively. Furthermore, in order to ensure convergence of the Newton solver for ε = 1/(512π) we decrease the timestep size k = 10 −6 . To be able to directly compare with the results for ε = 1/(512π), we take the values of the realization of the noise generated with step size k = 10 −5 , which was used in the other simulations, and to obtain values at the intermediate time levels we employ linear interpolation in time. We observe that the results in Figure 4   Mullins-Sekerka model are more sensitive to the mesh size, i.e., the graph for the mesh with h = 2 −7 differs significantly from the remaining results. For the mesh with h min = 2 −8 the results for the stochastic Mullins-Sekerka model are in good agreement with the results for the stochastic Cahn-Hilliard model. We note that for values smaller than ε = 1/(128π) we do not observe significant improvements of the approximation of the stochastic Mullins-Sekerka problem. This is likely caused by the discretization errors in the numerical approximation of the stochastic Mullins-Sekerka model which, for small values of ε, are greater than the approximation error w.r.t. ε in the stochastic Cahn-Hilliard equation.
From the above numerical results we conjecture that for ε ↓ 0 the solution of the stochastic Cahn-Hilliard equation with a non-vanishing noise term (γ = 0) converges to the solution of a stochastic Mullins-Sekerka problem (6.4). Formally, the stochastic Mullins-Sekerka problem (6.4) can be obtained as a sharp-interface limit of a generalized Cahn-Hilliard equation where the noise is treated as a deterministic function G 1 (t) = gẆ (t), cf. (2.3) in [3] and (1.12) in [4].
To examine the robustness of previous results with respect to adaptive mesh refinement we recompute the previous problems with the noise (6.2) using the adaptive mesh refinement algorithm with h max = 2 −6 and h min = π 4 ε. The stochastic Mullins-Sekerka model is computed with h max = 2 −6 and the mesh is refined along the interface Γ with mesh size h min = 2 −8 .
We note that with adaptive mesh refinement the results differ from those computed using uniform meshes, since the noise (6.2) is mesh dependent. For instance, in the regions with coarse mesh the noise (6.2) is not properly resolved. The computed results with the adaptive mesh refinement can be interpreted as replacing the additive noise (6.2) with a multiplicative type noise that has lower intensity when u ≈ ±1. The presented computations contain an additional "geometric" factor in the numerical error that is due to the fact that the mesh is adapted according to the position of the interface, as well as due to the fact that the adaptive mesh refinement algorithm for the Mullins-Sekerka problem is different. Nevertheless, the results are still in good agreement with the stochastic Mullins-Sekerka problem, see Figure 5. In particular we observe that the convergence for smaller values of ε is more obvious for the zero level-set of the solution than in the case of uniform meshes. In Figure 5 we also include a graph ('ftilde' in pink) which was computed using a modification of scheme (6.1) with f (X j h ), ψ h replaced by f (X j h , X j−1 h ), ψ h wheref (X j h , X j−1 h ) = 1 2 (|X j h | 2 − 1)(X j h + X j−1 h ); for equal time-step size the modified scheme provides worse approximation of the Mullins-Sekerka problem due to numerical damping.