# Sharp-interface problem of the Ohta-Kawasaki model for symmetric diblock copolymers

The Ohta-Kawasaki model for diblock-copolymers is well known to the scientific community of diffuse-interface methods. To accurately capture the long-time evolution of the moving interfaces, we present a derivation of the corresponding sharp-interface limit using matched asymptotic expansions, and show that the limiting process leads to a Hele-Shaw type moving interface problem. The numerical treatment of the sharp-interface limit is more complicated due to the stiffness of the equations. To address this problem, we present a boundary integral formulation corresponding to a sharp interface limit of the Ohta-Kawasaki model. Starting with the governing equations defined on separate phase domains, we develop boundary integral equations valid for multi-connected domains in a 2D plane. For numerical simplicity we assume our problem is driven by a uniform Dirichlet condition on a circular far-field boundary. The integral formulation of the problem involves both double- and single-layer potentials due to the modified boundary condition. In particular, our formulation allows one to compute the nonlinear dynamics of a non-equilibrium system and pattern formation of an equilibrating system. Numerical tests on an evolving slightly perturbed circular interface (separating the two phases) are in excellent agreement with the linear analysis, demonstrating that the method is stable, efficient and spectrally accurate in space.

## Authors

• 1 publication
• 2 publications
• 5 publications
• 7 publications
• 2 publications
• 1 publication
06/19/2021

### Sharp-interface limits of the diffuse interface model for two-phase inductionless magnetohydrodynamic fluids

In this paper, we propose and analyze a diffuse interface model for indu...
01/16/2022

### Numerical study on viscous fingering using electric fields in a Hele-Shaw cell

We investigate the nonlinear dynamics of a moving interface in a Hele-Sh...
01/13/2021

### How do degenerate mobilities determine singularity formation in Cahn-Hilliard equations?

Cahn-Hilliard models are central for describing the evolution of interfa...
07/23/2021

### Circular arc polygons, numerical conformal mappings, and moduli of quadrilaterals

We study numerical conformal mappings of planar Jordan domains with boun...
07/24/2021

### Nonlinear simulation of vascular tumor growth with chemotaxis and the control of necrosis

In this paper, we develop a sharp interface tumor growth model to study ...
04/29/2020

### On a nonlocal Cahn-Hilliard model permitting sharp interfaces

A nonlocal Cahn-Hilliard model with a nonsmooth potential of double-well...
12/04/2021

### Analysis and application of a lower envelope method for sharp-interface multiphase problems

We introduce and analyze a lower envelope method (LEM) for the tracking ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The Ohta-Kawasaki (OK) model ohta1986equilibrium was originally derived by Takao Ohta and Kyozi Kawasaki to investigate mesoscopic phase separation in block copolymers. The phase separation in copolymeric substances results in the formation of two distinct regions, each rich in a particular ingredient. Domains of various shape may emerge in the system under various ratios of molecular weight of the two species. It is necessary to investigate such systems as the resulting properties are different from those observed in multiphase systems of single monomer types. The model has garnered strong interest since its emergence and has been connected to areas beyond which it was originally proposed. Examples of applications include problems in condensed matter physics and biological systems abels2018model .

In the original work of Ohta and Kawasaki ohta1986equilibrium , an energy functional was proposed to investigate the phenomenon of phase separation where both attractive (short-range) and repulsive (long-range) forces play their part in determining the configurations. The evolution equation corresponding to the functional and its steady version was first mentioned in nishiura1995 , where a connection was made between Hele-Shaw (HS) flow equations and the time-dependent OK problem. In this paper, we present a formal derivation of the corresponding sharp-interface limit using matched asymptotic expansions, and show that the limiting process leads to an HS-type moving interface problem. This allows us to recast the long-time evolution of the OK problem as a modified HS problem and focus our attention to the latter to obtain insight into the original pattern formation problem. The analytical solutions are ruled out owing to the complicated geometry and we investigate the problem mainly using numerical approaches.

The boundary integral method is a preferred choice as a numerical method for HS-type problems because it entails dimension reduction, i.e., the problem defined on a domain becomes a problem defined on the domain boundary. However, the equations of dynamics constitutes stiff equations due to the surface tension acting at the fluid-fluid interface, and without the special numerical techniques described in HLS , it is practically impossible to perform long-time numerical simulations. Several references have used this technique with great success and we refer the interested reader to zhao2020pattern ; zhao2018computation ; ShuwangPRL ; ShuwangPhysicaD ; hou2001boundary . We also note that our equations differ from the traditional HS equations ShuwangJCP in a few subtle ways. In the original HS model, the far-field boundary condition is of Neumann type which very naturally corresponds to injection/removal of the fluid. Our problem, on the other hand, is driven by a Dirichlet type boundary condition in the far-field. This renders the constraint on the integral of velocity to be different in our case. We also note that the far-field boundary is at a finite distance from the origin in our case while in the classical HS problems, the radius of the far-field boundary is infinite.

The main contribution of this paper can be summarized as follows: starting with a rescaled formulation of the OK equation, we present a matched asymptotic analysis in the long-time limit that governs the dynamics of the emerging interfaces and this leads to modified HS equations of the OK model. We then prescribe a transformation that converts the HS equations from the Poisson equation to the Laplace equation and transform the interfacial and far-field boundary conditions accordingly. The equations are then investigated using a linear analysis. We prescribe a boundary integral formulation for the Laplace equation using free-space Green’s function and we investigate the boundary integral equations numerically as the analytical solutions are known in very limited cases. The numerical methods allow us to investigate the steady-state configuration for various patterns hitherto not explored in detail. Throughout our computation, we demonstrate high accuracy which is a trademark of boundary integral computations. Nonlinear computations indicate that the interface morphologies depend strongly on the mass flux into the system before the system reaching equilibrium. Simulations of multiple equilibrating interfaces show complicated interactions between phase domains including interface alignment and coarsening.

This paper is organized as follows: In Section 2, we give a formulation for the boundary value problem of the OK equation in a rescaled form that is suitable for the asymptotic analysis using matched asymptotic expansions, which is carried out in Section 3. In Section 4, the analytical solutions of the problem are discussed. Numerical methods on the boundary integral equations, the spatial discretization of the integral equations using spectrally accurate quadrature rules, the dynamical equations, and the small-scale decomposition are discussed in Section 5. The interface is updated based on these methods. Finally, we present results of numerical simulations in Section 6 and summarize our findings in Section 7.

## 2 Formulation of the Ohta-Kawasaki phase-field model

In the framework of density functional theory, the OK problem in its dimensionless form it is ohta1986equilibrium ; nishiura1995

 FOK[ϕ]=∫Ω12(∇ϕ)2+F(ϕ)−F(ϕ−)+α2ψ(ϕ−¯ϕ)dxdy. (1)

In a domain , is the density difference, , at position and at time , where the overbar denotes the average of a quantity, e.g.

 ¯ϕ≡1|Ω|∫Ωϕdxdy. (2)

is given by the solution of the Poisson problem,

 −Δψ =ϕ−¯ϕ on Ω, (3a) ∂ψ∂n∂Ω =0 on ∂Ω, (3b) ¯ψ =0, (3c)

where the last condition is introduced to enforce the uniqueness of . Here, we use for the double-well free energy the form

 F(ϕ)=14ϕ4−12ϕ2 (4)

which has two minima at . The chemical potential is obtained by the first variation of the functional

 μ =−Δϕ+(ϕ3−ϕ)−αψ, (5a) which yields the flux j =−∇μ. (5b) The system is closed via mass conservation ∂ϕ∂t =−∇⋅j (5c) together with boundary and initial conditions j⋅n∂Ω =0,∂ϕ∂n∂Ω=0on ∂Ω, (5d) ϕ(x,0) =ϕinit(x). (5e)

Derivations of the Ohta-Kawasaki phase-field model using the gradient flow approach can be found in, e.g., zhang2006periodic ; parsons2012numerical ; le2010convergence .

## 3 The sharp-interface limit

For diblock copolymers, the long-time interface formation during phase separation that sets the small-scale related to the interface width is directly connected to the parameter , via ohnishi1999analytical ; choksi_2001 . It is thus convenient to rescale the Ohta-Kawasaki model to this regime via , , , , and . After dropping the tildes, the rescaled free energy can be written as

 FOK[ϕ]=∫Ω12ε(∇ϕ)2+ε−1(F(ϕ)−F(ϕ−))+12ψ(ϕ−¯ϕ), (6)

and thus the corresponding phase-field model

 ∂ϕ∂τ =Δμ, (7a) μ =−εΔϕ+ε−1(ϕ3−ϕ)−ψ, (7b) −Δψ =ϕ−¯ϕ, (7c) ∂ϕ∂n∂Ω =0,∂ψ∂n∂Ω=0,∂μ∂n∂Ω=0on ∂Ω, (7d) ϕ(x,0) =ϕ0(x). (7e)

Due to the small parameter multiplying the Laplace operator in the chemical potential, the problem is singularly perturbed as . While such problems have been considered before with different methods henry_2001 ; le_2008 ; nishiura1995 , we investigate this “outer” problem through matched asymptotic expansions, where asymptotic approximations for the outer problem are matched to approximations of a corresponding “inner” problem in the neighborhood of the sharp interface. Our investigation follows a similar method applied by pego_1989 for the Cahn-Hilliard equations. We assume , , and have the asymptotic expansions, , , and . Substitution into (7) yields the asymptotic problems for up to order ,

 O(ε0):∂τϕ0=Δμ0,O(ε1):∂τϕ1=Δμ1,O(ε2):∂τϕ2=Δμ2. (8)

Similarly for ,

 O(ε−1):0 =F′(ϕ0), (9a) O(ε0):μ0 =F′′(ϕ0)ϕ1+ψ0, (9b) O(ε1):μ1 =F′′(ϕ0)ϕ2+12F′′′(ϕ0)ϕ21−Δϕ0+ψ1, (9c)

and ,

 O(ε0):−Δψ0=ϕ0−¯ϕ,O(ε1):−Δψ1=ϕ1,O(ε2):−Δψ2=ϕ2. (10)

On the fixed boundary , the rescaled boundary conditions are

 ∂ϕi∂n∂Ω=0,∂μi∂n∂Ω=0,∂ψi∂n∂Ω=0,on∂Ωfori=0,1,2,…

To derive the inner problems, it is convenient to introduce a parametrization of the free interface via the arc length , and

, the normal inward-pointing vector along the free boundary, so that any point in the thin

-region around can be expressed by

 x(τ,s,z)=r(τ,s)+εzν(τ,s),

where is the distance along the inward normal direction from the sharp interface , given by

 ν(τ,s)=(−∂sr2,∂sr1),t(τ,s)=(∂sr1,∂sr2).

The relation between the derivatives of a quantity defined in inner coordinates and the derivatives in outer coordinates can be expressed as a product of matrices, see A and dreyer_wagner_2005 .

Similar to the outer problem, we assume that the inner asymptotic expansions for , , and are given by , , and . After application of the coordinate transformations to the governing equations, we obtain asymptotic subproblems for , and for the inner region. These problems are solved and matched to the outer solutions. The details of the arguments, the matching conditions for the asymptotic analysis, are carried out in A, resulting in the sharp-interface problem

 ϕ0 =±1, (11a) −Δψ0 =ϕ0−¯ϕ in Ω, (11b) Δμ0 =0 in Ω±, (11c) μ0 =σκ−ψ0 on Γ, (11d) V =12[∂μ0∂n] on Γ, (11e) ∂μ0∂n∞ =0,∂ψ0∂n∞=0 on ∂Ω, (11f)

where is the surface tension and a domain, with and the the regions where and , respectively, and is the interface between them. The normal to the latter pointing from to is called . We will, more specifically, denote by the exterior and the interior domain. The boundary of is denoted by and the jump of across the interface is given by

 [∂μ0∂n]=∂μ+0∂n−∂μ−0∂n.

Finally, the value of can be expressed as

 σ=1ϕ+−ϕ−∫ϕ+ϕ−√2(F(ϕ)−F(ϕ−)) dϕ. (12)

For the derivation of the boundary integral formulation, it is convenient to reformulate the sharp-interface problem in terms of the variable

 u:=ψ0+μ0. (13)

We consider a bounded domain where , the outer domain, and , the inner domain, are open sets of and is the moving interface separating the exterior domain and the interior domain . The interior domain is a disjoint union of finitely many open, connected components and thus The outer boundary of is denoted by . A schematic diagram of the problem is given in Fig. (1). The sharp-interface model is the following problem:

 −Δu =1−2χΩ− in Ω∖Γ, (14a) u =σκ on Γ, (14b) ∂u∂n∞ =0 on Γ∞, (14c) V =12[∂u∂n] on Γ, (14d)

where is an unknown function,

is the characteristic function of the set

, is the curvature of boundary , is the surface tension parameter, the operator is the normal derivative where denotes the normal directed from to . While the function is continuous, the derivative of suffers a jump across the interface and is given by , where and are the solutions of the OK problem in the exterior and interior domains respectively. The interface moves due to the velocity .

To eliminate the source term in the field equation and recast the problem in terms of the Laplace equation, we introduce a new function defined as

 w=u+(1−2χΩ−)4|x|2, (15)

where . Then the functions and are replaced by and in and respectively. The boundary condition Eq. (14b) on splits into conditions on and as follows:

 w− =σκ−|x|24, (16) w+ =σκ+|x|24. (17)

We also transform the far-field boundary condition Eq. (14c) to

 ∂w+∂n∞=12x∞⋅n∞, (18)

where is a point on the outer boundary and is the outward normal at . The normal velocity of the interface separating the interior and the exterior domain becomes

 V=12[∂w∂n]−12x⋅n, (19)

where, as in Eq. (14d), .

## 4 Analytical solution of original equations

It is not possible to find analytical solutions of the OK equations for arbitrary geometry and multiply connected regions. However, for simplified cases, like when is a circular domain centered at origin and a circular domain of smaller radius and centered at zero, it is possible to find an analytical solution. In such a case Barua , the solution inside is obtained as

 u−=14(x2+y2−R2)+σR. (20)

Similarly, in the exterior domain, the solution of the boundary value problem of the Poisson equation in coordinates is given by

 u+(r)=−r24+(R2∞2)logr+σR+R24−R2∞2logR. (21)

In steady state, the interface between the two domains does not move () and Eq. (14d) requires the normal derivative of to be continuous. From Eq. (20) and (21), we get

 ∂u∂n∣∣∣R− =R2, (22a) ∂u∂n∣∣∣R+ =−R2+R2∞2R. (22b)

Equating the two gives an additional relation between the radii of the interior and the total domain,

 R∞=√2R, (23)

which simply states that the area of the interior and exterior domains are equal, as expected for a symmetric diblock copolymer configuration in steady state.

The solution of the OK equations can be extended further via linear analysis on a domain with the shape of a slightly perturbed circle of the form

 r(t,R,θ)=R(t)+δ(t)coskθ,0≤θ<2π, (24)

where is the radius of the circle and is a small perturbation with . Thus, by continuity of the problem, we expect , at least for , where is possibly a short period of time. In this case, it is easier to work with the transformed equations and we presume that the solution in polar coordinates is given by

 (25)

where is the zeroth order solution and is the first order solution. A straightforward computation yields the zeroth order solution as

 w−0 =σR−R24, (26a) w+0 =R2∞2logr+σR+R24−R2∞2logR. (26b)

Next we compute the first order corrections and in this case, is of the form where

 A−=σ(k2−1)Rk+2−12Rk−1. (27)

The function is of the form where

 A+ = RkR2k+R2k∞[σ(k2−1)R2+R2−R2∞2R], (28) B+ = RkR2k∞R2k+R2k∞[σ(k2−1)R2+R2−R2∞2R]. (29)

Once the functions and are available up to first order, we may proceed to calculate the velocity of the interface as

 V≈˙r=˙R+˙δcoskθ (30)

where the “dot” on the respective variables indicate derivative with respect to time. The expression on the right of Eq. (30) captures the interface velocity up to first order. We equate the right hand side of Eq. (30) to the right hand side of Eq. (19) and obtain

 ˙R =R2∞/4R−R/2, (31) ˙δ =[−R2∞/R2+k(t2−t3)/2−kt1/2−1/2]δ. (32)

where

 t1 =σ(k2−1)/R3−1/2, (33) t2 =p1R2k−1/(R2k+R2k∞), (34) t3 =p1R2k∞/(R(R2k+R2k∞)), (35) p1 =σ(k2−1)/R2+R/2−R2∞/(2R). (36)

These solutions are used later on to validate our numerical methods.

## 5 Numerical methods

In this section, we describe the numerical methods including the derivation of the boundary integral equation, its solution, and methods to update the interface. The switch from differential equation to boundary integrals results in a dimension reduction as the original PDE problem should be solved over a domain while the integral equations only have to be solved on the boundary.

### Mathematical preliminaries

We observe that the interface , on which we have to solve the integral equation, is a union of disjoint, smooth, and closed curves where is the boundary of the region We assume that each interface is represented by

 ∂Ω−k={x(α,t)=(x(α,t),y(α,t)):0≤α<2π}, (37)

where the function is analytic and -periodic in the parameter . The local tangent and the normal vectors to the interface are

 s=(xα,yα)/sαandn=(yα,−xα)/sα (38)

respectively, where and are the derivatives w.r.t. to and is the local variation of arc length. If we introduce the angle tangent to the interface, then we may write and the curvature

### Boundary integral formulation

The introduction of the function in Eq. (15) allows us to transform the Poisson equation in the original problem to the Laplace equation. We further wish to recast the latter using boundary integral formulation. Consider the free space Green’s function . We then write the solution to the interior problem as a combination of single layer and double layer potential, i.e.,

 w−(x)=∫Γ{∂w−(x′)∂n(x′)G(x,x′)−w−(x′)∂G∂n(x′)}ds′, (39)

for As , we have

 12(σκ−|x|24)=∫Γ{∂w−(x′)∂n(x′)G(x,x′)−w−(x′)∂G∂n(x′)}ds′. (40)

Similarly for the exterior problem,

 w+(x)=~w∞−∫Γ{∂w+(x′)∂n(x′)G(x,x′)−w+(x′)∂G∂n(x′)}ds′, (41)

for , where is an unknown to be solved. As , we have

 12(σκ+|x|24)=~w∞−∫Γ{∂w+(x′)∂n(x′)G(x,x′)−w+(x′)∂G∂n(x′)}ds′. (42)

Adding equations (40) and (42) together, we have

 σκ=~w∞−∫Γ2VG(x,x′)ds′−∫Γ(x′⋅n′)G(x,x′)ds′+∫Γ|x′|22∂G∂n(x′)ds′. (43)

Eq. (43) is the boundary integral equation that we solve numerically. An additional equation is needed to complete the problem. To this end, we integrate in and in , and we then use the divergence theorem to get and . Subtracting these two equations and using equation (19), we get

 J=∫ΓVds=12Atotal−A−, (44)

where is the total area enclosed by and is the area enclosed by . We solve for and the normal velocity using equations (43) and (44). The physical meaning of in the integral equation is evident: It is the value of at corresponding to the flux given in the right hand side of Eq. (44). Our formulation thus allows us to investigate the (unknown) Dirichlet condition at the far-field corresponding to a (known) Neumann condition.

### Solving the integral equations

The boundary integral equation (43) in equal arc length parameter is given by

 ~w∞ −∫Γ2V(x(α′))G(x(α),x(α′))sα(α′)dα′ =σκ+∫Γ(x(α′)⋅n(α′))G(x(α),x(α′))sα(α′)dα′ −∫Γ|x(α′))|22∂G(x(α),x(α′))∂n(x(α′))sα(α′)dα′. (45)

This along with Eq. (44) should be solved to find the velocity of the interface as well as . We use the Nyström method to discretize the integral equations using highly accurate quadrature rules on the various integrals in Eq. (45). We discretize each of the curves using marker points using equal arc length parametrization where . We choose for some positive integer . Next, we investigate the smoothness of the various integrals in Eq. (45).

#### Double-layer potential

The kernel of the integral does not a have a singularity as with Thus, an application of trapezoidal or alternating point quadrature is enough to ensure spectral accuracy SI . One may also apply the hybrid Gauss-trapezoid quadrature rules derived using the Euler-Maclaurin formula, as suggested in alpert1999hybrid .

#### Single-layer potential

The second integrals, both in the left and right hand side of Eq. (45), possess a logarithmic singularity and cannot be handled by trapezoidal rule as it is only second-order accurate. However, the integration can be performed by first splitting the log kernel as

 (46)

and then by applying the additive rule of integration. The kernel of the integration is still singular at , but the use of a Hilbert transform HLS or quadrature referred in kress1995numerical results in spectral accuracy. In this work we use the method suggested in HLS . The kernel of second integration has a removable singularity at and can be evaluated via alternating point quadrature rule.

The overall discretization of the integral equation gives rise to a dense system of linear equations comprising of equations, where is the number of connected components of and is the number of marker points on the boundary of each component. We have an additional unknown in the form of . We solve this system using an iterative GMRES SS technique. The GMRES requires only the (dense) matrix-vector multiplication routine and this is the most time consuming part of the iterative solver. Since our matrix is dense, the routine is completed by operations. The cost of matrix-vector multiplication operation can be reduced by the application of a parallel matrix-vector multiplication. It can also be reduced to by the use of fast summation algorithms Greengard1 ; feng2014parallel ; Lindsay1 . We do not use any preconditioner in the solver.

### Evolution of domain interfaces

The discretization of the integral equation gives rise to a stiff system of ODEs as the motion of the interface is curvature driven HLS . The time explicit methods result in a stability constraint where is the spatial resolution. Moreover, the Lagrangian marker points can come close to each other during the course of evolution. To circumvent these problems, we implement the small scale decomposition technique due to Hou et. al. HLS . This special temporal scheme reduces the stiffness requirement to . The scheme also prevents two points from coming too close to each other by distributing the markers on the interface using equal arc length frame and then maintaining the same at all time by the addition of a tangential velocity at every step of calculation.

#### Dynamics of the interface

Once the velocity is obtained for each marker point, we do not update Eq. (19) directly. Instead, the dynamics of the problem is recast in terms of the lengths of the interfaces and the angle that the tangent to the marker point makes with the positive -axis. First, we add a tangent velocity to the interface where is given by

 T(α,t)=T(0,t)+∫α0s′ακ′Vdα′−α2πκ′Vdα′. (47)

After adding the tangential velocity, the motion of the interface is given by

 ddtx(α,t)=V(α,t)n+T(α,t)s. (48)

The addition of the tangential velocity does not change the shape of the interface; however, it is crucial for maintaining the equal arc length distribution of the marker points throughout the computation and prevents the clustering problem. Once the equal arc length distribution is taken care of, we pose the dynamics of the problem with the following two equations,

 Lit =∫2π0θiαVi(α,t)dα, (49) θit =2πLi(−Viα+Tiθiα),i=1,…,M. (50)

The subscripts and denote derivatives with respect to these variables. We use an additional superscript to indicate the interface for which the equations are written. We obtain one equation for for each of the domains, while we get one equation for for every marker point on the boundaries of the domains. Thus, we must solve ordinary differential equations in total. It should be noted that the interface can be fully recovered from and by integrating the relation

 (51)

#### Small-scale decomposition and updating the interface

The stiffness of the original problem propagates to Eq. (50), while Eq. (49) is non-stiff. The latter can be integrated explicitly, but the solution technique for the -equation is far from trivial. This equation is solved using small-scale decomposition (SSD), an idea which has been successfully used in a number of problems in the domain of, e.g., HS flow, micro-structure evolution Shuwang2 ; leo2000microstructural , vesicle wrinkling sohn2010dynamics , and dynamics of an epitaxial island li2011boundary . In problems driven by Laplace-Young boundary conditions, the critical factor in the numerical computation is the curvature of the interface. It introduces higher derivatives in the dynamical equations and results in severe stability constraints. For example, the analysis of the equations of motion reveals HLS that, at small spatial scales, where denotes the periodic Hilbert transform of and therefore Eq. (50) becomes

 θt=σs3αH[θααα]+N(α,t), (52)

where the term In the last equation and the subsequent ones, we suppress in the superscript to keep our notation simple, but its presence should be understood. SSD reveals that the part gives rise to a stiffness condition . The same analysis shows that the term is non-stiff.

We identify that in Fourier space, the dominant term on the right hand side of the Eq. (52) diagonalizes and the equation becomes

 ^θt=−σ|n|3s3α^θ(k,t)+^N(k,t). (53)

We time-integrate the -equation in Fourier space with a semi-implicit time-stepping algorithm HLS . Using an integrating factor, we obtain

 ddt⎛⎜⎝e−σ|n|3s3α^θt⎞⎟⎠=e−σ|n|3s3α^N(k,t). (54)

Then, we use a second-order Adams-Bashforth (AB2) method to discretize Eq. (54) as

 ^θn+1(k) =ek(tn,tn+1)^θn(k) +Δt2(3ek(tn,tn+1)^Nn(k)−ek(tn−1,tn+1)^Nn−1(k)), (55)

where the subscript/superscript denotes numerical solution at and we define

 ek(tn,tn+1)=exp(−σ|k|3∫tn+1tndts3α(t)). (56)

To evaluate the term , we first integrate the non-stiff Eq. (49) using AB2 which gives

 Ln+1=Ln+Δt2(3Mn−Mn−1), (57)

with . Also, , and we apply the trapezoidal rule to evaluate integrals in and as

 ∫tn+1tndts3α(t) ≈Δt2⎛⎝1(snα)3+1(sn+1α)3⎞⎠, (58) ∫tn+1tn−1dts3α(t) ≈Δ⎛⎝12(sn−1α)3+1(snα)3+12(sn+1α)3⎞⎠. (59)

The AB2 method depends on two previous values, and therefore, we initiate the computation at time using Euler’s method to obtain the relevant quantities at . In the subsequent time-steps, we use the AB2 method as two previous time-step values are always known. The accumulation of noise is a problem JLL ; therefore, we employ a cutoff filter to prevent the accumulation of round-off error Krasny1 and a 25th-order Fourier filter to damp the higher, nonphysical modes and suppress the error due to aliasing.

## 6 Numerical Results

In this section, we discuss the results of our numerical simulations. We first compare the results of nonlinear simulation with linear analysis and then demonstrate the spatio-temporal accuracy of our code. Finally, we compute several interesting cases where the domain has different initial configuration. In all our simulations, we set the surface tension parameter to .

### 6.1 Comparison of results of linear analysis and nonlinear simulation

The evolution of a perturbed circular interface is investigated, with the initial interface at given by

 R+δcos4θ=2+0.01×cos4θ, (60)

and we choose . The simulation is carried out up to a time . Evolution of and against time are shown in Fig. 2, using results from the nonlinear simulation and the linear analysis (Eqs. (31) and(32)). The plots indicate excellent match between the two in the beginning thus validating our numerical methods. Once becomes large, we observe disagreement between the results of the linear analysis and the nonlinear simulation, especially in the evolution of . It is evident from the plots that the linear system over-predicts the growth of the mode. This simulation confirms that the linear solution holds for a short time span and the fully nonlinear simulation is needed to predict the evolution over a longer time.

Fig. 3 shows the evolution of the interface, where the innermost contour corresponds to the shape at . For all simulations up to this point, we used a GMRES tolerance of . The filters are also set to this tolerance.

### 6.2 Spatio-Temporal convergence

Figs. 4(a) and 4(b) show the spatio-temporal accuracy of our numerical simulation using initial shape defined in Eq. (60) and with other parameters unchanged. Note that our numerical method is spectrally accurate is space and second-order accurate in time. In Fig. 4(a), we demonstrate the spectral accuracy of our code by plotting the maximum of

 −log10∣∣x(t,N)−x(t,Nf=1024)∣∣

for values , and at time . is chosen so that the results are very accurate in time. Observe that even with , the results match up to . This indicates very a rapid decay of error with and confirms the spectral accuracy of our code.

In Fig. 4(b), we plot the maximum of for and three values of ,