DIFFUSED SOLUTE-SOLVENT INTERFACE WITH ... - Semantic Scholar

SIAM J. APPL. MATH. Vol. 75, No. 5, pp. 2072–2092

c 2015 Society for Industrial and Applied Mathematics 

DIFFUSED SOLUTE-SOLVENT INTERFACE WITH POISSON–BOLTZMANN ELECTROSTATICS: FREE-ENERGY VARIATION AND SHARP-INTERFACE LIMIT∗ BO LI† AND YUAN LIU‡ Abstract. A phase-field free-energy functional for the solvation of charged molecules (e.g., proteins) in aqueous solvent (i.e., water or salted water) is constructed. The functional consists of the solute volumetric and solute-solvent interfacial energies, the solute-solvent van der Waals interaction energy, and the continuum electrostatic free energy described by the Poisson–Boltzmann theory. All these are expressed in terms of phase fields that, for low free-energy conformations, are close to one value in the solute phase and another in the solvent phase. A key property of the model is that the phase-field interpolation of dielectric coefficient has the vanishing derivative at both the solute and the solvent phases. The first variation of such an effective free-energy functional is derived. Matched asymptotic analysis is carried out for the resulting relaxation dynamics of the diffused solute-solvent interface. It is shown that the sharp-interface limit is exactly the variational implicitsolvent model that has successfully captured capillary evaporation in hydrophobic confinement and corresponding multiple equilibrium states of underlying biomolecular systems as found in experiment and molecular dynamics simulations. Our phase-field approach and analysis can be used to possibly couple the description of interfacial fluctuations for efficient numerical computations of biomolecular interactions. Key words. variational implicit-solvent model, diffused solute-solvent interface, Poisson– Boltzmann theory, matched asymptotic analysis, sharp-interface limit AMS subject classifications. 35J, 35K, 49S, 92C DOI. 10.1137/15M100701X

1. Introduction. The solvation of charged molecules such as DNA and proteins in aqueous solvent (i.e., water or salted water) is a fundamental biological process. Implicit-solvent models provide efficient predictions of the solvation free energies and equilibrium conformations of underlying molecular systems [35]. In such a model, the solvent molecules and ions are treated implicitly and their effects are coarse-grained. In a large class of implicit-solvent models, the solute-solvent interface that separates the solvent region from solutes is used to calculate the surface energy and electrostatic energy. The latter is often done by solving Poisson’s or the Poisson–Boltzmann equation with the dielectric coefficient close to 1 and 80 in the solute and solvent regions, respectively [7, 16, 27, 38]. In a variational implicit-solvent model (VISM) [19, 20, 44], one determines the solvation free energies and stable equilibrium conformations by minimizing a macroscopic free-energy functional of all possible solute-solvent interfaces, i.e., dielectric boundaries. Let us denote by Ω ⊂ R3 the region of an underlying solvation system; cf. Figure 1. Assume that there are total N solute atoms located at x1 , . . . , xN inside Ω and ∗ Received by the editors February 4, 2015; accepted for publication (in revised form) July 24, 2015; published electronically September 15, 2015. http://www.siam.org/journals/siap/75-5/M100701.html † Department of Mathematics and Quantitative Biology Graduate Program, University of California, San Diego, La Jolla, CA 92093-0112 ([email protected]). This author’s work was supported by the U.S. National Institutes of Health (NIH) through grant R01GM096188 and the U.S. National Science Foundation (NSF) through grant DMS-1319731. ‡ Department of Mathematics, Fudan University, Yangpu District, Shanghai, 200433, People’s Republic of China, and Department of Mathematics, University of California, San Diego, La Jolla, CA 92093-0112 ([email protected]). This author’s work was supported by the China Scholar Council.

2072

DIFFUSED SOLUTE-SOLVENT INTERFACE

Ωp

2073

xi

Γ Ωw Ω Fig. 1. A schematic view of a solvation system with an implicit solvent. A solute-solvent interface Γ separates the solute region Ωp (p stands for protein) from the solvent region Ωw (w stands for water). Dots represent solute atoms located at xi carrying partial charges Qi (i = 1, . . . , N ).

carrying partial charges Q1 , . . . , QN , respectively. Assume also that there are M ionic species in the solvent, and denote by c∞ j and qj = zj e the bulk concentration and charge for the jth ionic species, respectively, where zj is the valence and e is elementary charge. Let Γ be a closed and smooth surface, a possible solute-solvent interface, that encloses all x1 , . . . , xN and that divides Ω into two regions: the solute region Ωp (p stands for protein) and the solvent region Ωw (w stands for water). The VISM solvation free-energy functional is then defined for all such solute-solvent interfaces Γ by  U (x) dx F [Γ] = P Vol(Ωp ) + γ0 Area(Γ) + ρw Ωw    εΓ + − |∇ψ|2 + ρf ψ − χw V (ψ) dx. (1.1) 2 Ω The first term here is the energy of creating a solute cavity in the solvent against the pressure difference P between the solvent liquid and solute vapor. The second term is the surface energy with γ0 being the constant surface tension. For simplicity, we neglect the Tolman correction to the surface tension [40]. The third term is the energy of the van der Waals (vdW)-type interaction between the solute atoms xi (1 ≤ i ≤ N ) and continuum solvent, where ρw is the constant bulk solvent density. The interaction potential U is often defined by U (x) =

N  i=1

where each (i) ULJ (r)

 = 4εi

(i)

ULJ (|x − xi |),

σi 12  σi 6 − r r



2074

BO LI AND YUAN LIU

is a Lennard-Jones potential, and all the parameters εi of energy and σi of length are given. The last term is the electrostatic free energy, where ψ is the electrostatic potential, ρf is a fixed charge density approximating the εΓ is the coefficient of

permittivity, N solute point charges i=1 Qi δxi with δxi being the Dirac delta function at xi , χw is the characteristic function of the solvent region Ωw , and V (ψ) describes the ionic contribution to the electrostatic interaction. Here, ions are assumed to be in an equilibrium macroscopically; and the equilibrium ionic concentrations are given by the Boltzmann distributions through the electrostatic potential ψ [7, 27]. The potential ψ is the solution to a boundary-value problem of the dielectric-boundary Poisson– Boltzmann equation (1.2) (1.3)

∇ · εΓ ∇ψ − χw V  (ψ) = −ρf on ∂Ω, ψ = ψ∞

in Ω,

where ψ∞ is a given function on the boundary ∂Ω. The coefficient εΓ is defined by  εp ε0 if x ∈ Ωp , (1.4) εΓ (x) = if x ∈ Ωw , εw ε0 where εp and εw are the dielectric coefficients (relative permittivities) of the solute and solvent regions, respectively, and ε0 is the vacuum permittivity. In general, εp ≈ 1 and εw ≈ 80. In the classical Poisson–Boltzmann theory, V (ψ) = β

−1

M 

−βqj ψ

e c∞ −1 , j

j=1

where β = 1/(kB T ) is the inverse thermal energy with kB the Boltzmann constant and T absolute temperature. Different forms of V (ψ) can be used for the linearized or size-modified Poisson–Boltzmann equation [26, 45]. More complicated forms V = V (ψ, ∇ψ) involving both ψ and ∇ψ can be used to model effects such as the ionic concentration dependent dielectric response but can also be more complicated for implementation [2, 23, 29]. Cheng and co-workers [9, 10, 12, 41, 44] developed a highly accurate, compact coupling interface method for solving the dielectric-boundary Poisson–Boltzmann equation and a robust level-set method to minimize the VISM functional (1.1). The normal velocity in the level-set formulation is the (normal component of) effective boundary force defined to be the negative variational derivative −δΓ F [Γ] : Γ → R of the freeenergy functional F [Γ] with respect to the location change of interface Γ. Let n denote the unit normal to the interface Γ pointing from the solvent region Ωw to solute region Ωp . We have [5, 9, 10, 11, 28, 44]   2 1 ∂ψ 1 1 −δΓ F [Γ] = P + 2γ0 H − ρw U + − εΓ 2 εp ε0 εw ε0 ∂n 1 2 (1.5) + (εw ε0 − εp ε0 ) |∇Γ ψ| + V (ψ) on Γ, 2 where H is the mean curvature that is positive if Ωp is a sphere and ∇Γ = (I −n⊗n)∇, with I the identity matrix, is the surface gradient along Γ. While successful initially [9, 12, 22, 41, 44], the sharp-interface VISM needs to be improved in several aspects. The most critical one is to include the description of

DIFFUSED SOLUTE-SOLVENT INTERFACE

2075

fluctuations, both around the solute-solvent interface and in the bulk solvent. Such fluctuations are particularly crucial in the transition of one equilibrium conformation to another and in sampling different states to accurately predict the free energies of underlying biomolecular systems. It is certainly possible to describe interfacial fluctuations within a sharp-interface framework. But several implementational difficulties can arise. For instance, the extension of normal velocity in the level-set method can be hard for a moving fluctuating interface. Moreover, fluctuations can nucleate and coalesce small bubbles (water regions) inside solutes, making it hard to solve the dielectric-boundary Poisson–Boltzmann equation in a sharp-interface formulation. In seeking an alternative approach, we notice that initial theoretical and computational studies of interfacial fluctuations using a diffused-interface approach seem promising [3, 17, 24]. We therefore propose in this work a diffused-interface approach to the solvation of charged molecules. In slightly different language, this is a phase-field approach, as it is well appreciated that the solvent-solute interface in a biomolecular system resembles a liquid-vapor interface, and the solvent and solute can be regarded as two different phases [6, 42]. The phase-field approach has been widely used in studying interface problems arising in many scientific areas, such as materials physics, complex fluids, biomembranes, and cell motility; cf., e.g., [1, 13, 18, 25, 31, 37] and the references therein. Our phase-field model is governed by the free-energy functional     ξ 1 2 2 φ dx + γ0 |∇φ| + W (φ) dx + ρw (φ − 1)2 U dx Fξ [φ] = P ξ Ω Ω 2 Ω   ε(φ) + |∇ψφ |2 + ρf ψφ − (φ − 1)2 V (ψφ ) dx (1.6) − 2 Ω for all phase fields φ : Ω → R, where the electrostatic potential ψφ is the solution to the boundary-value problem of the Poisson–Boltzmann equation with a diffused dielectric boundary (1.7)

∇ · ε(φ)∇ψφ − (φ − 1)2 V  (ψφ ) = −ρf

(1.8)

ψφ = ψ∞

in Ω,

on ∂Ω.

Here ξ > 0 is a small number. With a low free energy, a phase field φ should be close to two values, say, 0 and 1, respectively, in the solvation region Ω, except a thin transition layer that represents the dielectric boundary. The sets {φ ≈ 0} and {φ ≈ 1} represent the solvent and solute regions, respectively. With such a convention, the first term in (1.6) corresponds to that in (1.1), describing the volumetric energy. The second term in (1.6), in which W (φ) = 18φ2 (1 − φ)2 , approximates the surface energy [30, 33, 39]. The third term in (1.6) corresponds to that in (1.1), describing the energy of solute-solvent interaction. Finally, the last term in (1.6) corresponds to that in (1.1), describing the electrostatic free energy. In this term, ε(φ) is defined to be a smooth function of φ that interpolates the piecewise constant dielectric coefficient: ε(0) = εw ε0 and ε(1) = εp ε0 . We shall require that ε (0) = ε (1) = 0. Note that the electrostatic potential ψφ depends on φ. We study our proposed phase-field solvation free-energy functional (1.6) and the corresponding relaxation dynamics. Our main results are as follows:

2076

BO LI AND YUAN LIU

(1) We prove the well-posedness of the boundary-value problem of Poisson– Boltzmann equation (1.7) and (1.8) and obtain some bounds for the electrostatic potential; cf. Theorem 2.1. (2) We derive, first intuitively and then rigorously, the first variation of the phasefield free-energy functional (1.6):  1  δφ F [φ] = 2P φ − γ0 ξΔφ − W (φ) + 2ρw (φ − 1)U ξ 1 − ε (φ)|∇ψ|2 − 2(φ − 1)V (ψ); 2 cf. Theorem 3.1. (3) We consider the relaxation dynamics ∂t φ = −δφ F [φ] for φ = φ(x, t) with t denoting the time variable. Using the method of matched asymptotic analysis [4, 14, 15, 34, 36], we show that, in the sharp-interface limit, the normal velocity of the solute-solvent interface is exactly the boundary force (1.5) in the original sharp-interface variational approach; cf. Theorem 5.1. Several remarks are in order. First, our diffuse-interface model (1.6) differs from our previous one [30, 43] significantly in that we use the Poisson–Boltzmann equation rather than the Coulomb-field approximation [8, 41], which requires no solution of partial differential equations, to describe the electrostatic interaction. Second, our relaxation dynamics ∂t φ = −δφ F [φ] is a nonconservative dynamics with respect to the phase-field function φ. In fact, the relaxation process can change the volume of solute region {φ ≈ 1}. The amount of ions, however, is controlled through their bulk concentrations c∞ j (j = 1, . . . , M ) that are input parameters. Third, the surface energy term described by the gradient-square and double-well terms in our phasefield energy functional (1.6) coincides with the part in an effective Hamiltonian for the large-scale solvent density in the Lum–Chandler–Weeks theory of solvation [32]. The rest of the paper is organized as follows: In section 2, we prove the existence and uniqueness of, and provide bounds for, the solution to the Poisson–Boltzmann equation with a diffused dielectric boundary that is given by a phase field. In section 3, we derive the first variation of the phase-field free-energy functional. In section 4, we use the matched asymptotic analysis to study the relaxation dynamics in fast and regular time scales. In section 5, we continue our matched asymptotic analysis for a slow time scale to show that the sharp-interface limit of the relaxation dynamics is exactly the same as in the original sharp-interface variational model. 2. Poisson–Boltzmann electrostatics with a diffused dielectric boundary. We make the following assumptions: (A1) The solvation region Ω is a bounded, connected, open set with a sufficiently smooth boundary ∂Ω. All x1 , . . . , xN are given points in Ω. All P , γ0 , ρw , εp , and εw are positive constants with εp < εw . The function ρf ∈ L∞ (Ω). The function ψ∞ ∈ W 1,∞ (Ω). (A2) The function U : R3 \ {x1 , . . . , xN } → R is continuously differentiable and bounded below. Moreover, limx→xi U (x) = +∞ for each i (1 ≤ i ≤ N ) and limx→∞ U (x) = 0. (A3) The function ε ∈ C 1 (R). There exist positive numbers εmin , εmax , and εmax such that εmin ≤ ε(φ) ≤ εmax and |ε (φ)| ≤ εmax for all φ ∈ R. Moreover, ε(0) = εw ε0 , ε(1) = εp ε0 , ε (0) = ε (1) = 0, and ε (φ) = 0 if 0 < φ < 1. (A4) The function V ∈ C 2 (R). It is strictly convex. Moreover, V (0) = mins∈R V (s) < V (s) for any s = 0, V (±∞) = +∞, and V  (±∞) = ±∞.

DIFFUSED SOLUTE-SOLVENT INTERFACE

2077

Here are two examples of the function ε = ε(φ) with all the desired properties. Example 1. ⎧ ε ε if φ ≤ 0, ⎪ ⎪ ⎨ w 0 tan(π(φ−1/2)) + εw ε0 e− tan(π(φ−1/2)) εp ε0 e ε(φ) = if 0 < φ < 1, ⎪ etan(π(φ−1/2)) + e− tan(π(φ−1/2)) ⎪ ⎩ if φ ≥ 1. εp ε0 In fact, this is a C ∞ -function. Moreover, for any integer m ≥ 1, ε(m) (φ) = 0 for all φ ∈ (0, 1) and ε(m) (φ) = 0 for all φ ∈ (0, 1). Example 2. ⎧ if φ ≤ 0, ⎪ ⎨ εw ε0 1 1 ε(φ) = (εw ε0 − εp ε0 ) cos πφ + (εw ε0 + εp ε0 ) if 0 < φ < 1, ⎪ 2 ⎩ 2 εp ε0 if φ ≥ 1. In what follows, we shall denote by C a generic constant that may depend on Ω, ρf , ψ∞ , and ε but is independent of φ ∈ H01 (Ω). We shall also denote by C(φ) a generic constant that can depend on φ, in addition to Ω, ρf , ψ∞ , and ε. The following theorem is our main result in this section. Theorem 2.1. For any φ ∈ H01 (Ω) there exists a unique weak solution ψφ ∈ 1 H (Ω) to the boundary-value problem (1.7) and (1.8); i.e., ψφ = ψ∞ on ∂Ω in the sense of trace and     2  ε(φ)∇ψφ · ∇η + (φ − 1) V (ψφ ) η dx = ρf η dx ∀η ∈ H01 (Ω). (2.1) Ω

Ω

Moreover, ψφ ∈ L∞ (Ω) and

ψφ H 1 (Ω) ≤ C 1 + φ L2 (Ω) ,

χ{φ=1} ψφ L∞ (Ω) ≤ C,  

ψφ L∞ (Ω) ≤ C 1 + φ 2H 1 (Ω) .

(2.2) (2.3) (2.4)

Proof. We divide our proof in three steps. Step 1. Let ψˆφ ∈ H 1 (Ω) be the unique weak solution to the boundary-value problem −∇ · ε(φ)∇ψˆφ = ρf in Ω and ψˆφ = ψ∞ on ∂Ω. This means that ψˆφ is the unique function in H 1 (Ω) such that ψˆφ = ψ∞ on ∂Ω and   (2.5) ε(φ)∇ψˆφ · ∇η dx = ρf η dx ∀η ∈ H01 (Ω). Ω

Ω

Since 0 < εmin ≤ ε(φ) ≤ εmax for all φ ∈ R, we have by the standard elliptic theory (cf. Theorems 8.3 and 8.16 in [21]) that ψˆφ ∈ H 1 (Ω) ∩ L∞ (Ω), and (2.6)

ψˆφ H 1 (Ω) + ψˆφ L∞ (Ω) ≤ C.

Step 2. We define I : H01 (Ω) → R ∪ {+∞} by   ε(φ) 2 2 ˆ I[u] = |∇u| + (φ − 1) V (u + ψφ ) dx 2 Ω

∀u ∈ H01 (Ω).

2078

BO LI AND YUAN LIU

It is easy to see that inf u∈H01 (Ω) I[u] is finite. By the direct method in the calculus of variations, there exists a unique minimizer u ∈ H01 (Ω) of I over H01 (Ω). By our assumption (A4) on the function V and by (2.6), there exists λ > 0, independent of φ, such that V  (λ + ψˆφ ) > 1 and V  (−λ + ψˆφ ) < −1 a.e. on Ω. We prove that |u| ≤ λ

(2.7) Define uλ : Ω → R by

⎧ ⎨ −λ u(x) uλ (x) = ⎩ λ

a.e. on {φ = 1}.

if u(x) < −λ, if − λ ≤ u(x) ≤ λ, if u(x) > λ.

Since I[uλ ] ≥ I[u] and |∇uλ | ≤ |∇u| a.e. on Ω, we have   ε(φ) 2 ˆ |∇u|2 dx (φ − 1) V (u + ψφ ) dx = I[u] − Ω Ω 2  ε(φ) |∇uλ |2 dx ≤ I[uλ ] − Ω 2  = (φ − 1)2 V (uλ + ψˆφ ) dx. Ω

Consequently, it follows from the convexity of V and our choice of λ that   0 ≥ (φ − 1)2 V (u + ψˆφ ) dx − (φ − 1)2 V (uλ + ψˆφ ) dx Ω Ω   = χ{u>λ} (φ − 1)2 V (u + ψˆφ ) − V (λ + ψˆφ ) dx Ω    + χ{uλ} (φ − 1)2 V  (λ + ψˆφ )(u − λ)dx Ω  + χ{uλ} (φ − 1)2 (u − λ)dx + χ{uλ} (φ − 1) | |u| − λ | dx Ω

≥ 0. This leads to (2.7). The minimizer u of I : H01 (Ω) → R ∪ {+∞} is the weak solution to the Euler– Lagrange equation ∇ · ε(φ)∇u − (φ − 1)2 V  (u + ψˆφ ) = 0

in Ω,

i.e., (2.8)

   ε(φ)∇u · ∇η + (φ − 1)2 V  (u + ψˆφ ) η dx = 0 Ω

∀η ∈ H01 (Ω).

DIFFUSED SOLUTE-SOLVENT INTERFACE

2079

This is true by (2.7) and the Lebesgue dominated convergence theorem if η ∈ H01 (Ω)∩ L∞ (Ω). For a general η ∈ H01 (Ω), this is also true, since (φ − 1)2 V  (u + ψˆφ ) ∈ L4 (Ω) by (2.7) and since H01 (Ω) ∩ L∞ (Ω) is dense in H01 (Ω). The convexity of V implies that the weak solution u ∈ H01 (Ω) defined by (2.8) is unique. If we choose η = u in the equation of (2.8), we obtain by (2.6), (2.7), and Poincar´e’s inequality that      

u 2H 1 (Ω) ≤ C (φ − 1)2 V  (u + ψˆφ ) u dx ≤ C (1 + φ2 ) dx, {φ=1}

leading to (2.9)

Ω



u H 1 (Ω) ≤ C 1 + φ L2 (Ω) .

It now follows from the regularity theory (cf. Theorem 8.16 in [21]), (2.6), (2.7), and the imbedding H 1 (Ω) → L4 (Ω) that (2.10)

   

u L∞ (Ω) ≤ C (φ − 1)2 V  (u + ψˆφ ) L2 (Ω) ≤ C 1 + φ 2L4 (Ω) ≤ C 1 + φ 2H 1 (Ω) .

Step 3. Let ψφ = u + ψˆφ . Then ψφ = ψ∞ on ∂Ω. Moreover, (2.5) and (2.8) imply (2.1). The uniqueness of ψφ follows from the convexity of V and a usual argument. The estimates (2.2)–(2.4) follow from (2.6), (2.9), (2.7), and (2.10). . 3. First variation of free-energy functional. We calculate the first variation of the free-energy functional Fξ defined by (1.6). The first variation of each of the first three parts in the functional Fξ can be obtained by routine calculations. So, we focus on the electrostatic part   ε(φ) (3.1) Fele,ξ [φ] := |∇ψφ |2 + ρf ψφ − (φ − 1)2 V (ψφ ) dx, − 2 Ω where ψφ is the unique solution to the boundary-value problem (1.7) and (1.8). Heuristically, if we denote by δφ the “derivative” with respect to φ and by δφ the variation of φ, and “differentiate” both sides of the above equation, then we obtain by the chain rule and product rule for differentiation that    ε (φ) δφ|∇ψφ |2 − ε(φ)∇ψφ · ∇δφ ψφ + ρf δφ ψφ δφ Fele,ξ [φ]δφ = − 2 Ω − 2(φ − 1)δφV (ψφ ) − (φ − 1)2 V  (ψφ )δφ ψφ dx    ε (φ) |∇ψφ |2 − 2(φ − 1)V (ψφ ) δφ dx − = 2 Ω     + ε(φ)∇ψφ · ∇δφ ψφ + (φ − 1)2 V  (ψφ )δφ ψφ dx ρf δφ ψφ dx − Ω   Ω ε (φ) = |∇ψφ |2 − 2(φ − 1)V (ψφ ) δφ dx, − 2 Ω where in the last step we used the weak form (2.1) of the Poisson–Boltzmann equation with η = δφ ψφ . We then identify δφ Fele,ξ [φ] as δφ Fele,ξ [φ] = −

ε (φ) |∇ψφ |2 − 2(φ − 1)V (ψφ ). 2

2080

BO LI AND YUAN LIU

We remark that −δφ Fele,ξ [φ] is an approximation of the electrostatic part of the boundary force −δΓ F [Γ] in (1.5) if φ ≈ 0 in the solvent region and φ ≈ 1 in the solute region. In fact, the V -terms are similar, since the perturbation δφ is localized at the interface. With our notation, the unit normal n ≈ −∇φ/|∇φ|. Moreover, ε(φ) ≈ εΓ . Therefore, ε (φ) ε (φ) ε (φ) |∇ψφ |2 = |∇ψφ · n|2 + |∇Γ ψφ |2 2 2 2  1 1 d ε (φ) 2 |∇Γ ψφ |2 =− |ε(φ)∇ψφ · n| + 2 dφ ε(φ) 2   1 1 1 1 2 ≈− − |εΓ ∇ψφ · n| + (εp ε0 − εw ε0 )|∇Γ ψφ |2 . 2 εp ε0 εw ε0 2 The final result is exactly the corresponding part in (1.5). We now give a rigorous justification of our derivation. We recall that if G : H01 (Ω) → R is a functional and if φ, η ∈ H01 (Ω), then the first variation of G at φ in the “direction” η is defined to be  d  G[φ + tη] − G[φ] =  G[φ + tη] δφ G[φ][η] = lim t→0 t dt t=0 if the limit exists. Theorem 3.1. Let φ ∈ H01 (Ω). Let ψφ ∈ H 1 (Ω) be the weak solution of the corresponding boundary-value problem (1.7) and (1.8). Let η ∈ L∞ (Ω) ∩ H01 (Ω). Then   1 (3.2) δφ Fele,ξ [φ][η] = − ε (φ)|∇ψφ |2 − 2(φ − 1)V (ψφ ) η dx. 2 Ω To prove the theorem, we need the following lemma. Lemma 3.2. Let φ ∈ H01 (Ω) and η ∈ L∞ (Ω) ∩ H01 (Ω). Let t ∈ R with |t| ≤ 1. Let ψφ and ψφ+tη be the weak solutions to the boundary-value problems of Poisson– Boltzmann equation (1.7) and (1.8) corresponding to φ and φ + tη, respectively. Then there exists a constant C = C(φ, η) > 0 that may depend on both φ and η such that

ψφ+tη − ψφ H 1 (Ω) ≤ C(φ, η)|t|

if |t| ≤ 1.

Proof. By the weak form for ψφ and that for ψφ+tη (cf. (2.1)), we have for any ζ ∈ H01 (Ω) that     2  ε(φ)∇ψφ · ∇ζ + (φ − 1) V (ψφ ) ζ dx = ρf ζ dx, Ω Ω    ε(φ + tη)∇ψφ+tη · ∇ζ + (φ + tη − 1)2 V  (ψφ+tη ) ζ dx = ρf ζ dx. Ω

Ω

It follows from these two equations with ζ = ψφ+tη − ψφ , the property of ε, and the convexity of V that εmin ∇(ψφ+tη − ψφ ) 2L2 (Ω)  ≤ ε(φ)∇(ψφ+tη − ψφ ) · ∇(ψφ+tη − ψφ ) dx Ω   = ε(φ)∇ψφ+tη · ∇(ψφ+tη − ψφ ) dx − ε(φ)∇ψφ · ∇(ψφ+tη − ψφ ) dx Ω

Ω

DIFFUSED SOLUTE-SOLVENT INTERFACE

2081

 =−

Ω

+ Ω =−  − −

Ω

 ≤−  −

Ω

Ω

Ω

[(ε(φ + tη) − ε(φ)] ∇ψφ+tη · ∇(ψφ+tη − ψφ ) dx   (φ + tη − 1)2 − (φ − 1)2 V  (ψφ+tη ) (ψφ+tη − ψφ ) dx

(φ − 1)2 [V  (ψφ+tη ) − V  (ψφ )] (ψφ+tη − ψφ ) dx

Ω

Ω

[(ε(φ + tη) − ε(φ)] ∇ψφ+tη · ∇(ψφ+tη − ψφ ) dx  ε(φ + tη)∇ψφ+tη · ∇(ψφ+tη − ψφ ) dx − ε(φ)∇ψφ · ∇(ψφ+tη − ψφ ) dx

[(ε(φ + tη) − ε(φ)] ∇ψφ+tη · ∇(ψφ+tη − ψφ ) dx   (φ + tη − 1)2 − (φ − 1)2 V  (ψφ+tη ) (ψφ+tη − ψφ ) dx

≤ εmax |t| η L∞ (Ω) ∇ψφ+tη L2 (Ω) ∇(ψφ+tη − ψφ ) L2 (Ω)

+ 2|t| η(φ − 1) L2 (Ω) + t2 η 2 L2 (Ω) V  (ψφ+tη ) L∞ (Ω) ψφ+tη − ψφ L2 (Ω) . Since |t| ≤ 1, we have by (2.2) and (2.4) that both ψφ+tη H 1 (Ω) and V  (ψφ + tη) L∞ (Ω) are bounded by some constants that can depend on φ and η. An application of Poincar´e’s inequality then implies the desired inequality. Proof of Theorem 3.1. For any t with |t|  1, we denote by ψφ+tη the weak solution to the boundary-value problem (1.7) and (1.8) with φ+tη replacing φ. Setting ψφ+tη − ψφ ∈ H01 (Ω) as the test function in the corresponding weak formulation for ψφ+tη , we obtain   ρf (ψφ+tη − ψφ ) dx = ε(φ + tη)∇ψφ+tη · ∇(ψφ+tη − ψφ ) dx Ω Ω  + (φ + tη − 1)2 V  (ψφ+tη )(ψφ+tη − ψφ ) dx Ω

1 = ε(φ + tη) |∇ψφ+tη |2 + |∇ψφ+tη − ∇ψφ |2 − |∇ψφ |2 dx 2 Ω  + (φ + tη − 1)2 V  (ψφ+tη )(ψφ+tη − ψφ ) dx. Ω

Consequently, by the definition of the functional Fele,ξ (cf. (3.1)), we have Fele,ξ [φ + tη] − Fele,ξ [φ]   ε(φ) ε(φ + tη) |∇ψφ+tη |2 + |∇ψφ |2 + ρf (ψφ+tη − ψφ ) = − 2 2 Ω 2 2 − (φ + tη − 1) V (ψφ+tη ) + (φ − 1) V (ψφ ) dx  1 = ε(φ + tη)|∇ψφ+tη − ∇ψφ |2 dx Ω 2  1 [ε(φ + tη) − ε(φ)] |∇ψφ |2 dx − 2 Ω   (φ + tη − 1)2 − (φ − 1)2 V (ψφ ) dx − Ω

2082

BO LI AND YUAN LIU

 − (3.3)

Ω

(φ + tη − 1)2 [V (ψφ+tη ) − V (ψφ ) − V  (ψφ+tη )(ψφ+tη − ψφ )] dx

=: Λ1 (t) − Λ2 (t) − Λ3 (t) − Λ4 (t),

where Λi (t) (i = 1, . . . , 4) denote the corresponding integrals. By Lemma 3.2, we have    Λ1 (t)  εmax 2   (3.4)  t  ≤ 2|t| ψφ+tη − ψφ H 1 (Ω) ≤ C(φ, η)t → 0 as t → 0. Since η ∈ L∞ (Ω), we have for a.e. x ∈ Ω that 1 [ε(φ + tη) − ε(φ)] = ε (φ + θtη)η → ε (φ)η t

as t → 0,

where θ = θ(x) ∈ [0, 1]. The Lebesgue dominated convergence theorem then implies that  Λ2 (t) 1  (3.5) lim = ε (φ)|∇ψφ |2 η dx. t→0 t Ω 2 Similarly, since ψφ ∈ L∞ (Ω) by Theorem 2.1,  Λ3 (t) (3.6) lim = 2(φ − 1)V (ψφ )η dx. t→0 t Ω Let a ∈ R, and define g(b) = V (b) − V (a) − V  (b)(b − a) for all b ∈ R. Note that g(a) = 0 and g  (b) = −V  (b)(b − a). By Taylor’s expansion, g(b) = g(a) + g  (a + μ(b − a))(b − a) = −μV  (a + μ(b − a))(b − a)2 , where μ = μ(b) ∈ [0, 1]. With a = ψφ and b = ψφ+tη , we obtain by (2.4), the Cauchy– Schwarz inequality, the imbedding H 1 (Ω) → L4 (Ω), and Lemma 3.2 that     Λ4 (t)  C(φ, η)  ≤ (φ + tη − 1)2 |ψφ+tη − ψφ |2 dx  t  t Ω C(φ, η) ≤

(φ + tη − 1)2 L2 (Ω) ψφ+tη − ψφ 2L4 (Ω) t C(φ, η)

(φ + tη − 1)2 L2 (Ω) ψφ+tη − ψφ 2H 1 (Ω) ≤ t (3.7) →0 as t → 0. Now (3.2) follows from (3.3)–(3.7). The next theorem provides the formula of first variation of the solvation freeenergy functional Fξ defined in (1.6). It follows from Theorem 3.1 and routine calculations. We thus omit the proof. Theorem 3.3. Let φ ∈ H01 (Ω) and η ∈ L∞ (Ω) ∩ H01 (Ω) be such that   2 (φ − 1) U dx < ∞ and η 2 U dx < ∞. Ω

Ω

DIFFUSED SOLUTE-SOLVENT INTERFACE

2083

Let ψφ ∈ H 1 (Ω) be the unique weak solution to the boundary-value problem (1.7) and (1.8) corresponding to φ. Then    γ0 γ0 ξ∇φ · ∇η + 2P φ + W  (φ) + 2ρw (φ − 1)U δφ Fξ [φ][η] = ξ Ω  ε (φ) − |∇ψφ |2 − 2(φ − 1)V (ψφ ) η dx. 2 If, in addition, φ ∈ H 2 (Ω), then the integral of ∇φ · ∇η becomes that of −Δφ η. Therefore, we shall denote (3.8) δφ Fξ [φ] = 2P φ − γ0



ε (φ) 1  |∇ψφ |2 − 2(φ − 1)V (ψφ ) ξΔφ − W (φ) + 2ρw (φ − 1)U − ξ 2

and call it the first variation of Fξ at φ. 4. Relaxation dynamics: Fast and regular time scales. We now consider the relaxation dynamics ∂t φ = −δφ Fξ [φ]. By Theorem 3.3, this dynamics is described by the following initial-boundary-value problem for the phase field φ = φ(x, t) and the corresponding electrostatic potential ψ = ψ(x, t):  1  ∂t φ = −2P φ + γ0 ξΔφ − W (φ) − 2ρw (φ − 1)U ξ 1  2 + ε (φ)|∇ψ| + 2(φ − 1)V (ψ), (4.1) 2 (4.2) ∇ · ε(φ)∇ψ − (φ − 1)2 V  (ψ) = −ρf , (4.3) (4.4)

φ=0 on ∂Ω, ψ = ψ∞ on ∂Ω,

together with some initial condition for φ. If the solutions φ = φ(x, t) and ψ = ψ(x, t) are smooth, then we have by (1.6), (3.8), the chain rule, and integration by parts that  d Fξ [φ(·, t)] = ∂t φ δφ Fξ [φ] dx dt Ω    − ε(φ)∇ψ · ∇∂t ψ − ρf ∂t ψ + (φ − 1)2 V  (ψ)∂t ψ dx Ω = − (∂t φ)2 dx ≤ 0, Ω

where the vanish of the second integral follows from the weak formulation for (4.2) with the test function ∂t ψ. Therefore, the energy decays, explaining the terminology “relaxation dynamics.” Assume 0 < ξ  1. We expect that, after a short time period, the solution φ will take the value close to 0 or 1 in most of the region Ω, due to the fast reaction described by the term −(1/ξ)W  (φ). Moreover, in the regions {φ ≈ 0} and {φ ≈ 1}, the leading-order term of the electrostatic potential should solve the boundary-value problem (4.2) and (4.4), as the electrostatic relaxation is instantaneous: there is no ∂t ψ term. In the subsequent long time period, the relaxation dynamics is mainly the motion of the interface that separates the two regions {φ ≈ 0} and {φ ≈ 1}. It is reasonable to assume that the early fast time scale is determined by the basic reaction-diffusion equation 1 ∂t φ = ξΔφ − W  (φ), ξ

2084

BO LI AND YUAN LIU

as again the electrostatic relaxation is instantaneous. Perturbing the unstable constant equilibrium solution φ0 (x) = 1/2 by δ exp (ik · x + ωt) with |δ|  1, we find the most unstable mode kc = 0 and the corresponding fastest growth rate ω(kc ) = O(1/ξ). Therefore, the fast time scale is t/ξ. The next time scales are the regular time scale t, a slow time scale ξt, and so on. Let us first consider the fast time scale T = t/ξ. We assume that (4.5)

φ(x, t) = φ0 (x, T ) + ξφ1 (x, T ) + ξ 2 φ2 (x, T ) + · · · ,

(4.6)

ψ(x, t) = ψ0 (x, T ) + ξψ1 (x, T ) + ξ 2 ψ2 (x, T ) + · · · ,

where φi = φi (x, T ) and ψi = ψi (x, T ) (i = 0, 1, . . .) are smooth and bounded functions satisfying the boundary conditions (4.7)

φi = 0

(i ≥ 0),

ψ0 = ψ∞ ,

and ψi = 0 (i ≥ 1)

on ∂Ω.

Substitute φ and ψ in (4.1) by these expansions, and match terms with the same orders of ξ. At the leading order O(ξ −1 ), we obtain ∂T φ0 = −γ0 W  (φ0 ). This equation has two stable fixed points 0 and 1. So, (4.8)

for any x ∈ Ω, φ0 (x, T ) → 0 or 1 exponentially as T → ∞.

Similarly, placing our expansions into (4.2), we obtain the leading-order O(1)-equation ∇ · ε(φ0 )∇ψ0 − (φ0 − 1)2 V  (ψ0 ) = −ρf .

(4.9)

This is exactly (4.2), as expected. More equations corresponding to other higher powers of ξ can be obtained, but they provide less useful information. We therefore need to examine the next time scale, the regular time scale. We assume that, at the regular time scale, the solutions φ and ψ have the expansions φ(x, t) = φ0 (x, t) + ξφ1 (x, t) + ξ 2 φ2 (x, t) + · · · , ψ(x, t) = ψ0 (x, t) + ξψ1 (x, t) + ξ 2 ψ2 (x, t) + · · · , where φi = φi (x, t) and ψi = ψi (x, t) (i = 0, 1, . . .) are smooth and bounded functions that satisfy the boundary conditions (4.7). Note that these are functions of (x, t), different from those of (x, T ) in (4.5) and (4.6). Substituting φ and ψ in (4.1) by these expansions and matching the terms with the same orders of ξ, we obtain the first two equations

W  (φ0 ) = 0, (4.10) O ξ −1 : O(1) : (4.11)

∂t φ0 = −2P φ0 − γ0 W  (φ0 )φ1 − 2ρw (φ0 − 1)U 1 + ε (φ0 )|∇ψ0 |2 + 2(φ0 − 1)V (ψ0 ). 2

Note that (4.2) does not involve any time derivatives. Thus, substituting φ and ψ in (4.2) by their expansions, we obtain the same equation (4.9) in the leading order. Equation (4.10) implies that φ0 = 0, 1, or 1/2. By (4.8), we must have φ0 = 0 or 1. Since φ0 = 0 on ∂Ω, we have φ0 = 0 near ∂Ω. By assumption (A2) on the function U in section 2 and (4.11), φ0 cannot be identically 0 in Ω. Thus, both {φ0 = 0} and {φ0 = 1} are nonempty. But φ0 is a continuous function. This means that our

DIFFUSED SOLUTE-SOLVENT INTERFACE

2085

above expansions for φ and ψ are valid only in the outer region, i.e., the union of − + Ω+ ξ (t) := {x ∈ Ω : φ(x, t) ≈ 1} and Ωξ (t) := {x ∈ Ω : φ(x, t) ≈ 0}. The region Ωξ (t) − contains all the points xi (1 ≤ i ≤ N ). The boundary of the region Ωξ (t) contains the boundary of Ω. We assume the inner region, a thin transition layer of width O(ξ) that is a complement to the outer region, centers around the closed and smooth interface Γξ (t) = {x ∈ Ω : φ(x, t) = 1/2} enclosing Ω+ ξ (t). This interface is the perturbation of a closed and smooth interface Γ(t) that is also inside the inner region and that does not depend on ξ. We assume that all the principal curvatures of these interfaces are O(1) with respect to ξ → 0. For a given point x in the inner region, we denote by s(x, t) the signed distance from x to Γξ (t), the positive interior and negative exterior of the interface. This is a smooth function depending on ξ such that |∇s(x, t)| = 1. Let Pt x ∈ Γξ (t) be the projection of x onto Γξ (t), defined by |s(x, t)| = x − Pt x . (We use the Euclidean norm and distance.) The projection Pt depends on ξ. For ξ > 0 small enough, one can expect that the projection Pt x ∈ Γξ (t) is unique. We shall assume so. Consequently, the vector x − Pt (x) is normal to Γξ (t) at Pt x. We shall show that, to the leading order with respect to ξ → 0, the interface Γ(t), or equivalently the interface Γξ (t), does not move at this time scale. To this end, let us introduce a local coordinate system for the inner region [14, 15]. Consider an open neighborhood D(t) in R3 that covers an O(ξ)-neighborhood of a part of Γξ (t). Let Q(t) ⊂ R2 be a smooth domain and Φ(·, t) : Q(t) → D(t) ∩ Γξ (t) be a smooth parameterization of the patch of interface D(t) ∩ Γξ (t) such that the coordinate lines yi = 0 (i = 1, 2) are exactly the lines of principal curvature and the coordinates y1 and y2 of y = (y1 , y2 ) ∈ Q(t) are the corresponding arc lengths of these lines on Γξ (t). Then any point x ∈ D(t) of the inner region can be expressed uniquely as (4.12)

x = Φ(y, t) + ξz(t)n(y, t),

where Φ(y, t) is exactly the projection of x onto Γξ (t), z(t) = s(x, t)/ξ, and n(y, t) = ∇x s(x, t) is the unit normal to the interface Γξ (t) at the projection Φ(y, t), pointing from the exterior to the interior of the interface Γξ (t). We shall denote by J = J(y, t) the Jacobian matrix of the mapping Φ(·, t) : Q(t) → D(t) ∩ Γξ (t). This is a 3 × 2 matrix at each point y = (y1 , y2 ) with the jth column being the vector ∂yj Φ(y, t) (j = 1, 2). These columns are orthonormal vectors and are tangent to the surface Γξ (t). Note that all D(t), Φ(y, t), J(y, t), z(t), and n(y, t) can depend on ξ. Both Φ(y, t) and J(y, t) are O(1) with respect to ξ → 0. Note also that the z component of the (z, y) coordinate is a global variable, independent of the parameterization Φ(·, t) of the patch of interface D(t) ∩ Γξ (t). The relation (4.12) determines locally y = y(x, t) and z = z(x, t) as smooth functions of x and t. In particular, z(x, t) = s(x, t)/ξ. Let v = v(y, t) and H = H(y, t) be the normal velocity and mean curvature of the point Φ(y, t) ∈ Γξ (t) defined by v(y, t) = ∂t Φ(y, t) · n(y, t) and H(y, t) = ∇y · n(y, t)/2, respectively. We have for any x related to (z, y) by (4.12) that (4.13)

∇x z(x, t) = ξ −1 n(y, t),

(4.14)

Δx z(x, t) = 2ξ −1 H(y, t) + O(1),

(4.15)

∂t z(x, t) = −ξ −1 v(y, t),

(4.16)

∇x y(x, t) = J(y, t) + O(ξ).

Equation (4.13) follows from n(y, t) = ∇x s(x, t) and ξz(x, t) = s(x, t). Equation (4.14) follows from the fact that ∂yj n is the jth column of J multiplied by the jth principal

2086

BO LI AND YUAN LIU

curvature; cf., equation (2.3) in [15]. Fixing x and replacing z(t) by s(x, t) in (4.12), differentiating both sides of this equation with respect to t, and then dotting with the unit vector n = n(y, t), we then obtain (4.15). Finally, Equation (4.16) follows from a formula for ∇x y in Lemma 2.1 of [15]. With these, we further obtain by a series of calculations that, for any smooth functions f = f (x, t) and f˜ = f˜(z, y, t) that satisfy f (x, t) = f˜(z, y, t) with x and (z, y) related by (4.12), (4.17) (4.18) (4.19)

∇x f = (J∇y + ξ −1 n ∂z )f˜ + O(ξ),

2 f˜ + O(1), Δx f = ξ −1 2H∂z + ξ −2 ∂zz

∂t f = ∂t − ξ −1 v ∂z f˜.

We now assume that the solutions φ = φ(x, t) and ψ = ψ(x, t) have the following expansions in the intersection of D(t) and the inner region [4, 34, 36]: φ(x, t) = φ˜0 (z, y, t) + ξ φ˜1 (z, y, t) + ξ 2 φ˜2 (z, y, t) + · · · , ψ(x, t) = ψ˜0 (z, y, t) + ξ ψ˜1 (z, y, t) + ξ 2 ψ˜2 (z, y, t) + · · · , where x and (z, y) are related by (4.12), and all φ˜i = φ˜i (z, y, t) and ψ˜i = ψ˜i (z, y, t) (i = 0, 1, . . . ) are smooth and bounded functions. Substitute φ and ψ in (4.1) with these inner expansions, apply (4.17)–(4.19), and group terms of the same orders in terms of powers of ξ. Equating the terms of the order O(ξ −2 ), we get ε (φ˜0 )(∂z ψ˜0 )2 = 0. Since φ0 = 0 or 1 in the outer region, we can assume that the set of points in the inner region at which φ˜ = 0 or 1 has the Lebesgue measure zero. Thus, by our assumption (A3) in section 2, ε (φ˜0 ) = 0 a.e. in the inner region. Consequently, ∂z ψ˜0 = 0 in the inner region. With this, we obtain at the next order O(ξ −1 ) that   2 ˜ (4.20) −v ∂z φ˜0 = γ0 ∂zz φ0 − W  (φ˜0 ) . By matching the inner and outer solutions, we have φ˜0 (z, y, t) → 0 if z → −∞ and φ˜0 (z, y, t) → 1 if z → +∞. Moreover, since φ˜0 is smooth and bounded, we have also that ∂z φ˜0 → 0 as z → −∞ or ∞. Therefore, multiplying both sides of (4.20) by ∂z φ˜0 and integrating over z from −∞ and ∞, we obtain z=∞   ∞ 2  1 2 ˜ ˜ ˜ ∂z φ0 dz = γ0 (∂z φ0 ) − W (φ0 )  −v = 0. 2 −∞ z=−∞ Hence v = v(y, t) = 0. Therefore, at the leading order of ξ, the interface does not move at the regular time scale. However, in the next long time period measured by a slow time scale, the interface still moves. Note that all these results depend only on the global variable z but not on the local variable y, i.e., not on the parameterization Φ(·, t) : Q(t) → D(t) ∩ Γξ (t). 5. Relaxation dynamics: Slow time scale and sharp-interface limit. We now consider the slow time scale τ = ξt and obtain the leading order of the normal velocity for this interface motion. Based on the previous analysis, we assume at each τ the system region Ω is divided by a closed and smooth interface Γ(τ ) into two regions Ω− (τ ) and Ω+ (τ ), outside and inside of Γ(τ ), respectively. We assume the principal curvatures at any point of Γ(τ ) are of order O(1) with respect to ξ → 0. We have φ ≈ 0 in Ω− (τ ) and φ ≈ 1 in Ω+ (τ ). Moreover, all xi ∈ Ω+ (τ ) (1 ≤ i ≤ N ). Hence

DIFFUSED SOLUTE-SOLVENT INTERFACE

2087

Ω+ (τ ) approximates the solute region Ωp and Ω− (τ ) approximates the solvent region Ωw . We assume that in the outer region—the region outside an O(ξ) neighborhood of Γ(τ )—the solutions φ and ψ have the expansions φ(x, t) = φ0 (x, τ ) + ξφ1 (x, τ ) + ξ 2 φ2 (x, τ ) + · · · , ψ(x, t) = ψ0 (x, τ ) + ξψ1 (x, τ ) + ξ 2 ψ2 (x, τ ) + · · · , where φi = φi (x, τ ) and ψi = ψi (x, τ ) with τ = ξt are smooth and bounded functions that satisfy the boundary conditions (4.7). Plugging these into (4.1) and comparing the terms of the same orders of ξ, we obtain O(ξ −1 ) :

(5.1)

0 = W  (φ0 ).

Similarly, we have by (4.2) that (5.2)

O(1) :

∇ · ε(φ0 )∇ψ0 − (φ0 − 1)2 V  (ψ0 ) = −ρf .

By (5.1) and the boundary condition for φ0 , we get φ0 = 0 in the Ω− (τ ) part of the outer region and φ0 = 1 in the Ω+ (τ ) part of the outer region. In particular, as ξ → 0, φ0 converges to 1 and 0, inside and outside Γ(τ ), respectively. We now consider the solutions φ and ψ in the inner region, an O(ξ)-neighborhood of Γ(τ ). As before, we introduce a local coordinate system for the inner region. Let s(x, τ ), Pτ (x), D(τ ), Q(τ ), Φ = Φ(y, τ ), J = J(y, τ ), z = z(τ ), n = n(y, τ ), H = H(y, τ ), and v = v(y, τ ) all be the same as those in section 4 with τ replacing t and Γ(τ ) replacing Γξ (t), respectively. We relate a point x ∈ D(τ ) in the inner region to the local coordinates (z, y) by (5.3)

x = Φ(y, τ ) + ξz(τ )n(y, τ ).

We still have (4.13), (4.14), and (4.16) with t replaced by τ. Equation (4.15) is changed now to ∂t z(x, τ ) = −v(x, τ ). If f (x, τ ) and f˜(z, y, τ ) are two smooth functions and f (x, τ ) = f˜(z, y, τ ) with x and (z, y) related by (5.3), then (4.17) and (4.18) still hold true; but (4.19) is changed to ∂t f = (ξ∂τ − v ∂z ) f˜. We assume now the following local expansions in the inner region: φ(x, t) = φ˜0 (z, y, τ ) + ξ φ˜1 (z, y, τ ) + ξ 2 φ˜2 (z, y, τ ) + · · · , ψ(x, t) = ψ˜0 (z, y, τ ) + ξ ψ˜1 (z, y, τ ) + ξ 2 ψ˜2 (z, y, τ ) + · · · , where φ˜i = φ˜i (z, y, τ ) and ψ˜i = ψ˜i (z, y, τ ) (i = 0, 1, . . . ) are smooth and bounded functions, and x and (z, y) are related by (5.3). We also have by (5.3) that ˜0 (y, τ ) + O(ξ), U (x) = U (Φ(y, τ ) + ξz(τ )n(y, τ )) = U ˜0 (y, τ ) = U (Φ(y, τ )). If x = Φ(y, τ ) ∈ Γ(τ ), then U ˜0 (y, τ ) = U (x). where U Plug these expansions into (4.1) and equate the terms with the same power of ξ to obtain (5.4) (5.5)

O(ξ −2 ) : O(ξ −1 ) :

0 = ∂z ψ˜0 , 2 ˜ 0 = ∂zz φ0 − W  (φ˜0 ),

2088

BO LI AND YUAN LIU

O(1) :

  2 ˜ φ1 − W  (φ˜0 )φ˜1 −v ∂z φ˜0 = − 2P φ˜0 + γ0 2H φ˜0 + ∂zz ˜0 − 2ρw (φ˜0 − 1)U   1 + ε (φ˜0 ) |∇Γ(τ ) ψ˜0 |2 + (∂z ψ˜1 )2 + 2(φ˜0 − 1)V (ψ˜0 ). 2

(5.6)

In obtaining (5.4), we use the fact that φ0 = 0 or 1 in the outer region to assume that φ˜0 ∈ (0, 1) and hence ε (φ˜0 ) = 0 a.e. in the inner region. Equation (5.4) is used to obtain (5.5) and (5.6). By (5.4), ψ˜0 (z, y, τ ) = ψ˜0 (0, y, τ ) is the value of ψ˜0 at the point Φ(y, τ ) ∈ Γ(τ ). Note that the surface gradient term ∇Γ(τ ) ψ˜0 in (5.6) is an O(ξ)-approximation of J∇y ψ˜0 by (4.16) and the fact that columns of J are orthogonal to n: ∇Γ(τ ) ψ˜0 = (I − n ⊗ n)∇x ψ˜0 = (I − n ⊗ n)Dx y∇y ψ˜0 = (I − n ⊗ n)J∇y ψ˜0 + O(ξ) = J∇y ψ˜0 + O(ξ). Similarly, plugging the inner expansions of φ and ψ into (4.2), we obtain by (5.4) and the fact that columns of J = J(y, τ ) are orthogonal to the normal n = n(y, τ ) that   (5.7) ∂z ε(φ˜0 )∂z ψ˜1 = 0. O(ξ −1 ) : Fix x = Φ(y, τ ) ∈ Γ(τ ). Let us employ the following matching conditions [4, 34, 36]: (5.8) (5.9) (5.10)

lim φ˜0 (z, y, τ ) = φ± 0 (x, τ ),

z→±∞

lim ψ˜0 (z, y, τ ) = ψ0± (x, τ ),

ψ˜1 (z, y, τ ) = ψ1± + z∇ψ0± · ∇s (x, τ ) + o(1) z→±∞

as z → ±∞,

± where φ± 0 (x, τ ) = limα→0± φ0 (x+α∇s(x, τ )) and ψi (x, τ ) = limα→0± ψi (x+α∇s(x, τ )). Since φ0 = 0 in Ω− (τ ) and φ0 = 1 in Ω+ (τ ), we have by (5.5) and (5.8) that

(5.11)

e3z − e−3z 1 . φ˜0 (z, y, τ ) = φ˜0 (z) = + 2 2(e3z + e−3z )

We now derive the following matching condition that will be used later: (5.12)

∂z ψ˜1 (z, y, τ ) = ∇ψ0± (x, τ ) · ∇s(x, τ ) + o(1)

as z → ±∞.

Let us rewrite the matching condition (5.10) as ψ˜1 (z, y, τ ) = ψ1± (x, τ ) + z∇ψ0± (x, τ ) · ∇s(x, τ ) + g ± (z, y, τ ) for some smooth functions g + = g + (z, y, τ ) (z > 0) and g − = g − (z, y, τ ) (z < 0) such that g ± (z, y, τ ) = o(1) as z → ±∞, respectively. Thus (5.13)

∂z ψ˜1 (z, y, τ ) = ∇ψ0± (x, τ ) · ∇s(x, τ ) + ∂z g ± (z, y, τ ).

By (5.7), we have ∂z ψ˜1 (z, y, τ ) = h(y, τ )/ε(φ˜0 (z)) for some smooth and bounded function h = h(y, τ ) independent of z. This and (5.13) imply that ∂z g ± are bounded as z → ±∞. Moreover, with these two equations and (5.11), we obtain    h(y, τ )  6ε |h(y, τ )|  2 ±   2  ∂zz g (z, y, τ ) = ∂zz  ≤ max e−6z . ψ˜1 (z, y, τ ) = ∂z ε2max ε(ψ˜0 (z)) 

2089

DIFFUSED SOLUTE-SOLVENT INTERFACE

Consequently, for any A2 > A1 > 0,    A2   A2       

   2 2 2     + + + ∂z g +  ∂zz g +  dz → 0 dz  ≤ 2 ∂z ∂z g  ∂z g (A2 ) − ∂z g (A1 )  =    A1 A1 as A1 , A2 → ∞. Therefore, (∂z g + )2 → a+ as z → ∞ for some a+ ≥ 0. We must have a+ = 0, since g + → 0 as z → ∞. Therefore ∂z g + → 0 as z → ∞. Similarly ∂z g − → 0 as z → −∞. Hence (5.12) follows from (5.13). We now show that ψ0 solves the Poisson–Boltzmann equation with the dielectric boundary Γ(τ ). By (5.4), ψ˜0 is independent of z. Therefore, (5.2), the matching condition (5.9), and the boundary condition ψ0 = ψ∞ on ∂Ω allow us to determine ψ0 in the outer region. Moreover, the matching condition (5.9) implies that ψ˜0 = ψ0+ = ψ0− on Γ(τ ). Integrating both sides of (5.7) over z from −∞ to ∞, we obtain by (5.12) that εw ε0 ∇ψ0− · n = εp ε0 ∇ψ0+ · n on Γ(τ ). These, together with (5.2) and (4.7), imply that ψ0 is the solution to the elliptic interface problem ⎧ εw ε0 Δψ0 − V  (ψ0 ) = −ρf in Ω− (τ ), ⎪ ⎪ ⎪ ⎪ ε Δψ = −ρ in Ω+ (τ ), ε ⎨ p 0 0 f − + on Γ(τ ), ψ0 = ψ0 (5.14) ⎪ ⎪ on Γ(τ ), ⎪ εw ε0 ∇ψ0− · n = εp ε0 ∇ψ0+ · n ⎪ ⎩ ψ0 = ψ∞ on ∂Ω. This is equivalent to the boundary-value problem of the Poisson–Boltzmann equation (1.2) and (1.3) with ψ replaced by ψ0 [27]. We finally derive the motion law for the interface Γ(τ ). For notational simplicity, let us skip writing the variables y and τ for a moment. Multiplying both sides of (5.6) by φ˜0 and integrating over z from −∞ to ∞, we obtain that  ∞  ∞  ∞   ˜ ˜ (5.15) −vS = R1 (z)φ0 (z) dz + R2 (z)φ0 (z) dz + R3 (z)φ˜0 (z) dz, −∞

−∞

−∞

where by (5.11) 

(5.16)

 ∞ 36e12z 6x dz = dx = 1, 6z + 1)4 (e (x + 1)4 −∞ −∞ 0 ˜0 + 2(φ˜0 − 1)V (ψ˜0 ), R1 (z) := −2P φ˜0 − 2ρw (φ˜0 − 1)U   2 ˜ φ1 − W  (φ˜0 )φ˜1 , R2 (z) := γ0 2H∂z φ˜0 + ∂zz   1 R3 (z) := ε (φ˜0 ) |∇Γ(τ ) ψ˜0 |2 + (∂z ψ˜1 )2 . 2

S :=

∞

  2 φ˜0 (z) dz =

∞

˜0 and V (ψ˜0 ) are independent of z, we have Since U  ∞  ∞  ∞    ˜0 R1 (z)φ˜0 (z) dz = −2P φ˜0 φ˜0 dz + 2 V (ψ˜0 ) − ρw U φ˜0 − 1 φ˜0 dz −∞ −∞ −∞ z=∞   2 z=∞  2 ˜0 φ˜0 − 1  = −P φ˜0  + V (ψ˜0 ) − ρw U  z=−∞

(5.17)

z=−∞

= −P + ρw U − V (ψ0 ),

˜0 = U and that ψ˜0 = ψ0 by (5.4) and (5.9). where we have on Γ(τ ) that U

2090

BO LI AND YUAN LIU

By integration by parts and the fact that φ˜0 (±∞) = φ˜0 (±∞) = 0, we have  ∞  ∞  2 ˜ ˜ ∂zz φ1 φ0 dz = φ˜0 φ˜1 dz. −∞

−∞

Therefore, by (5.16) and (5.5), we have   ∞ (5.18) R2 (z)φ˜0 (z)dz = 2γ0 H + γ0 −∞

∞

−∞

  ∂z φ˜0 − W  (φ˜0 ) φ˜1 dz = 2γ0 H.

Now, since ψ˜0 = ψ0− = ψ0+ on Γ(τ ), we have ∇Γ(τ ) ψ˜0 = ∇Γ(τ ) ψ0− = ∇Γ(τ ) ψ0+ = ∇Γ(τ ) ψ0 . Consequently,   ∞  1 ∞  ˜  R3 (z)φ˜0 (z) dz = ε (φ0 ) |∇Γ(τ ) ψ˜0 |2 + (∂z ψ˜1 )2 φ˜0 dz 2 −∞ −∞  1 ∞  ˜ 1 (5.19) ε (φ0 )(∂z ψ˜1 )2 φ˜0 dz. = (εp ε0 − εw ε0 )|∇Γ(τ ) ψ0 |2 + 2 2 −∞ 2 ˜ Rewrite (5.7) as ε(φ˜0 )∂zz ψ1 = −ε (φ˜0 )φ˜0 ∂z ψ˜1 . With this and (5.12), we have z=∞  ∞   1 ∞  ˜ 1 2 ˜ ε (φ0 )(∂z ψ˜1 )2 φ˜0 dz = ε(φ˜0 )(∂z ψ˜1 )2  − ε(φ˜0 )∂z ψ˜1 ∂zz ψ1 dz 2 −∞ 2 −∞ z=−∞ 1 1 = εp ε0 |∇ψ0+ · n|2 − εw ε0 |∇ψ0− · n|2 2 2 ∞

+ −∞

ε (φ˜0 )(∂z ψ˜1 )2 φ˜0 dz.

Note that the last term on the right-hand side doubles the term on the left-hand side. This and the third equation in (5.14) then lead to  1 ∞  ˜ 1 1 ε (φ0 )(∂z ψ˜1 )2 φ˜0 dz = εw ε0 |∇ψ0− · n|2 − εp ε0 |∇ψ0+ · n|2 2 −∞ 2 2     1 1 1 εΓ(τ ) ∇ψ0 · n2 , (5.20) = − 2 εw ε0 εp ε0 where εΓ(τ ) is defined in (1.4) with Γ(τ ) replacing Γ. Finally, it follows from (5.15)–(5.20) that v = P − 2γ0 H − ρw U     1 1 1 εΓ(τ ) ∇ψ0 · n2 + 1 (εw ε0 − εp ε0 )|∇Γ(τ ) ψ0 |2 + V (ψ0 ). − + 2 εp ε0 εw ε0 2 Since our normal n points from the interior to the exterior and mean curvature H = ∇· n/2, opposite to that in the sharp-interface model (cf. the description above (1.5)), this is exactly the normal velocity vn = −δΓ F [Γ] in the sharp-interface model (cf. (1.5)). We summarize our analysis in the following Theorem 5.1. Assume as above the existence of a closed and smooth interface Γ(τ ), the outer and inner expansions for the solutions φ and ψ to the initial-boundaryvalue problem (4.1)–(4.4), and the corresponding matching conditions. Then the following are true: (1) As ξ → 0, the leading order of the phase-field function φ converges to 1 in Ω+ (τ ) and 0 in Ω− (τ ), respectively.

DIFFUSED SOLUTE-SOLVENT INTERFACE

2091

(2) The leading order of the electrostatic potential is the unique solution to the boundary-value problem of the Poisson–Boltzmann equation (1.2) and (1.3). (3) As ξ → 0, the normal velocity of the interface Γ(τ ) converges to the boundary force −δΓ F [Γ] in (1.5), with Γ(τ ) replacing Γ, as in the sharp-interface model. Acknowledgments. The authors thank Professor John Lowengrub, Dr. Jiayi Wen, and Professor Yanxiang Zhao for many helpful discussions. REFERENCES [1] D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Diffuse-interface methods in fluid mechanics, Annu. Rev. Fluid Mech., 30 (1998), pp. 139–165. [2] D. Ben-Yaakov, D. Andelman, and R. Podgornik, Dielectric decrement as a source of ion-specific effects, J. Chem. Phys., 134 (2011), 074705. [3] R. Ben´ıtez and L. Ram´ırez-Piscina, Sharp-interface projection of a fluctuating phase-field model, Phys. Rev. E, 71 (2005), 061603. [4] G. Caginalp and P. C. Fife, Dynamics of layered interfaces arising from phase boundaries, SIAM J. Appl. Math., 48 (1988), pp. 506–518. [5] Q. Cai, X. Ye, J. Wang, and R. Luo, Dielectric boundary forces in numerical Poisson– Boltzmann methods: Theory and numerical strategies, Chem. Phys. Lett., 514 (2011), pp. 368–373. [6] D. Chandler, Interfaces and the driving force of hydrophobic assembly, Nature, 437 (2005), pp. 640–647. [7] J. Che, J. Dzubiella, B. Li, and J. A. McCammon, Electrostatic free energy and its variations in implicit solvent models, J. Phys. Chem. B, 112 (2008), pp. 3058–3069. [8] H. B. Cheng, L.-T. Cheng, and B. Li, Yukawa-field approximation of electrostatic free energy and dielectric boundary force, Nonlinearity, 24 (2011), pp. 3215–3236. [9] L.-T. Cheng, J. Dzubiella, J. A. McCammon, and B. Li, Application of the level-set method to the implicit solvation of nonpolar molecules, J. Chem. Phys., 127 (2007), 084503. [10] L.-T. Cheng, B. Li, and Z. Wang, Level-set minimization of potential controlled Hadwiger valuations for molecular solvation, J. Comput. Phys., 229 (2010), pp. 8497–8510. [11] L.-T. Cheng, B. Li, M. White, and S. Zhou, Motion of a cylindrical dielectric boundary, SIAM J. Appl. Math., 73 (2013), pp. 594–616. [12] L.-T. Cheng, Y. Xie, J. Dzubiella, J. A. McCammon, J. Che, and B. Li, Coupling the level-set method with molecular mechanics for variational implicit solvation of nonpolar molecules, J. Chem. Theory Comput., 5 (2009), pp. 257–266. [13] J. B. Collins and H. Levine, Diffuse interface model of diffusion-limited crystal growth, Phys. Rev. B, 31 (1985), pp. 6119–6122. [14] S. Dai and Q. Du, Motion of interfaces governed by the Cahn–Hilliard equation with highly disparate diffusion mobility, SIAM J. Appl. Math., 72 (2012), pp. 1818–1841. [15] S. Dai and K. Promislow, Geometric evolution of bilayers under the functionalized Cahn– Hilliard equation, Proc. Royal Soc. A, 469 (2013), p. 20120505. [16] M. E. Davis and J. A. McCammon, Electrostatics in biomolecular structure and dynamics, Chem. Rev., 90 (1990), pp. 509–521. [17] A. Dembo and T. Funaki, Lectures on Probability Theory and Statistics, Lecture Notes in Math. 1869, Springer–Verlag, Berlin, 2005. [18] Q. Du, C. Liu, and X. Wang, A phase field approach in the numerical study of the elastic bending energy for vesicle membranes, J. Comput. Phys., 198 (2004), pp. 450–468. [19] J. Dzubiella, J. M. J. Swanson, and J. A. McCammon, Coupling hydrophobicity, dispersion, and electrostatics in continuum solvent models, Phys. Rev. Lett., 96 (2006), 087802. [20] J. Dzubiella, J. M. J. Swanson, and J. A. McCammon, Coupling nonpolar and polar solvation free energies in implicit solvent models, J. Chem. Phys., 124 (2006), 084905. [21] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, 2nd ed., Springer, Berlin, 1983. [22] Z. Guo, B. Li, J. Dzubiella, L.-T. Cheng, J. A. McCammon, and J. Che, Heterogeneous hydration of p53/MDM2 complex, J. Chem. Theory Comput., 10 (2014), pp. 1302–1313. [23] M. M. Hatlo, R. van Roij, and L. Lue, The electric double layer at high surface potentials: The influence of excess ion polarizability, Europhys. Lett., 97 (2012), 28010. [24] A. Karma and W. J. Rappel, Phase-field model of dendritic sidebranching with thermal noise, Phys. Rev. E, 60 (1999), pp. 3614–3625.

2092

BO LI AND YUAN LIU

[25] J. S. Langer, Models of pattern formation in first-order phase transitions, in directions in condensed matter physics, G. Grinstein and G. Mazenko, eds., World Scientific, Singapore, 1986, pp. 165–186. [26] B. Li, Continuum electrostatics for ionic solutions with nonuniform ionic sizes, Nonlinearity, 22 (2009), pp. 811–833. [27] B. Li, Minimization of electrostatic free energy and the Poisson–Boltzmann equation for molecular solvation with implicit solvent, SIAM J. Math. Anal., 40 (2009), pp. 2536–2566. [28] B. Li., X. Cheng, and Z. Zhang, Dielectric boundary force in molecular solvation with the Poisson–Boltzmann free energy: A shape derivative approach, SIAM J. Applied Math., 71 (2011), pp. 2093–2111. [29] B. Li, J. Wen, and S. Zhou, Mean-field theory and computation of electrostatics with ionic concentration dependent dielectrics, Commun. Math. Sci., to appear. [30] B. Li and Y. Zhao, Variational implicit solvation with solute molecular mechanics: From diffuse-interface to sharp-interface models, SIAM J. Appl. Math., 73 (2013), pp. 1–23. ¨ tz, and A. Voigt, Phase-field modeling of the dynamics of multi[31] J. S. Lowengrub, A. Ra component vesicles: Spinodal decomposition, coarsening, budding, and fission, Phys. Rev. E, 79 (2009), 031926. [32] K. Lum, D. Chandler, and J. D. Weeks, Hydrophobicity at small and large length scales, J. Phys. Chem. B, 103 (1999), pp. 4570–4577. [33] L. Modica, The gradient theory of phase transitions and the minimal interface criterion, Arch. Rational Mech. Anal., 98 (1987), pp. 123–142. [34] R. L. Pego, Front migration in the nonlinear Cahn–Hilliard equation, Proc. Roy. Soc. London Ser. A, 422 (1989), pp. 261–278. [35] B. Roux and T. Simonson, Implicit solvent models, Biophys. Chem., 78 (1999), pp. 1–20. [36] J. Rubinstein, P. Sternberg, and J. B. Keller, Fast reaction, slow diffusion, and curve shortening, SIAM J. Appl. Math., 49 (1989), pp. 116–133. [37] D. Shao, W.-J. Rappel, and H. Levine, Computational model for cell morphodynamics, Phys. Rev. Lett., 105 (2010), 108104. [38] K. A. Sharp and B. Honig, Electrostatic interactions in macromolecules: Theory and applications, Annu. Rev. Biophys. Chem., 19 (1990), pp. 301–332. [39] P. Sternberg, The effect of a singular perturbation on nonconvex variational problems, Arch. Rational Mech. Anal., 101 (1988), pp. 209–260. [40] R. C. Tolman, The effect of droplet size on surface tension, J. Chem. Phys., 17 (1949), pp. 333–337. [41] Z. Wang, J. Che, L.-T. Cheng, J. Dzubiella, B. Li, and J. A. McCammon, Level-set variational implicit solvation with the Coulomb-field approximation, J. Chem. Theory Comput., 8 (2012), pp. 386–397. [42] A. P. Willard and D. Chandler, The molecular structure of the interface between water and a hydrophobic substrate is liquid-vapor like, J. Chem. Phys., 141 (2014), 18C519. [43] Y. Zhao, Y. Kwan, J. Che, B. Li, and J. A. McCammon, Phase-field approach to implicit solvation of biomolecules with Coulomb-field approximation, J. Chem. Phys., 139 (2013), 024111. [44] S. Zhou, L.-T. Cheng, J. Dzubiella, B. Li, and J. A. McCammon, Variational implicit solvation with Poisson–Boltzmann theory, J. Chem. Theory Comput., 10 (2014), pp. 1454– 1467. [45] S. Zhou, Z. Wang, and B. Li, Mean-field description of ionic size effects with non-uniform ionic sizes: A numerical approach, Phys. Rev. E, 84 (2011), 021901.