As soft condensed matter, amphiphiles are known to form various structures in aqueous environments. An amphiphilic molecule usually consists of a hydrophilic head and a hydrophobic tail connected by a covalent bond. Consequently, amphiphiles spontaneously arrange themselves in water in such a way that the hydrophobic tails are segregated from water, protected by the hydrophilic heads. Soft matter systems tend to self-assemble into lower-dimensional structures such as surfaces, curves and points, giving rise to sheet-like membranes, polymer networks and colloidal dispersions, respectively [1, Pages 107 and 108]. One particularly important example is the bilayer membrane formed by lipids in water, which exhibits both rigidity and fluidity in that the membrane resists deformation while allowing rapid lateral diffusion of lipid molecules within each monolayer. The elasticity of the membrane is very different from those of solid materials such as aluminum foil and plastic film. The membrane is soft, which is a crucial property for biological cells and artificial liposomes. The typical energies required to bend a membrane are small enough for thermal fluctuations at room temperature to be important [2,Pages 3 and 4]. Indeed, the bending elasticity of membranes is only a high-order effect, as will be discussed in Section 3.2.
At the macroscopic level, the Helfrich energy introduced in 1973 proved to be a successful continuum model for describing the elasticity of the lipid bilayer membrane (see Appendix A). In this model, the membrane is treated as a two-dimensional surface of zero thickness, with its energy given by the surface integral of a quadratic function in the principle curvatures. However, the actual membrane is of a bilayer structure and nonzero thickness (usually a few nanometers). In order to gain a detailed knowledge at the microscopic level, an atomistic molecular dynamics simulation was carried out in 1992 [3]. Albeit accurate, such a simulation was time-consuming and thus restricted to a relatively small spatio-temporal scale (60 nm 3 × 0.2 ns), rendering the physical processes of interest out of reach. 1 Since our primary interest is not in individual atoms but rather in the collective behaviors of large numbers of atoms, it is natural to group several neighboring atoms into a single bead, leading to the so-called coarse-graining methods, which reduce the degrees of freedom and accelerate the computation [2,[5][6][7].
On an even coarser scale, a smooth density function is often used to represent the spatial distribution of each type of atoms or atom groups at the mesoscopic level. In 1986, Ohta and Kawasaki derived a density functional theory from statistical physics to explain the mesoscopic periodic patterns formed by diblock copolymers [8]. This theory was later generalized to triblock copolymers [9,10] as well as mixtures of diblock copolymers and homopolymers [11]. The latter generalization was recently shown to be capable of modeling the lipid bilayer membrane along with its fusion processes [12], and we study a special case of such a generalization in this paper.
We focus on the sharp interface limit (also known as the strong segregation limit). In our energy, we let and denote the regions occupied by the hydrophobic tails and hydrophilic heads of the lipids, respectively. The rest of the space is occupied by water. As their names suggest, the hydrophobic tails are insoluble in water, while the hydrophilic heads are assumed to be miscible with water. Therefore, the interfacial tension exists only on the interface of , but not on the -water interface. An additional Coulombic term accounts for the covalent bonding between the hydrophobic tails and hydrophilic heads. This energy has been studied in the small mass regime in connection with the spherical micelle formed by amphiphilic surfactants in water [13]. In the large mass regime, although the question remains largely open, it is believed that the energy minimizers might resemble the lipid bilayer membrane [13,Page 4]. In fact, a variant of this energy (which makes use of the 1-Wasserstein distance) has been proposed to model the lipid bilayer membrane at the mesoscopic level [14,15].
For two subsets , of R satisfying ∩ = ∅, we study the minimization problem of the following energy In this paper we are mainly concerned with the three-dimensional (3-D, i.e., = 3) case which is physically most relevant, and we will also consider the 2-D case ( = 2) for comparison purposes. Note that the 1-D case ( = 1) with ( -) = -| -|∕2 and = 1 has been solved in [16, Section 3]: any local minimizer consists of one or multiple non-overlapping bilayer(s), such as two bilayers 0 (where 0 represents a layer of water of arbitrary thickness, and any layer is twice as thick as any layer), and the global minimizer selects the local minimizer whose layers are of thickness close to 3 √ 12∕ . We expect those 1-D results can be generalized to ≠ 1. Intuitively, we can imagine that and uniformly carry equal amounts of positive and negative charges, respectively, so that the total electrostatic potential energy arising from the electrostatic interactions is given by ( , ). We define the associated electrostatic potential as
According to [13,Equation (2.6)], by noticing -= 1 -1 ∕ , we can rewrite as
⃝ The global minimum of has a lower bound and an upper bound. As → ∞, both bounds scale linearly with [16,Section 4].
Next we give a brief review of existing studies for various .
In the large mass regime ≫ 1, only some qualitative properties of the energy minimizers are known in the literature, e.g., there is a uniform bound on the mean curvature of for the global minimizer for ⩾ 1 [13, Page 9]. The exact global minimizer is unknown and is conjectured in several works (see [13,Figure 2(c) and Page 9] and [17,Bottom of Page 78]) to resemble a planar bilayer membrane, cut off at large distance.
In the literature there is a variant of the energy (1), with the nonlocal term replaced by the 1-Wasserstein distance. For this variant, as → ∞, it is energetically preferable for and to form a closed bilayer membrane with an approximately uniform thickness [14,15]. In [17,Bottom of Page 21] van Gennip mentioned a failed attempt to generalize this result from such a variant to the problem (1).
In the limit of → 0, the problem (1) reduces to the so-called surface charge model [18,Section 6]. For the global minimizer, the negative charge is concentrated on the surface of , and acts like a Faraday shield canceling out the positive charge carried by , i.e., we have = 0 on R 3 ∖ . As mentioned above, for = 1 and ≫ 1 the global minimizer should locally resemble the bilayer membrane shown in Fig. 3-a, where a layer of of approximately uniform thickness surrounds . For ⩽ 1 and ≫ 1, it is natural to expect that the global minimizer takes on a similar bilayer membrane structure, with the thickness of the layer converging to 0 as → 0.
In the limit of → ∞, the mass of becomes infinitely large, and the negative charge carried by becomes infinitely dilute, therefore the problem (1) reduces to the liquid drop model [19] where we only consider the positive charge carried by . It is widely believed that the infimum of the energy in the liquid drop model is attained by a ball for ≪ 1, and is approached by many equally large distant balls for ≫ 1.
Therefore, for ≫ 1, it is natural to expect the global minimizer to be a micelle for ≪ 1, and to be the union of many equally large non-overlapping micelles for ≫ 1.
We now summarize our main contributions and then outline the overall structure of the rest of this paper. As mentioned in Section 1.2, regarding the global minimizer for = 1 and ≫ 1, several works (e.g., [13, Figure 2(c) and Page 9] and [17, Bottom of Page 78]) documented the belief that is approximately a large disk with a radius of order √ and with a thickness of order 1, and that is of an approximately uniform thickness surrounding , as can be seen from our numerical simulations in Fig. 9. However, our numerical results indicate that such a disk-shaped membrane has higher energy than a liposome (see Definition 2.1) for = 1 and ≫ 1. In Fig. 9, the disk-shaped membrane is slightly thicker near its rim, and thus the rim carries an energy penalty of order √ , proportional to the circumference. Therefore for the disk-shaped membrane, the energy-to-mass ratio ∕ converges with order 1∕ √ , which is consistent with [16,Theorem 8]. Meanwhile, according to Corollary 2.2, the convergence rate for the liposome is 1∕ . Note that the infimum of ∕ is attained in the limit of → ∞ [13, Proposition 8.3]. Our work suggests that the liposome is the global minimizer for = 1 and ≫ 1, since a sphere minimizes the Helfrich energy in Proposition 3.10.
In the existing literature, the equal mass case = 1 has been the focal point of most studies (e.g., [13,[16][17][18]20,21], except for a 1-D study [11]). In this paper we consider the general cases ∈ (0, ∞). As mentioned in the paragraph following (1), the 1-D results should be qualitatively the same for different . But as we will illustrate in Fig. 3, as increases, the optimal morphology in 3-D should undergo transitions from the bilayer membrane to cylindrical micelle to spherical micelle. As mentioned in Section 1.2, in the large mass regime, it is natural to expect that for sufficiently large , the global minimizer consists of many spherical micelles which are scattered and non-overlapping, as depicted in Fig. 3-c. Meanwhile, for intermediately large , it is natural to expect the global minimizer to resemble a cylindrical micelle, see Fig. 3-b.
In Sections 2.1 and 3.1, we use asymptotic analysis to study the energy minimizer among radially symmetric liposome candidates for > 0 and ≫ 1. An important finding is that the inner layer is slightly thicker but has slightly less mass compared to the outer layer. Interestingly, we notice that such slight differences also exist in a variant model (at least in 2-D), where the Coulombic nonlocal term is replaced by the 1-Wasserstein distance, as mentioned in Remark 3.6-3 ⃝. If we impose equal mass on the inner and outer layers at the expense of optimality, then the energy in our model will increase on the same order as the bending energy, as mentioned in Remark 3.6-2 ⃝. For a general bilayer membrane which is not necessarily radially symmetric, up to a suitable rescaling, we can consider its Gamma-limit with vanishing thickness as → ∞. Under the conjecture that its energy converges to the Helfrich energy, in Proposition 3.10 we calculate the bending and Gaussian moduli using our asymptotic results for 2-D and 3-D radially symmetric liposome candidates (note that a 2-D circular liposome may be viewed as a 3-D cylindrical bilayer which resembles a very long tube [22, Middle of Figure 1]). Our calculation reveals that for small , the Gaussian modulus is positive, and thus the optimal structure may resemble triply periodic minimal surfaces (TPMS) which exist ubiquitously in copolymer systems and biological specimens (see Remark 3.11).
In Fig. 1 we summarize the conjectured candidates for 3-D global minimizers in the large mass regime, with the parameters to be defined later. For < 1 , we believe that the local structure of the global minimizer resembles a bilayer membrane, and that its global structure is a liposome ( 0 < < 1 ) or approximately a TPMS ( < 0 ). For 1 < < 2 , its local structure resembles a cylindrical micelle, and its global structure resembles a circular ring. For > 2 , the global minimizer consists of many scattered droplets of spherical micelles.
In Section 5, we present numerical evidence for our conjectures. The numerical simulations are based on a phase-field reformulation of the sharp interface model (1). In order to justify such a reformulation, we prove a Gamma-convergence result in Section 4. The novelty of the proof is that our phase-field energy is degenerate: only one of the two order parameters is penalized by the Dirichlet energy, and the potential well has non-isolated minimizers. The rest of this paper is organized as follows. In Section 2, we restrict ourselves to the simplest radially symmetric case and derive some asymptotic results. In order to better present our asymptotic results, in Section 3 we rescale the energy and propose some conjectures. In Section 4 we present our phase-field reformulation and prove its Gamma-convergence to the sharp interface model. In Section 5 we present some numerical simulations. In Section 6, we conclude with remarks about future directions. In Appendix A we provide some background knowledge on the Helfrich and Willmore energies. In Appendix B, we present the detailed calculations in the radially symmetric case. In Appendix C, we recall the results in a variant model where the Coulombic nonlocal term is replaced by the 1-Wasserstein distance, in order to make it convenient for readers to draw comparisons.
In this section, we only consider radially symmetric candidates, thus reducing the variational problem to a finite-dimensional optimization problem which can be analyzed asymptotically. We obtain the asymptotic expansions of the minimum energy and the optimal layer thickness, thus laying the foundation for further studies of non-radially-symmetric cases. In Sections 2.1 and 2.2, we consider two types of candidates, namely liposome and micelle, as depicted in Fig. 2, and then draw comparisons in Section 2.3.
Liposome candidates are concentric rings or shells as shown on the left of Fig. 2. More specifically, we have the following definition.
) for some 0 < 0 < 1 < 2 < 3 , where ( ) is an -dimensional ball of radius centered at the origin .
⃝ satisfies the following mass constraints
With and fixed, as → ∞, the minimizer of ( , ) among the liposome candidates satisfies
with the following asymptotics
)
) , for = 3.
Under the additional assumption that the inner and outer layers have the same mass [17,Equation (5.16)], i.e., 3 -2 = 1 -0 , with and fixed, as → ∞, the minimizer of ( , ) among the liposome candidates still satisfies (6), with the following asymptotics
)
) , for = 3.
Remark 2.4.
⃝ The additional assumption in Corollary 2.3 was used in [17,Pages 83 and 84]. Albeit not the optimal choice, it was thought to be only slightly non-optimal. Our calculations show that the leading-order term of the energy remains the same under this additional assumption, but that the next-order term (corresponding to the bending energy) becomes 6 to 21 times as large, depending on and . In other words, the bending energy can be decreased by dropping this additional assumption and allowing the inner and outer layers to have different masses.
Q. Du et al.
⃝ From (6) we know that +1 -converges to a constant as → ∞. This can serve as an illustration that the lipid bilayer has an intrinsically preferred thickness, which is consistent with [16,Remark 4]. In the large mass regime ≫ 1, the bilayer should be of thickness of order 1.
Definition 2.5. In Definition 2.1, if 0 = 1 = 0, then ( , ) is called a micelle candidate. (Note that for = 2, in some literature [24] the terminology core-shell is used instead.) Proposition 2.6. The energy (1) of a micelle candidate is ( , ) = Per + ( , ), where
Proof. Proposition 2.6 is a consequence of Proposition B.1 with 0 = 1 = 0 and
. □ Corollary 2.7. For a micelle candidate, if = 2, then ∕ attains its minimum at = 4
, and min
, and min
Remark 2.8.
⃝ Due to the screening property (see Proposition 1.1-3 ⃝), for the disjoint union of non-overlapping components, the energy is simply the sum of the energy for each connected component. Therefore, in the large mass regime ≫ 1, we can construct a candidate consisting of many non-overlapping micelles, each of which has a mass close to the optimal mass given in Corollary 2.7 (similar to the construction mentioned in [16,Page 96]). Such a candidate will attain the optimal energy-to-mass ratio asymptotically as → ∞. In fact, our energy in Proposition 2.6 is of a similar form to [25,Equation (6.1)]. Using the same approach as [25, Proof of Lemma 6.2], we can prove that the optimal way to allocate the mass is as follows: let (assuming > 0) be the mass of each connected component, then all but one (say, { } ≠1 ) must be equal, with 1 = (1). Therefore, in the large mass limit → ∞, the optimal way to allocate the mass yields asymptotically the same energy-to-mass ratio as Corollary 2.7.
⃝ According to [16,Theorem 8], we can regard a 2-D micelle as the limiting case of an infinitely long 3-D cylindrical micelle, with asymptotically the same energy-to-mass ratio.
⃝ We have also considered the candidates with 0 = 0 and 1 > 0, and they always have higher energy-to-mass ratios than the liposome candidates in the large mass limit, according to our calculations. Therefore, they are omitted in our discussions.
We now compare the candidates that are considered in Sections 2.1 and 2.2. By comparing the leading-order term of ∕ in Corollaries 2.2 and 2.7, we are led to the belief that in the large mass regime in 3-D, as increases, the preferred morphology should be successively bilayer membrane, cylindrical micelle, and spherical micelle, as shown in Fig. 3. This picture is qualitatively consistent with the predictions by other theories [23,Section 1.3.2.3]. Conjecture 2.9. Let = 3. For any fixed > 0 and > 0 (where 0 is defined in Remark 3.11-2 ⃝, for reasons that will be mentioned therein), we have inf , ( ,
Here 1 ≈ 1.81696 is the unique nonzero root of + 2 ∕3 = ( + 1) ln( + 1), and 2 ≈ 3.64572 is the unique nonzero root of 5 ( ( +1) ln( +1) -
) .
Remark 2.10.
⃝ The piecewise function ( ) is plotted on the left of Fig. 3, where 1 and 2 are also marked. As increases, the minimum energy is attained on different branches indicated by different colors, and the preferred morphology should transition from bilayer membrane to cylindrical micelle to spherical micelle, which are illustrated by Fig. 3-a, b and c, respectively.
⃝ In [17, Section 5.3], van Gennip also compared the energy among the bilayer membrane, cylindrical micelle, and spherical micelle candidates. However, van Gennip only considered the case of = 1, and concluded that the bilayer membrane has the lowest energy among all the candidates. Our calculations reveal that cylindrical micelle and spherical micelle may have lower energy for ≠ 1, and that amphiphiles can self-assemble into not only sheet-like membranes, but also polymer networks and colloidal dispersions.
In Section 2, we have briefly presented some asymptotic results as → ∞. In order to make it convenient and rigorous to translate our asymptotic results within the Gamma-convergence framework, in this section we rescale the energy functional so that the minimizer converges in Radon measure to lower-dimensional structures (surfaces, curves and points). In Section 3.1, we elaborate on the asymptotic results of the liposome candidates for the rescaled energy, and in Section 3.2, we propose some conjectures about Gamma-convergence of the rescaled energy. Definition 3.1. Analogous to [14], we define a rescaled version of (1):
where can be chosen from {1, 2, 3} as needed, = supp( ) and = supp( ). Moreover, Proof. Use the scaling properties of the perimeter and Newtonian potential. □ Remark 3.3. According to Proposition 3.2, is equivalent to up to a rescaling, and thus the results in Proposition 1.1 also apply to : for any , , > 0 and for any , the global minimizer of exists and satisfies the screening property.
Throughout this subsection, we let = 1. We consider liposome candidates given by = ( 2 )∖ ( 1 ) and =
where ( ) is an -dimensional ball of radius centered at the origin (see also Definition
and has the following
) .
Corollary 3.5 (of Propositions 3.2 and B.4). Under the additional assumption that the inner and outer layers have the same mass [17,Equation (5.16)], i.e., 3 -2 = 1 -0 , let = 1, with and fixed, as → 0, the minimizer of among the liposome candidates has the following asymptotics if = 2,
and has the following asymptotics if = 3,
Remark 3.6.
⃝ Due to the curvature of the bilayer, the inner layer is slightly thicker than the outer layer (i.e., 1 -0 > 3 -2 ), and is penalized for bending on the second order (i.e., the 2 term corresponding to the bending energy).
⃝ Under the equal mass assumption 3 -2 = 1 -0 , the second-order term of is 6 to 21 times as large as that of the optimal liposome candidate whose inner layer has slightly less mass than the outer layer.
⃝ Under the equal mass assumption, the difference in the thickness between the inner and outer layers is asymptotically three times that of the optimal liposome candidate, in both 2-D and 3-D. This relation is also true at least in 2-D even if the Coulombic nonlocal term is replaced by the 1-Wasserstein distance (see Appendix C). It is therefore natural to wonder how universal this relation can be.
We now propose some conjectures about the Gamma-limits of the energy functional ′ ∶= ( -( )
∕ 2 as → 0. Throughout this subsection, is chosen to be the codimension of the expected geometry of the global minimizer as shown in Fig. 3:
where 1 and 2 are defined in Conjecture 2.9.
Conjecture 3.7. Let = 2 and → 0. For ∈ (0, 1 ), the Gamma-limit of ′ is the elastica functional defined for closed 2 , 2 curves in R 2 in the sense of Radon measure, similar to the elastica functional W mentioned in [14,Theorem 4.1]. For ∈ ( 1 , ∞), the Gamma-limit of ′ is a mass partition functional defined for weighted Dirac delta point measures, similar to the mass partition functional E 2d 0 mentioned in [25,Theorem 6.1] (see also the mass partition functional 0 mentioned in [24,Theorem 4.2]). Conjecture 3.8. Let = 3 and → 0. For ∈ ( 0 , 1 ), the Gamma-limit of ′ is a quadratic form in the principal curvatures defined for closed surfaces in the sense of Radon measure, similar to [15,Conjecture 2.4]. For ∈ ( 1 , 2 ), the Gamma-limit of ′ is the elastica Q. Du et al. functional defined for closed 2 , 2 curves in R 3 in the sense of Radon measure, which is a 3-D generalization of [14,Theorem 4.1]. For ∈ ( 2 , ∞), the Gamma-limit of ′ is a mass partition functional defined for weighted Dirac delta point measures, similar to the mass partition functional E 3d 0 mentioned in [25,Theorem 4.3]. If the first statement in Conjecture 3.7 and the second statement in Conjecture 3.8 are correct, then Proposition 3.9 tells us that the global minimizer should resemble a circle. Proposition 3.9. Let be a closed 2 , 2 curve in R 3 with a prescribed length, and let be its curvature, then the elastica functional ∫ 2 d is minimized when is a circle.
Proof. The length of is given by ∫ 1d and is assumed to be fixed. By Cauchy-Schwarz inequality we have ∫
, where the second inequality is due to [26,Theorem 7.2.3] and becomes an equality if and only if is a planar convex curve. □ Proposition 3.10. Assuming the first statement in Conjecture 3.8 is correct, i.e., for = 3 and ∈ ( 0 , 1 ), the Gamma-limit of ′ is of the following quadratic form (cf. Appendix A for Helfrich and Willmore energies)
where = ( 1 + 2 )∕2 (with 1 and 2 being principal curvatures, positive if is a sphere), = 1 2 , and 1 , 2 ∈ R. If the recovery sequence is given by a bilayer with the middle layer of uniform thickness (similar to the construction in [15,Section 4]), then we have
Proof. According to Proposition B.2 and the Proof of Proposition B.1, the liposome candidates shown in Corollary 3.4 satisfy the screening property. For any fixed , there exists a unique satisfying the screening property [18, Remark 4.2], which is a necessary condition for a minimizer (Proposition 1.1-3 ⃝). Therefore, for the Radon measure of a sphere in 3-D, the recovery sequence is by assumption given by radially symmetric liposome candidates in 3-D. According to [16,Theorem 8], the 2-D liposome candidate can be regarded as a 3-D cylindrical bilayer which resembles a very long tube [22, Middle of Figure 1]. For a sphere in 3-D with radius , we have = 1∕ 2 and = 1∕ . According to Corollary 3.4, we have
For a cylinder in 3-D with length and radius , we have = 0 and = 1∕ (2 ). Similarly, we have
) .
In Corollary 3.4, the radius of the 2-D liposome candidate can be approximated by 4 ). Therefore the coefficients 1 and 2 can be determined. □ Remark 3.11.
⃝ By the Gauss-Bonnet formula, for closed surfaces in the same homotopy class (i.e., closed surfaces of the same genus ), we have ∫ d = 4 (1-), therefore ( 8) can be reduced to the Willmore energy ∫ 2 d as long as no topological change occurs.
2
We are now ready to explain the requirement > 0 imposed in Conjectures 2.9 and 3.8. If such conjectures were also true for < 0 , then we would still be able to prove (9) in Proposition 3.10. However, for < 0 we have 2 > 0. Therefore, the quadratic form in Proposition 3.10 is not positive semi-definite, and it is energetically more favorable to have zero mean curvature and negative Gaussian curvature. As pointed out by [27, Page 143], if 2 becomes positive, then it is preferable for a surface to have large genus , or many ''holes''. For example, triply periodic minimal surfaces shown in Fig. 4 possess many ''holes''. Since the quadratic form is indefinite, it is desirable for those triply periodic minimal surfaces to have infinitesimal lattice constants. However, we expect an equilibrium to be attained at a finitely small lattice constant (on the order of ) due to counteracting higher-order terms. ⃝ It surprised us to find that 2 can be positive for relatively small . However, with hindsight we may come up with the following intuitive but non-rigorous explanation. For relatively large , bending of the bilayer is penalized because there will be a slight difference in the thickness between the inner and outer layers, as mentioned in Remark 3.6-1 ⃝. For relatively small , such penalty may be relatively small because the layers have vanishing thickness as → 0. Meanwhile, the middle layer may prefer saddle-splay deformations, because the layer becomes more ''spread out'' in this way and thus the Coulombic repulsion within the layer may decrease. ⃝. More specifically, our numerical results indicate that a gyroid-like minimizer has lower energy-to-mass ratio than the planar bilayer for = 0.6 < 0 , and thus the latter cannot be a global minimizer for ≫ 1. For = 1 > 0 , our numerical results suggest the opposite. 5 ⃝ Various cubic bicontinuous structures resembling triply periodic minimal surfaces can be observed in copolymer systems and biological specimens (see [28,29] and [30, Section 4.1]), e.g., in the endoplasmic reticulum, Golgi apparatus and mitochondria [1, Bottom of Page 144].
Our finding provides a plausible explanation for such phenomena in the parameter regime < 0 , and also demonstrates the ability of the Ohta-Kawasaki energy to capture such aspects of amphiphile self-assembly. Remark 3.12.
1 ⃝ Conjectures 3.7 and 3.8 states that our nonlocal problem Gamma-converges to a local problem. Our rationale behind those conjectures is as follows. According to the screening property (Proposition 1.1-3 ⃝), each connected component of the minimizer satisfies the charge neutrality condition, and the total energy is just the sum of the energy for each connected component. Therefore we can consider each connected component separately. For each connected component, although the problem is nonlocal, the screening property may cause such nonlocality to decay sufficiently fast as → 0.
2 ⃝ Analogously, the screening property is also the reason that the electrostatic potential energy of a crystal is an extensive property. The electrostatic interactions are nonlocal, but an extensive property is local, i.e., it is proportional to the size of the crystal. In a variant model where the nonlocal term is given by the 1-Wasserstein distance, a similar phenomenon occurs: the nonlocal (or global) problem converges to the elastica bending energy which is a local problem [14, Section 9.3].
In order to provide some evidence for our conjectures in Section 3.2, we propose a phase-field reformulation which will be used for the numerical simulations in Section 5. We also prove the Gamma-convergence of our phase-field reformulation to a Gamma-limit, which is shown to be equivalent to the original problem (1).
The phase-field method [31] is a useful tool to study the motions of interfaces. The basic idea is to use a narrow but diffuse interface in place of the sharp interface, and the thickness of the interfacial layer is controlled by a diffuseness parameter . The interface is implicitly given by the level set of a smooth function, so there is no need to explicitly track the interface. Our phase-field reformulation is very similar to a previous work by two of us [19, Section III.A].
On a bounded domain ⊂ R with | | ≥ + , we define the following phase-field functional
where P is the diffuse interface version of the perimeter
with the diffuseness parameter > 0, and the double-well potential given by (see Fig. 5)
which is to penalize violations of the three conditions:
Since |∇ | 2 is not penalized in our energy, we choose a degenerate well for instead of the classical double well ( -2 ) 2 , in order to prevent from being trapped in the local minima 0 and 1.
We choose a function which strictly increases on [0, 1] satisfying (0) = 0 and (1) = 1. The simplest choice is the identity function ( ) ∶= , which will be used in this section for our proof of the Gamma-convergence. However, for our simulations in Section 5, we will choose to be a nonlinear function in order to achieve the numerical efficiency as mentioned in [19, Section III.A]. As we will see later in Propositions 4.3 and 4.14, if ( , ) is a minimizer of E for ≪ 1, then and should be approximately 0 or 1 at most places in . Therefore, the specific choice of should have a vanishing effect on the minimizer of E as → 0.
With penalty coefficients 1 , 2 > 0, the mass constraint term C is defined as follows
)
The nonlocal term N is defined as follows
where the nonlocal kernel is the Green's function of negative Laplacian -. We consider two types of boundary conditions for . The first type is no boundary conditions, under which is the Green's function in the free space (i.e., the so-called fundamental solution). The second type is periodic boundary conditions, in which case we require to be a rectangle ( = 2) or a rectangular cuboid ( = 3), and according to [25, Equations (2.1) and (2.2)],
is the periodic extension of the fundamental solution plus a continuous function, satisfying ∫ = 0. In Section 4.3.2, we justify that the above two types of boundary conditions are equivalent. In fact, for = 1 such equivalence has already manifested itself in the striking similarity between Theorems 1 and 2 in [16] (note that ⋯ therein is equivalent to ⋯ , since an infinitesimal block can be appended to the end without energy penalty, as long as -0 interfaces are not penalized).
In this subsection, we aim to prove the Gamma-convergence of E as → 0 and 1 , 2 → ∞ (where the subscript is used to emphasize the dependence of E on ).
As mentioned before, in our proof of the Gamma-convergence, we let be the identity function. As a first step, we let 1 , 2 → ∞, so that the admissible function space of E is restricted to
, and thus the constraint term C vanishes. For convenience, in most parts of our proof of the Gamma-convergence, we do not explicitly incorporate the mass constraints, because it should be straightforward except in the proof of the limsup inequality (Proposition 4.7), where we do provide more details of the mass constraints. In addition, we have
where = * ( -∕ ), i.e., (⃗ ) = ∫ (⃗ -⃗ )
for any ⃗ ∈ . If we impose no boundary conditions on (i.e., is the fundamental solution), then satisfies -= -∕ ; if we impose periodic boundary conditions on , then satisfies -= constant + -∕ and ∫ = 0, where the constant is zero as long as ∫ = ∫ . Q. Du et al. We assume that ⊂ R is an open bounded set with Lipschitz boundary. We extend the domain of definition of E to 2 ( ) × 2 ( ):
We want to prove that E Gamma-converges to the following functional defined on 2 ( ) × 2 ( ): Our proof of the Gamma-convergence of E to E 0 follows closely that of the classical Cahn-Hilliard energy functional given in [32]. The novelty of our proof lies in the degeneracy of the potential shown in Fig. 5, i.e., has non-isolated minimizers ( , ) ∈ {0} × [0, 1]. Furthermore, |∇ | 2 is penalized in our energy while |∇ | 2 is not. Sternberg [33, Section 2] also considered the case where has non-isolated minimizers, but both |∇ | 2 and |∇ | 2 are penalized there. In the rest of this section, we consider an arbitrary sequence of diffuseness parameters { ∶ ∈ N} ⊂ R >0 converging to 0. For brevity we abbreviate the subscript and write E , which should be understood as E for some ∈ N. We also write lim →0 E instead of lim →∞ E .
Definition 4.1. A sequence of functions { ∶ ∈ N} ⊂ 1 ( ) is called equi-integrable (i.e., uniformly integrable) if: for any > 0, there exists an > 0, such that
Remark 4.2. In Definition 4.1, note that { 1 } is always equi-integrable. Indeed, we can approximate 1 using a simple function so that ∫ | 1 -| < ∕2, and then choose such that
Similarly, any finite subset of { ∶ ∈ N} is always equi-integrable, therefore we have an equivalent definition:
then there exists a subsequence { ( , ) } (for brevity, this subsequence is not relabeled) and ( 0 , 0 ) ∈ ( ; {0, 1}) × 2 ( ; [0, 1]) such that → 0 in 2 ( ), ⇀ 0 in 2 ( ), 0 0 = 0 a.e.
Step 1: We show that 2 + 2 is uniformly bounded in 1 ( ) and equi-integrable. According to (10), there exists > 0 such that ( , ) ≥ 2 + 2 whenever 2 + 2 ≥ . For any ≥ , we have
Step 2: We show that a subsequence of { } converges pointwise a.e. to some 0 ∈ ( ; {0, 1}) (by mimicking [32]). Denote According to [34, Theorem 5.5], there exists a subsequence { } (not relabeled) and some 0 ∈ ( ) such that → 0 in 1 ( ). A further subsequence { } (not relabeled) converges pointwise a.e. to 0 . Since ( ) is strictly increasing and continuous in , its inverse -1 is continuous. Define 0 = -1 ( 0 ), since = -1 ( ), we know that converges pointwise a.e. to 0 . Moreover, by Fatou's lemma [35,Lemma 1.83
. Therefore 1 ( 0 ) 2 = 0 a.e., that is, 0 ∈ {0, 1}
a.e. So 0 ∈ { (0), (1) } a.e. Since 0 ∈ ( ), we can write 0 = (0) (1-1 ) + (1) 1 with a set ⊂ of finite perimeter. Hence, 0 = 1 belongs to ( ; {0, 1}).
Step 3: We show that the subsequence of { } obtained in Step 2 converges to 0 in 2 ( ). This is guaranteed by Vitali's convergence theorem [35, Theorems 2.22, 2.24 and 2.29], in view of Step 1.
Step 4: We show that a subsequence of { } converges weakly in 2 ( ) to some 0 ∈ 2 ( ). According to Step 1, { } is uniformly bounded in 2 ( ). By [35,], it has a subsequence (not relabeled) such that ⇀ 0 in 2 ( ) for some 0 ∈ 2 ( ).
Step 5: We show that 0 ∈ 2 ( ; [0, 1]) and 0 0 = 0 a.e. In view of Steps 2, 3 and 4, we can assume that the sequence { ( , ) } satisfies ⇀ 0 in 2 ( ), → 0 in 2 ( ), and → 0 a.e. Our goal is to prove ∫ ( 0 , 0 ) = 0. Since ( , ) is continuous in , by Fatou's lemma, we have the following
According to (10), ( , ) is convex in , and we know
) , therefore we have
where the first summand on the right-hand side can be bounded by E ( , ) , and the second summand can be split into
where the first summand converges to 0 because ( 0 , 0 ) ∈ 2 ( ) and ⇀ 0 in 2 ( ), and the absolute value of the second summand can be bounded by
Since { } is uniformly bounded in 2 ( ), we know ‖ 0 -‖ 2 ( ) is uniformly bounded. Moreover, we have
To summarize, we have proved ∫ ( , 0 ) → 0 as → 0. Therefore ( 0 , 0 ) = 0 a.e. □
Proposition 4.4. For any 0 , 0 ∈ 2 ( ) and { }, { } ⊂ 2 ( ) such that both → 0 and ⇀ 0 in 2 ( ), we have
Proof. We assume that the above right-hand side is finite (otherwise there is nothing to prove). Without loss of generality (by extracting a subsequence if necessary), we further assume that the lim inf is actually a limit. By Proposition 4.3, we know 0 ∈ ( ; {0, 1}), 0 ∈ 2 ( ; [0, 1]), and 0 0 = 0 a.e. We can also assume E ( , ) < ∞ for all , and thus { } ⊂ 1,2 ( ).
Step 1: We first consider the case where = 0. Denote
Noticing 0 = (1) 0 and (1) = 1∕6, we obtain
Step 2: We now consider the case where > 0. We want to prove the following
.
Therefore we only need to estimate ∫ (⃗ -⃗ ) 2 d ⃗ .
Step 1: We first consider the case where is the fundamental solution with = 3. In such case, we have
Inspired by [36, Page 159], we denote
, where (⃗ ; ) is a ball of radius centered at ⃗ . Using the monotonicity of , we obtain for any ⃗ ∈ ,
Step 2: We now consider the case where is the fundamental solution with = 2. In such case, we have
Notice that for any ∈ R >0 , we have
where the first summand monotonically decreases and can be bounded in a similar manner to Step 1, while the second summand monotonically increases and can be bounded in terms of ln sup
uniformly bounded in ⃗ .
Step 3: We now consider the case where is equipped with periodic boundary conditions. In this case (⃗ -⃗ ) is periodic, therefore
According to [25, Equations (2.1) and (2.2)], is the sum of the fundamental solution and a regular part, therefore ∫
) has an upper bound independent of , because ⇀ 0 in 2 ( ). In addition, we have
where the last equality is again due to ⇀ 0 in 2 ( ), and the second equality is due to the Dominated Convergence Theorem. To see why the Dominated Convergence Theorem applies, notice that according to Lemma 4.5, | 0 (⃗ )| is uniformly bounded in ⃗ , and that | (⃗ )| is uniformly bounded in ⃗ and , because ‖ ‖ 2 ( ) is uniformly bounded in . □
Proposition 4.7. For any 0 , 0 ∈ 2 ( ), there exist { }, { } ⊂ 2 ( ) such that both → 0 and → 0 in 2 ( ), satisfying
Proof. We assume that the above left-hand side is finite (otherwise there is nothing to prove), therefore we have 0 ∈ ( ; {0, 1}), 0 ∈ 2 ( ; [0, 1]), and 0 0 = 0 a.e.
Step 1: We first consider the case where = 0 and = 1. We further assume = [-1, 1] and 0 = 1 [0,1] . Therefore 0 ( ) = 0 for ∈ [0, 1] a.e. We take
By Dominated Convergence Theorem, we have both → 0 and → 0 , pointwise a.e. and thus in 2 ( ). We can compute
In addition, we have
where the first summand has the following integral
and the second summand can be bounded as follows
where the inequality is due to 0 ≤ 1 a.e. Since increases monotonically, we have
To summarize, we have
where the right-hand side converges to 1 as → 0. Therefore E 0 ( 0 , 0 ) = 1 ≥ lim sup →0 E ( , ).
So far in our proof of the Gamma-convergence, we have not explicitly dealt with the mass constraints, which should be straightforward because the weak convergence in 2 ( ) implies the convergence of the Lebesgue measure. However, we need to provide more details here. To make sure ∫ = ∫ 0 , we rescale defined above, i.e., we redefine
In this way, we have
. Therefore, after the above modification, satisfies the mass constraint while E ( , ) increases only by order 1∕3 .
Step 2: We now consider the case where = 0 and > 1. We assume 0 = 1 for some open set ⊂ R with being a nonempty compact hypersurface of differentiability class 2 . We further assume H -1 ( ∩ ) = 0, where H is the Hausdorff measure. We can rewrite 0 = 0 ( ) where 0 = 1 [0,∞) , and is the signed distance function defined by
We take = ( ) and = ℎ ( ) 0 , where ( ) ∶=
By Lemma 4.8, we have
Step 1, we know ∫ 1)
We now show both → 0 and → 0 in 2 ( ). To this end, we use |∇ | = 1 and Lemma 4.8 again to obtain
→ 0 pointwise a.e., by Dominated Convergence Theorem, we have ∫ ( -0 ) 2 → 0. Similarly, we have ∫ ( -0 ) 2 → 0.
In order to guarantee the mass constraint ∫ = , we can make some technical modifications to the above in a similar way to the classical Cahn-Hilliard energy functional [32,Equation (26)]. Similar to Step 1, we can also modify to guarantee the mass constraint ∫ = .
Step 3: We now remove the regularity assumption imposed on in Step 2. By Lemma 4.9, there exists a sequence of open sets with being a nonempty compact hypersurface of differentiability class 2 and satisfying H -1 ( ∩ ) = 0, such that → 1 in 2 ( ) and Per → Per . According to Step 2, for each fixed we can find a sequence { , } ⊂ 1,2 ( ) and { , } ⊂ 2 ( ) such that both , → and , → (1-) 0 in 2 ( ) as → 0, in addition to
Therefore we have
By Lemma 4.10, we have a diagonal sequence { , } and { , } such that both , → 0 and , → 0 in 2 ( ), and
Step 4: We now the case where > 0. According to Step 2 in the proof of Proposition 4.4, we have lim →0 N( , ) = N( 0 , 0 ). □ Lemma 4.8. Coarea Formula for Lipschitz Functions [32, Theorem 1.14]. On an open set ⊂ R , let ∶ → R be Lipschitz continuous and let ∶ R → R be Borel measurable. If • is integrable, then ∫ R ( ) H -1 ( {⃗ ∈ ∶ (⃗ ) = } ) d = ∫ ( ) |∇ |. Lemma 4.9. Approximation of a set of finite perimeter [32, Lemma 1.15]. If a bounded open set ⊂ R has Lipschitz boundary, and ⊂ R has finite perimeter, then there exists a sequence of open sets with being a nonempty compact hypersurface of differentiability class 2 and satisfying H -1 ( ∩ ) = 0, such that → 1 in 2 ( ), Per → Per , and | | = | |. Lemma 4.10. Diagonalization Argument. If double-indexed sequences { , } ⊂ R and { , } ⊂ R satisfy lim sup →∞ lim sup →∞ , ≤ , lim sup →∞ lim sup →∞ , ≤ , for some constants , ∈ R. Then there exists a diagonal sequence → ∞ as → ∞ such that lim sup →∞ , ≤ , lim sup →∞ , ≤ . Proof. Define * 1 = min { ∈ N ∶ lim sup →∞ , ≤ +1 and lim sup →∞ , ≤ +1 } . We claim that * 1 exists. Otherwise for each ∈ N, we have ∈ A ∪ B, where A ∶= { ∈ N ∶ lim sup →∞ , > + 1 } and B ∶= { ∈ N ∶ lim sup →∞ , > +1 } . Therefore | a contradiction. From the definition of * 1 , we know that * 1 (defined as follows) exists. * 1 ∶= min { ∈ N ∶ sup ≥ * 1 , ≤ +2 and sup ≥ * 1 , ≤ +2 } . Recursively, for = [2, ∞) ∩ N, define * = min { > * -1 ∶ lim sup →∞ , ≤ +1∕ and lim sup →∞ , ≤ +1∕ } , Q. Du et al. and * = min { > * -1 ∶ sup ≥ * , ≤ +2∕ and sup ≥ * , ≤ +2∕ } . For each ≥ 1 and each ∈ [ * , * +1 ) ∩ N, define = * , then , ≤ sup ≥ * * , ≤ +2∕ , , ≤ sup ≥ * * , ≤ +2∕ . □ 4.2.4. Convergence of global minimizers Theorem 4.11. E 0 is the Gamma-limit of E as → 0. Proof. See Propositions 4.4 and 4.7 for the liminf and limsup inequalities. □ Proposition 4.12. E has a global minimizer.
Proof. We use the direct method in the calculus of variations. We can show that for any fixed , any minimizing sequence { ( , )
} of E are uniformly bounded in
1,2 ( ) × 2 ( ), and thus possesses a subsequence (not relabeled) satisfying ⇀ ∞ in 1,2 ( ) and ⇀ ∞ in 2 ( ) for some ∞ ∈ 1,2 ( ) and ∞ ∈ 2 ( ). By the Rellich-Kondrachov theorem, → ∞ both in 2 ( ) and pointwise (up to a further subsequence). Similar to Step 5 of the Proof of Proposition 4.3, we can show ( ∞ , ∞ ) ⩽ lim inf →∞ ( , ). Similar to Step 2 of Proposition 4.4, we can show N( ∞ , ∞ ) = lim →∞ N( , ). Finally, since ⇀ ∞ in 1,2 ( ) and → ∞ in 2 ( ), by the Cauchy-Schwarz inequality we know ‖∇ ∞ ‖ 2 ( ) ⩽ lim inf →∞ ‖∇ ‖ 2 ( ) . Therefore, we have proven E ( ∞ , ∞ ) ⩽ lim inf →∞ E ( , ), and therefore ( ∞ , ∞ ) is a global minimizer of E . □ Theorem 4.13. Let ( , ) ∈ 1,2 ( ) × 2 ( ) be a global minimizer of E with → 0, then any subsequence of { ( , ) } has a further subsequence (not relabeled) such that → 0 in 2 ( ) and ⇀ 0 in 2 ( ), where ( 0 , 0 ) ∈
Proof. According to [17, Bottom of Page 18], {E } is equi-coercive thanks to the compactness result (Proposition 4.3). Since E 0 is the Gamma-limit of E as → 0 (Theorem 4.11), according to [17, Theorem 1.4.5], E 0 has a global minimizer ( 0 , 0 ) such that E 0 ( 0 , 0 ) = lim inf →0 E ( , ). Thanks to the compactness result (Proposition 4.
We are now left with some justifications to make in order to establish the equivalence between the Gamma-limit E 0 and the original problem . Recall that given by ( 1) is a functional of two subsets of R , with being the fundamental solution. In the definition of E 0 ( , ), the function is allowed to take any value between 0 and 1. However, as we show in Section 4.3.1, the minimizer ( * , * ) of E 0 satisfies * ∈ {0, 1} a.e., which means that * can be regarded as an indicator function. Since E 0 is the Gamma-limit of E , the boundary conditions on should be inherited from E to E 0 , which can be either no boundary conditions (i.e., is the fundamental solution) or periodic boundary conditions. In Section 4.3.2, we show the equivalence between those two types of boundary conditions.
We can regard E 0 as a relaxed version of , since the second argument of E 0 is allowed to take intermediate values between 0 and 1. As we now explain, such relaxation does not affect the energy minimizers.
Proposition 4.14. For any fixed ∈ ( ; {0, 1}), if * is a global minimizer of N( , ) (given by (11)) among all the ∈ 2 ( ; [0, 1]) such that = 0 a.e. and ∫ = ∫ , then * must be an indicator function, i.e., * ∈ {0, 1} a.e.
Proof. Our proof is very similar to [18,Equation (3.11)] (see also [16,Lemma 4]). Given any fixed > 0, we want to prove According to [36,Theorem 9.9], we have * ∈ 2 , ( ) for all 1 < < ∞, and * = - * ∕ a.e. on . In particular, we have * < -∕ a.e. on , because * > > 0 a.e. on and = 0 a.e. on (due to * = 0 a.e. on ). We compute the following first-order variation:
Due to the minimality condition on
* , the above left-hand side is zero, and thus * is constant a.e. on . According to [36, Lemma 7.7], which asserts that the weak derivatives of a Sobolev function are zero a.e. on its level set, we have ∇ * = ⃗ 0 a.e. on , and consequently * = 0 a.e. on , which is a contradiction. □ Nonlinear Analysis 250 (2025) 113665 Fig. 6. Step 3 in the Proof of Proposition 4.15. 4.3.2. Justification for boundary conditions 4.3.2.1. Compactly supported minimizers.
We now establish that periodic boundary conditions and no boundary conditions on are equivalent for the purpose of energy minimization in the sharp interface limit, as long as is sufficiently large so that the minimizers are compactly supported in . The key is the screening property satisfied by the energy minimizers [13, Corollary 3.3], i.e., different connected components of a minimizer do not interact with each other. So far we have been ambiguously using for the Green's functions ofunder either periodic boundary conditions or no boundary conditions. From now on, we let ̃ denote the former, and let denote the latter. We use analogous notations for other symbols as well: e.g., is given by ( 1) with the nonlocal kernel being the fundamental solution, and ̃ is the periodic counterpart with the nonlocal kernel being the Green's function under periodic boundary conditions.
Step 1: We first prove ( ,
, where ̃ is the periodic counterpart of (i.e., ̃ is the nonlocal term with the nonlocal kernel being the Green's function under periodic boundary conditions). According to Proposition 1.1-3 ⃝, we have , = 0 in ∖( ∪ ), where , is the electrostatic potential associated with and , given by (2). Since , satisfies periodic boundary conditions and the Poisson's equation -, = 1 -1 ∕ , we have ̃ , = constant + , , where ̃ , is the solution to -̃ , = 1 -1 ∕ under periodic boundary conditions. According to (3), we have 2 ( ,
In Step 3 we will prove ̃ ̃ , ̃ = constant in ∖( ̃ ∪ ̃ ), so we can extend ̃ ̃ , ̃ to be this constant in R ∖ . Since ̃ , ̃ vanishes at infinity [18, Lemma 3.1], ̃ ̃ , ̃ -̃ , ̃ is harmonic and bounded in R , therefore ̃ ̃ , ̃ -̃ , ̃ = constant. Similar to Step 1, we can prove ̃ ( ̃ , ̃ ) = ( ̃ , ̃ ).
Step 3: Now let us prove ̃ ̃ , ̃ = constant in ∖( ̃ ∪ ̃ ). Our proof is similar to [18,Section 4]. According to [13, Top of Page 8], is a volume perturbation of the perimeter, and is a quasi-minimizer of the perimeter, so is of differentiability class ∞ (Proposition 1.1-2 ⃝). Similarly, ̃ is also a volume perturbation of the perimeter, so we can assume that ̃ and ̃ are open, and that ̃ is of differentiability class ∞ . If ̃ is not open, we can instead consider the points at which the Lebesgue density of ̃ is positive, similar to [18, Page 1141]. See Fig. 6 for visualization.
Using the optimality of ̃ for ̃ ( ̃ , ⋅ ), we can prove the following by contradiction:
Otherwise, we can decrease the energy by moving some negative charge from a neighborhood of ⃗ to a neighborhood of ⃗ . More precisely, since ̃ ̃ , ̃ is continuous and
) where is small enough, and ̃ ( ̃ , ̃ ) < ̃ ( ̃ , ̃ ).
On ∖( ̃ ∪ ̃ ), since ̃ ̃ , ̃ is harmonic and satisfies the periodic boundary conditions, its minimum and maximum are attained on the boundary
Without loss of generality, let us assume that ̃ is connected. Because ̃ ̃ , ̃ ⩾ on both ̃ ∖ ̃ and ̃ ∩ ̃ , and -̃ ̃ , ̃ = 1 on ̃ , by the strong minimum principle for superharmonic functions, we have ̃ ̃ , ̃ > on ̃ . Since ̃ satisfies the interior ball condition ( ̃ is of differentiability class ∞ ), the Hopf boundary point lemma [36,Lemma 3.4] states that the outer normal derivative of ̃ ̃ , ̃ at ⃗ is negative, which is a contradiction to the minimality of , because there is ⃗ outside of ̃ ∪ ̃ such that ⃗ is very close to ⃗ and ̃ ̃ , ̃ ( ⃗ ) < . □ 4.3.2.2. Non-compactly supported minimizers. Even if is not large enough and the global minimizer of ̃ is not compactly supported in , we believe the asymptotics in Conjecture 4.16 to hold true. Our rationale is that according to [16,Theorem 8], any lowerdimensional structure with zero dipole moment can be extended to a higher dimension using a radially symmetric cutoff function, with asymptotically the same energy-to-mass ratio. , but we no longer require ( ̃ , ̃ ) to be a global minimizer or compactly supported in . Additionally, we assume that ( ̃ , ̃ ) has zero dipole moment, i.e., ∫ ⃗
(
which can be achieved by selecting a suitable translational representative according to Lemma 4.19. Let ( ̃ , ̃ ) denote the juxtaposition of copies of ( ̃ , ̃ ) (there are cycles along the direction of each standard basis vector in R , similar to the cubic crystal structure). Then
Remark 4.17.
⃝ For = 3, Conjecture 4.16 has a well-known discrete variant in which ̃ and ̃ are replaced by Dirac delta functions representing point charges with zero surface area. When calculating the electrostatic potential energy-to-mass ratio of a crystal made up of many (but finite) unit cells, solid physicists and material scientists usually take a shortcut: they consider only one unit cell and use Ewald summation which implicitly assumes periodic boundary conditions on Poisson's equation [38, Left of Page 7888]. Such Ewald summation gives an asymptotically correct value for a large but finite crystal surrounded by vacuum as long as the dipole moment in the unit cell is zero [39, Equation (1.8)], in which case the net Coulomb interaction decays sufficiently fast so that the summation is absolutely convergent. In this sense, periodic boundary conditions and no boundary conditions are equivalent for zero dipole moment.
⃝ For a nonzero dipole moment, the correct electrostatic energy-to-mass ratio is the sum of an intrinsic part and an extrinsic part.
The intrinsic part is given by the above Ewald summation and is dependent on the unit cell but independent of the global shape of the crystal [40]. The extrinsic part depends on the shape of the crystal, and can be interpreted as the outcome of ''effective'' charges distributed on the surface of the crystal (those surface charges reproduce the total dipole moment produced by the unit cells in the bulk) [39,41]. In this sense, periodic boundary conditions and no boundary conditions are not equivalent for a nonzero dipole moment.
⃝ The above extrinsic part should produce a voltage across the crystal formed by polar unit cells. However, this voltage seems to be absent in everyday life except in pyroelectric materials (whose polarization depends on the temperature). After a sufficiently long time at a constant temperature, even pyroelectric materials will lose such a voltage, because external charges will build up on the surface of the crystal through leakage currents (conducted by the crystal itself or the ambient atmosphere), thus canceling out that voltage in a similar way to grounding or earthing. Therefore, it seems more realistic to assume that the electrostatic potential is zero on the surface of the crystal, which leads to the homogeneous Dirichlet boundary conditions. According to Lemma 4.18, the electrostatic potential energy is actually minimized by the homogeneous Dirichlet boundary conditions. Many physicists and chemists refer to homogeneous Dirichlet boundary conditions as ''tin foil'' or ''conducting'' boundary conditions (imagine a crystal wrapped in a tin foil, or submerged in a conducting medium). They also believe that periodic boundary conditions and homogeneous Dirichlet boundary conditions are equivalent (see [42, The following lemma shows that within a bounded and connected domain, among all the electrostatic potentials that solve the same Poisson's equation, the electrostatic potential energy is minimized by the one satisfying homogeneous Dirichlet boundary conditions. Lemma 4.18. Assuming that is connected. For any fixed ∈ 1 , 2 ( ), consider the functional ℎ ↦ ‖∇( -ℎ)‖ 2 ( ) defined on the set consisting of all harmonic functions on . Its global minimizer ℎ * is unique (up to addition of a constant) and satisfies Dirichlet boundary conditions ℎ * = on .
Proof. Let ℎ * be a global minimizer, we can consider a competitor ℎ * + ℎ, where ∈ R and ℎ is any harmonic function on . Since ℎ * is stationary, we have
On the one hand, if ℎ * -is constant on , then the above right-hand side vanishes, because ∫ ℎ∕ ⃗ = ∫ ℎ = 0. On the other hand, if ℎ * -is non-constant on , then we can find a solution to Laplace's equation ℎ = 0 under Neumann boundary conditions ℎ∕ ⃗ = constant + -ℎ * on , so that ℎ * is not stationary against the small perturbations along the direction ℎ. In summary, ℎ * -must be constant on . □
The following lemma shows that any ̃ and ̃ contained in can be periodically (or circularly) shifted so that the dipole moment vanishes, thus satisfying the assumption in Conjecture 4.16. Without loss of generality, we only present the 2-D case ( = 2). Lemma 4.19. For any fixed
Nonlinear Analysis 250 (2025) 113665
Therefore, there exists * 1 ∈ [0, 1] such that for any 2 ∈ R, we have ℎ 1 ( * 1 , 2 ) = ℎ 1 ( * 1 , 0) = 0. Similarly, there exists * 2 ∈ [0, 1] such that for any 1 ∈ R, we have ℎ 2 ( 1 , * 2 ) = 0. In summary we have ⃗ ℎ( * 1 , * 2 ) = ⃗ 0. □
In this section, we present some numerical simulations based on the phase-field reformulation proposed in Section 4.
The numerical methods used in the simulations are very similar to those adopted in a previous work for the liquid drop model [19,Section III]. The main difference is that two phase-field functions are used here, whereas only one was used in the previous work.
Similar to [19,Section III.A], we use the following nonlinear function for the definition of E in Section 4.1:
For the time-marching scheme, we still use the convex splitting scheme, but our choice of the convex splitting is different from [19,Equation (8)]. For the potential well defined by (10), our choice here is = 1 + 2 , where 1 ( , ) = 87 2 ∕2 + 27 + 27 2 and 2 = -1 . It is easy to check the convexity of 1 and the concavity of 2 on [-0.1, 1.1] × [-0.1, 1.1], which are sufficient to stabilize our numerical scheme. Similar to [19,Equation (5)], in order to find the local minimizers, we use the following 2 gradient flow which is called the penalized Allen-Cahn-Ohta-Kawasaki (pACOK) dynamics:
where 1 , 2 > 0 are called mobility coefficients. We use the semi-implicit time-marching scheme similar to [19,Equation (8)].
In most simulations, the phase-fields almost vanish near . According to Section 4.3.2, no boundary conditions (i.e., is the fundamental solution) and periodic boundary conditions on are equivalent for the purpose of energy minimization. Therefore we only need to consider periodic boundary conditions, which allow us to use the Fourier spectral method similar to [19, Section III.D].
Although the perimeter of is not penalized in (1), we may add a penalty term ∫ |∇ | 2 to E in order for to converge faster. In practice, the penalty coefficient of ∫ |∇ | 2 can be chosen to be around 1∕1250000 of that of ∫ |∇ | 2 . Note that doing so is optional, and it will not significantly alter the numerical results but only provide some regularization for .
In Sections 5.1 and 5.2, we present some possible local minimizers as well as some snapshots taken from pACOK dynamics, in order to provide numerical evidence for our conjectures about the Gamma-expansion of the sharp interface energy in 2-D and 3-D for < 1 (see the first statement in Conjecture 3.7 and the first statement in Conjecture 3.8, respectively).
In our 2-D simulations, we choose = 1, = 1500, = 5 × 10 -2 , 1 = 3 × 10 4 , 2 = 24 × 10 3 ∕5, 1 = 1, 2 = 5, = 25 × 10 -5 ∕2 (which is the time step), and the simulation domain is chosen to be [0, 13∕5] × [0, 13∕5] which is then discretized into 256 × 256 uniform grid points. We visualize the results using RGB images of 256 × 256 pixels, where the hydrophobic region is drawn in purplish pink, and the hydrophilic region is drawn in greenish yellow. More specifically, we assign an RGB triplet to every possible value of ( , ). We assign (211, 95, 183), (220, 220, 98) and (255, 255, 255) to (1, 0), (0, 1) and (0, 0), respectively. For other values of ( , ), the RGB triplet is linearly interpolated from the above three, and then truncated to [0, 255] 3 . The color at the position ⃗ is then determined by the RGB triplet assigned to ( (⃗ ), (⃗ ) ) .
We present 13 local minimizers in Fig. 7. Next to each local minimizer is ( , E∕ ). Each local minimizer is obtained by choosing a suitable initial value, and then numerically simulating the pACOK dynamics until the iteration converges. The reason we treat them as local minimizers is that they appear to be stable against random perturbations in our simulations. One could further investigate their stability by evaluating the second-order variation of the energy or using rigorously validated numerics (similar to [44]).
As we can see in Fig. 7, the local minimizer is a micelle for = 0.1, which is consistent with the 2-D generalization of Proposition 1. 1234⃝. As increases to 0.15 and 0.2, the micelle is no longer stable and deforms into a shape that is similar to the eye-mask shaped local minimizers in [19,Figure 8]. As further increases to 0.4 and beyond, as shown in the second row of Fig. 7, the local minimizer becomes more and more elongated, resembling a straight bilayer of approximately uniform thickness, except that near the two ends the layer is slightly thicker while the layer is slightly thinner. The two open ends cause the straight bilayer to have slightly higher energy than the local minimizer shown in the first row of Fig. 7, which resembles a liposome. The liposome and straight bilayer seem to prefer roughly the same thickness, which is consistent with Remark 2.4-2
⃝. Their energy-to-mass ratios also seem to converge to roughly the same constant as → ∞, and the convergence rate for the liposome seems to be second-order, which is consistent with Corollary 2.2. The convergence rate for the straight bilayer seems to be first-order, indicating that the two open ends carry asymptotically constant energy penalties. Our numerical results therefore suggest a Gamma-expansion very similar to [14,Equation (1.6)] for = 2 and = 1.
For the liposome local minimizers in Fig. 7, we note that the inner layer is slightly thicker than the outer layer, and that such a difference becomes less noticeable as increases, which are consistent with Remark 3.6-1 ⃝. Our numerical evidence indicates that for = 2 and = 1, the optimal liposome candidate (whose asymptotics is given in Corollaries 2.2 and 3.4) is indeed a local minimizer, and might even be a global minimizer. In fact, if the first statement in Conjecture 3.7 is true, i.e., for < 1 the secondorder term in the Gamma-expansion of our energy is the elastica functional, then the global minimizer should be approximately (or exactly) circular according to Proposition 3.9.
In order to provide some numerical evidence for the first statement in Conjecture 3.7 and verify the resistance of the bilayer to bending, we carry out two simulations of the pACOK dynamics, as shown in Fig. 8. In the first simulation, we choose = 1.163, and we choose the initial value to be a bilayer of approximately uniform thickness, which resembles a random non-convex closed curve. We observe that over time, the bilayer becomes convex and resembles an ellipse, and eventually becomes circular after a sufficiently long time. The terminal value is shown in the top-right of Fig. 8, and its energy-to-mass ratio is indicated by the black cross × in the bottom-right of Fig. 7. In the second simulation, we choose = 0.6, and we choose the initial value by making a small hole in the liposome local minimizer shown in the top-middle of Fig. 7. We observe that over time, the bilayer straightens and eventually converges to a straight bilayer, which is a rigid transformation of the straight bilayer shown in the middle of the second row in Fig. 7. The convergence is relatively slow in both simulations, which is not surprising if the elastica functional is the second-order term in the Gamma-expansion. Fig. 9. Top: stationary points obtained in numerical simulations. Each rectangle represents the cross-section of the simulation box. Next to each stationary point are the mass and the energy-to-mass ratio. Bottom: curve fitting to the energy-to-mass ratios of the above stationary points.
In our 3-D simulations, we choose = 500, 1 = 1, 2 = 4, and the simulation domain [0, ] × [0, ] × [0, ] is discretized into N × M × P uniform grid points. We visualize the results using the following two MATLAB commands:
isosurface(u+v,1/2); isosurface(u,1/2);
where the former is set to be greenish yellow and transparent, representing the boundary of ∪ , while the latter is set to be purplish pink and opaque, representing the boundary of . In order to visualize the inner structures, we also plot the cross-section with the cutting plane parallel to the front view and passing through the center of the simulation box. We present 8 stationary points in Fig. 9, each of which is the terminal value of the pACOK dynamics starting from a suitable initial value and evolving over a long period of time until the shape barely changes. We present two simulations of the pACOK dynamics in Fig. 10. In Figs. 9 and 10, we choose = 1, = 6 × 10 -2 , 1 = 5∕2 × 10 4 , 2 = 4 × 10 3 , and = 2.1 × 10 -5 . In the top rows of Figs. 9 and 10, we choose = = = 3.66 and M = N = P = 256. In the last row in Fig. 10 and in the left three columns of the last row in Fig. 9, we choose = = 4 = 3.66 and M = N = 4P = 256. In the rightmost column of the last row in Fig. 9, we choose ∕2 = ∕2 = 4 = 3.66 and M∕2 = N∕2 = 4P = 256.
The stationary points in Fig. 9 resemble liposomes and disk bilayers, and they prefer roughly the same thickness, which is consistent with Remark 2.4-2 ⃝. Near the rim of a disk bilayer, the layer is slightly thicker while the layer is slightly thinner, which is a typical manifestation of frustration. Therefore, as → ∞, the disk bilayer should have a radius of order √ and a thickness of order 1. Its rim should have a perimeter of order √ and thus carry an energy penalty of order √ . Consequently, its energy-to-mass ratio should converge to a constant with order 1∕2, which is confirmed by the bottom of Fig. 9. We can also see that the energy-to-mass ratio of the liposome converges to roughly the same constant with order 1, which is consistent with Corollary 2.2. Our numerical calculations show that the liposome has lower energy than the disk bilayer, so that the latter cannot be a global minimizer. The liposome seems to be a local minimizer in our simulation. However, we are not confident that the disk bilayer is also a local minimizer, although its shape remains almost unchanged after 3.5 × 10 7 iterations. This is because the disk bilayer resembles an open surface, and we expect that the perimeter of its rim is penalized on the first order in the Gamma-expansion, and that the bending energy is only a second-order effect (see Section 3). Therefore, although the disk bilayer has zero bending energy, its priority should be to close in on itself and form a closed surface for sufficiently large . In fact, we can see that for = 7 the disk bilayer in Fig. 9 has higher energy than the curved bilayer shown in the top-left of Fig. 10, because the latter has smaller rim perimeter despite larger bending energy. As mentioned above, for a bilayer resembling an open surface, we expect a first-order energy penalty associated with its rim. In order to gain more insights, we carry out two simulations in Fig. 10. In the first and second simulations, the initial value is chosen by making a hole in the liposome ( = 7) and the disk bilayer ( = 2.4) that are obtained in Fig. 9, respectively. We observe that the hole diminishes and vanishes over time (cf. Bottom of Fig. 8 where the hole enlarges). Therefore, our numerical results demonstrate the self-healing property of lipid bilayers in 3-D, which is consistent with experimental observations [45] and is essential to biological membranes.
In order to provide some numerical evidence for the first statement in Conjecture 3.8 (i.e., the Willmore energy appears in the second-order Gamma-expansion), we present two numerical simulations in Fig. 11. Recall that the Clifford torus and its image under a conformal transformation are non-isolated local minimizers of the Willmore energy (see Fig. 14). In the first and second simulations, we start from an initial value resembling a torus and a deformed torus, respectively, let it evolve according to the pACOK dynamics over a long period of time until the shape barely changes, and plot the terminal value in the left and right of and its image under a conformal map, respectively. The former seems to be a local minimizer, while the latter has slightly higher energy than the former and seems to be evolving very slowly in the direction of becoming the former. We think the reason why the latter is not a local minimizer is that diffuse interfaces are used in our phase-field simulations, or that the bilayer has nonzero thickness, so that higher-order terms in the Gamma-expansion destroy the non-isolated local minimality. In order to provide numerical evidence for Remark 3.11-2 ⃝(i.e., triply periodic minimal surfaces may be preferred over planar bilayer for small ), we present four local minimizers shown in Fig. 12. They are obtained as the terminal values of the pACOK dynamics starting from suitable initial values after a sufficiently long time. We choose = 3.5×10 -2 , 1 = 30×10 4 ∕7, 2 = 48×10 3 ∕7, = 21 × 10 -5 ∕40, and M = N = P = 512. We choose = 0.6 for the left two local minimizers, and choose = 1 for the right two local minimizers. We choose = = = 3.51 for the left three local minimizers, and choose = = = 3.6855 for the rightmost local minimizer. From left to right, we choose = 7.6886, 11.5742, 7.1262 and 11.8021, respectively, so that the respective energy-tomass ratios are locally minimized with respect to . From Fig. 12 we can see that for = 0.6, the gyroid-like local minimizer has lower energy-to-mass ratio than the planar bilayer, and vice versa for = 1. This combined with Conjecture 4.16 is consistent with Remark 3.11-2 ⃝. In theory, the energy-to-mass ratio of the planar bilayer should equal that of the liposome in the limit of → ∞, which is given by 3 √ 9 ( +1)∕8 in Corollary 2.2. With = 500, this constant equals 9.6549 and 10.4004 for = 0.6 and 1, respectively, which are slightly different than the respective numerical energy-to-mass ratios of the first and third local minimizers shown in Fig. 12, with the relative error being -0.35% and -0.25%, respectively. We believe that this error is due to the diffuse interfaces used in our simulations. In fact, the initial values in this simulation are interpolated from the terminal values obtained in another simulation with a coarser grid, where we chose = 7 × 10 -2 and M = N = P = 256, with other parameters being the same. In the simulation with a coarser grid, we obtained four local minimizers similar to those shown in Fig. 12, with their energy-to-mass ratios being 9.5211, 9.4955, 10.2989, and 10.3184, respectively from left to right. As → 0 and M , N , P → ∞, we expect that the first and third energy-to-mass ratios converge to their respective theoretical values, and that the second and fourth ones maintain their relative differences to the first and third ones, respectively.
In Fig. 3, we can see that as increases, the optimal morphology should transition from bilayer membrane to cylindrical micelle to spherical micelle. In the liquid drop model [19], a ball loses stability when its mass exceeds a threshold. Similarly, we expect that there exist two thresholds of , beyond which the bilayer membrane and cylindrical micelle lose stability, respectively. For = 1, the straight bilayer membrane is stable on any 2-D periodic strip [21, Figure 1], and we think that the stability analysis therein can be generalized to any ∈ (0, ∞), thus allowing us to determine the threshold of beyond which the bilayer membrane is unstable. Note that such a threshold should be higher than 1 in Fig. 3. As we mentioned in Remark 3.6, the inner layer of the optimal liposome is slightly thicker but has slightly less mass compared to the outer layer. We can intuitively explain this phenomenon in Fig. 13-b, where a mismatch occurs as soon as a lipid bilayer membrane is curved. The inner monolayer becomes slightly more densely packed, while the outer monolayer becomes slightly less crowded. Those changes will inevitably increase the energy. To alleviate such a problem, some lipids may be transferred from the inner monolayer to the outer monolayer, as shown in Fig. 13-d. Our results indicate that when a closed bilayer membrane deforms (e.g., from a sphere to an ellipsoid), the lipids in the outer monolayer should flow from low curvature areas to high curvature areas, and vice versa for the inner monolayer. According to the fluid mosaic model [46], the lipids in a bilayer membrane can move easily within each monolayer, and they can also, albeit relatively slowly, move from one monolayer to the other (a movement known as flip-flop) [47,48]. The flip-flop process can be facilitated by certain proteins known as flippase, floppase and scramblase. The flippase moves lipids from the outer monolayer to the inner monolayer (flipping), the floppase does the opposite (flopping), and the scramblase does both.
In this paper, we have studied a degenerate version of the Ohta-Kawasaki energy and demonstrated its remarkable ability to reproduce the fascinating phenomena exhibited by self-assembling amphiphiles. We have presented some asymptotic and numerical evidence for the partial localization property of our model. Partial localization is coined in [14, Section 1.2] and refers to the concentration to lower-dimensional structures. An important example is the bilayer membrane formed by lipids: it is thin along one direction, but relatively large along the other two. Such a structure is vital to the biological membranes of every living cell. Our study may help us better understand the formation of lipid bilayer membranes. In fact, there are only two terms in our energy: the perimeter term which models the immiscibility between water and hydrophobic tails, and a Coulombic nonlocal term which models the attractive force between the heads and tails. It turns out that our model still possesses similar properties even if this Coulombic term is replaced by the 1-Wasserstein distance [14, Section 9.5], which penalizes the heads for straying far from the tails. Neither our Coulombic term nor the 1-Wasserstein distance keeps track of which head is connected to which tail. Therefore, the specific structure of the lipid molecule (a head and a tail connected by a covalent bond) is not the essence of partial localization, although it is a practical way to enforce the long-range attractive force that is needed in our model.
For the variant model mentioned above (in which the nonlocal term is the 1-Wasserstein distance), only the case of = 1 has been considered in [14,15]. In view of the rich complexity exhibited in our model for various , it might be interesting to revisit this variant model in the general cases > 0. On this note, we also draw attention to the well-known linkage between the 2-Wasserstein distance and the nonlocal (negative) Sobolev space norm [49].
Our study is a step towards understanding the pattern formation phenomena from the viewpoint of energetic competition. It is the competition between the short-and long-range terms in the Ohta-Kawasaki energy that gives rise to various interesting mesoscopic periodic patterns that are commonly observed in block copolymers and many other systems [10,[50][51][52]. We show that in the degenerate case (i.e., only the interface of is penalized), the Ohta-Kawasaki energy is capable of reproducing the partial localization feature of self-assembling amphiphiles. It is natural to ask how far can our results be generalized. For example, it might be of mathematical interest to explore other variant models with the Euclidean perimeter replaced by the 1-perimeter [53], a fractional perimeter [54], or a general nonlocal perimeter [55], with the Coulomb potential replaced by a Yukawa potential [56], a Riesz potential [57], a fractional inverse Laplacian kernel [58,Appendix], or a general nonlocal kernel [59].
= 0 -0 ∕ with = * and 0 = * 0 , then