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 a domain , is the density difference, , at position and at time , where the overbar denotes the average of a quantity, e.g.
is given by the solution of the Poisson problem,
where the last condition is introduced to enforce the uniqueness of . Here, we use for the double-well free energy the form
which has two minima at . The chemical potential is obtained by the first variation of the functional
|which yields the flux|
|The system is closed via mass conservation|
|together with boundary and initial conditions|
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
and thus the corresponding phase-field model
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 ,
Similarly for ,
On the fixed boundary , the rescaled boundary conditions are
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
where is the distance along the inward normal direction from the sharp interface , given by
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
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
Finally, the value of can be expressed as
For the derivation of the boundary integral formulation, it is convenient to reformulate the sharp-interface problem in terms of the variable
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:
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
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:
We also transform the far-field boundary condition Eq. (14c) to
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
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
Similarly, in the exterior domain, the solution of the boundary value problem of the Poisson equation in coordinates is given by
Equating the two gives an additional relation between the radii of the interior and the total domain,
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
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
where is the zeroth order solution and is the first order solution. A straightforward computation yields the zeroth order solution as
Next we compute the first order corrections and in this case, is of the form where
The function is of the form where
Once the functions and are available up to first order, we may proceed to calculate the velocity of the interface as
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
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.
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
where the function is analytic and -periodic in the parameter . The local tangent and the normal vectors to the interface are
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.,
for As , we have
Similarly for the exterior problem,
for , where is an unknown to be solved. As , we have
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
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
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).
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 .
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
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
After adding the tangential velocity, the motion of the interface is given by
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,
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
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
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
We time-integrate the -equation in Fourier space with a semi-implicit time-stepping algorithm HLS . Using an integrating factor, we obtain
Then, we use a second-order Adams-Bashforth (AB2) method to discretize Eq. (54) as
where the subscript/superscript denotes numerical solution at and we define
To evaluate the term , we first integrate the non-stiff Eq. (49) using AB2 which gives
with . Also, , and we apply the trapezoidal rule to evaluate integrals in and as
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
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
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 ,