Ionic Size Effects: Generalized Boltzmann ... - Semantic Scholar
Ionic Size Effects: Generalized Boltzmann Distributions, Counterion Stratification, and Modified Debye Length Bo Li
∗
Pei Liu
†
Zhenli Xu
‡
Shenggao Zhou
§
July 18, 2013
Abstract Near a charged surface, counterions of different valences and sizes cluster; and their concentration profiles stratify. At a distance from such a surface larger than the Debye length, the electric field is screened by counterions. Recent studies by a variational mean-field approach that includes ionic size effects and by Monte Carlo simulations both suggest that the counterion stratification is determined by the ionic valence-tovolume ratios. Central in the mean-field approach is a free-energy functional of ionic concentrations in which the ionic size effects are included through the entropic effect of solvent molecules. The corresponding equilibrium conditions define the generalized Boltzmann distributions relating the ionic concentrations to the electrostatic potential. This paper presents a detailed analysis and numerical calculations of such a free-energy functional to understand the dependence of the ionic charge density on the electrostatic potential through the generalized Boltzmann distributions, the role of ionic valenceto-volume ratios in the counterion stratification, and the modification of Debye length due to the effect of ionic sizes.
1
Introduction
Electrostatic interactions between macromolecules and mobile ions in an aqueous solvent generate strong, long-ranged forces that play a key role in biological processes. In such interactions, ionic sizes or excluded volumes can affect many of the detailed chemical and ∗
Department of Mathematics and NSF Center for Theoretical Biological Physics, University of California, San Diego, 9500 Gilman Drive, Mail code: 0112, La Jolla, CA 92093-0112, USA. Email: [email protected]. † Department of Mathematics and Institute of Natural Sciences, Shanghai Jiao Tong University, 800 Dongchuan Rd., Shanghai, 200240, P. R. China. Email: [email protected]. ‡ Department of Mathematics, Institute of Natural Sciences, and Ministry of Education Key Laboratory in Scientific and Engineering Computing, Shanghai Jiao Tong University, 800 Dongchuan Rd., Shanghai, 200240, P. R. China. Email: [email protected]. § Department of Mathematics and NSF Center for Theoretical Biological Physics, University of California, San Diego, 9500 Gilman Drive, Mail code: 0112, La Jolla, CA 92093-0112, USA. Email: [email protected].
1
physical properties of an underlying biological system. For instance, differences in ionic sizes can affect how mobile ions bind to nucleic acids [3, 4, 32]. The size of monovalent cations can influence the stability of RNA tertiary structures [21]. The ionic size effect is more profound in the ion channel selectivity, see, e.g., [16, 24]. Detailed Monte Carlo simulations and integral equations calculations also confirm some of these experimentally observed properties due to the non-uniformity of ionic sizes [17, 28, 34, 36]. The classical Poisson–Boltzmann (PB) equation (PBE) is perhaps the most widely used and efficient model of the electrostatics in ionic solutions, particularly in biomolecular systems [2, 10, 12, 14, 27, 30]. Despite of its many successful applications, this mean-field theory is known to fail in capturing the ionic size effects and ion-ion correlations [15, 17, 28, 35, 36]. Numerical calculations based on the classical PBE often predict unphysically high concentrations of counterions near a charged surface [5, 38]. Recently, there has been a growing interest in studying the ionic size effects within the PB-like mean-field framework [5–7, 9, 11, 18–20, 22, 23, 25, 26, 31, 33, 38]. One of the main approaches is to include the additional entropic effect of solvent molecules in the free-energy functional of all local ionic concentrations, c1 = c1 (x), . . . , cM = cM (x) (M being the total number of ionic species), similar to that in the classical PB model. Given the volume vi of an ion of the ith species (1 ≤ i ≤ M ) and the volume v0 of a solvent molecule, one can define the local solvent concentration c0 = c0 (x) by " # M X vi ci (x) . (1.1) c0 (x) = v0−1 1 − i=1
The size-modified electrostatic free-energy functional is then given by [5, 20, 22, 26] # Z " M M X X 1 −1 F [c1 , · · · , cM ] = ρψ + β ci [ln (vi ci ) − 1] − µi ci dV. 2 i=0 i=1
(1.2)
P Here, ρ = M i=1 qi ci is the ionic charge density, where qi = zi e, and zi is the valence of an ion of the ithe species and e is the elementary charge, β = (kB T )−1 with kB the Boltzmann constant and T the temperature, and µi is the chemical potential of an ion of the ith species. The electrostatic potential ψ is determined by Poisson’s equation −∇ · ε∇ψ = ρ, together with some boundary conditions, where ε is the dielectric coefficient. One easily verifies that the equilibrium conditions δci F = 0 (i = 1, . . . , M ) are given by vi ln (v0 c0 ) − ln (vi ci ) = β (qi ψ − µi ) , v0
i = 1, . . . , M,
(1.3)
where ψ is the electrostatic potential corresponding to the equilibrium ionic concentrations c1 , . . . , cM . It is shown in [22] that this system of nonlinear algebraic equations has a unique solution ci = Bi (ψ) (i = 1, . . . , M ), defining the generalized Boltzmann distributions. If v0 = v1 = · · · = vM , then these distributions are [22, 23] ci =
−βqi ψ c∞ i e
1+
PM
j=1
−βqj ψ − 1) vj c∞ j (e
2
,
i = 1, . . . , M,
where c∞ i =
vi−1 eβµi , P βµj 1+ M j=1 e
i = 1, . . . , M.
Note that, for each i, c∞ i is the value of ci when ψ = 0. In general, if v0 , v1 , . . . , vM are not all the same, then explicit formulas for ci = Bi (ψ) (i = 1, . . . , M ) seem unavailable. In studying the ionic size effects in biomolecular reactions, Lu and Zhou [26] minimized the free-energy functional (1.2) by solving for a steady-state solution of a size-modified Poisson–Nernst–Planck system of partial differential equations. In [38], Zhou et al. developed a constrained optimization method for numerically minimizing the functional (1.2) with Poisson’s equation as a constraint. They found the remarkable stratification of ionic concentrations near a highly charged surface. Moreover, they found that the ionic valenceto-volume ratio is an important parameter that determines the stratification. Specifically, the counterion with the largest valence-to-volume ratio has the highest ionic concentration near the surface, the counterion with the second largest such ratio has also the second highest ionic concentration, and so on. Subsequently, Wen et al. [37] performed Monte Carlo simulations that confirm such a phenomenon. In this work, we further study the size-modified, mean-field free-energy functional (1.2). We begin with a description of two different forms of the electrostatic free-energy functional, one with and the other without the constraint of total number of ions for each ionic species. We show the equivalence of these two forms; cf. Theorem 2.2. We also derive the conditions for the equilibrium concentrations; cf. (1.3). We then focus on the equilibrium conditions (1.3) to analyze the generalized Boltzmann distributions.PIn particular, we characterize the behavior of the equilibrium ionic charge density ρi = M i=1 qi ci as a function of ψ. Our main results of analysis are as follows: (1) The concentration of ions with the largest and that with smallest valence-to-volume ratios monotonically decreases and monotonically increases, respectively, with respect to the electrostatic potential ψ. As ψ → −∞ (or ∞), the concentration of ions with the largest (or smallest) valence-to-volume ratio approaches to the inverse of its corresponding volume, while all other concentrations approach zero. (2) The total ionic charge density ρi = ρi (ψ) is a monotonic function of the potential ψ. If all ions are cations (i.e., all zj > 0), then this density is always positive, ρi (∞) = 0, and ρi (−∞) exists and is finite. If all ions are anions (i.e., all zj < 0), then this density is always negative, ρi (−∞) = 0, and ρi (∞) exists and is finite. If both cations and anions exist, then both ρi (−∞) and ρi (∞) exist and are finite, and the density reaches zero at a unique value of the potential ψ. All these are summarized in Figure 2 in Subsection 3.2. ˆ D , is derived with a weak potential limit. A formula (3) A size-modified Debye length, λ ˆ D is given in (3.23) in Subsection 3.3. of of λ Note that ionizable groups of biomolecules can dissociate in aqueous solution to produce anions or cations. So, studying an ionic system with only cations or only anions is of interest. In general, such a system can serve a model approximation system, as near a charged surface the concentration of coions is usually very low. We finally consider a charged spherical molecule of radius R0 in an ionic solution that 3
occupies a region defined by {R0 < |x| < R∞ } for some R∞ > R0 , and study the PBE, both classical and size-modified. In the case R∞ is finite, the detailed properties of the equilibrium electrostatic potential are described through two parameters: one is the prescribed surface charge density σ at the boundary |x| = R0 and the other is an effective surface charge density σ∞ at the boundary |x| = R∞ . This effective density is set by the total charge neutrality. When σσ∞ < 0 the potential ψ is monotonic. Otherwise, it decreases then increases, or it increases then decreases. These are summarized in Theorem 4.2 and illustrated in Figure 3 in Section 4. With the same setting of a charged spherical molecule, we numerically minimize the functional (1.2) and confirm our analysis; cf. Figure 4 in Section 5. If R∞ = ∞, then only in the case that both cations and anions exist and they neutralize in the bulk does the PBE (classical or size-modified) permit a unique solution, the electrostatic potential, that is radially symmetric and monotonic, and vanishes at infinity. We remark that considering a finite R∞ is of certain interest as it is known that the effect of finite size of charged system can be significant, cf. e.g., [29], and as it is practically useful in numerical computations. The rest of this paper is organized as follows: In Section 2, we describe the electrostatic free-energy functionals and derive the conditions of equilibrium concentrations. In Section 3, we present an analysis on the equilibrium conditions and on the stratification of counterions of multiple valences and sizes near a charged surface. In Section 4, we study a charged spherical molecule immersed in an ionic solution. In Section 5, we report numerical calculations to confirm our analysis. Finally, in Section 6, we draw conclusions.
2
The Free-Energy Functional and Equilibrium Conditions
We consider an ionic solution with M (M ≥ 1) ionic species, occupying a region Ω that is a bounded, open, and connected subset of R3 with a smooth boundary Γ 6= ∅. The boundary Γ is divided into two disjoint and smooth parts ΓD and ΓN (with D for Dirichlet and N for Neumann). We assume that a surface charge density is given on ΓN and the electrostatic potential is specified on ΓD . Our set up covers the case that ΓN = ∅ and ΓD = Γ and that ΓD = ∅ and ΓN = Γ. It also covers the important case of biomolecular solvation with an implicit solvent: The region Ω is the solvent region and the part of the boundary ΓN is the solute-solvent interface; cf. Figure 1. The surface charge density on ΓN effectively replaces the density of partial charges carried by solute atoms. For each i (1 ≤ i ≤ M ), we denote by ci = ci (x) the local ionic concentration of the ith PM ionic species at position x ∈ Ω. As before, ρ = i=1 qi ci is the ionic charge density. We also denote by vi (1 ≤ i ≤ M ) the volume of an ion of the ith species and by v0 the volume of a solvent molecule. The local concentration of the solvent, c0 = c0 (x), is then defined by (1.1). For convenience we denote c = (c1 , . . . , cM ). We consider minimizing the mean-field electrostatic free-energy functional F [c] =
Z
Ω
1 ρψ dV + 2
Z
M
ΓN
X 1 σψ dS + β −1 2 i=0
Z
4
ci [ln(vi ci ) − 1] dV − Ω
M Z X i=1
µi ci dV, Ω
(2.1)
Ω
+ -
-
+
+
+
-
-
+
-
+
+ +
ΓN
-
ΓD
+
+
+ +
-
+
+
-
+ -
+
-
+ +
+
+
+ +
-
+
-
+
-
-
Figure 1: A schematic view of a charged molecule immersed in an ionic solution. The molecule is represented by the inner sphere filled with symbols +.
∇ · ε∇ψ = −ρ ψ = ψ∞ with ε ∂ψ = σ ∂n
in Ω, on ΓD ,
(2.2)
on ΓN .
Here, the first two terms in (2.1) describe the electrostatic potential energy. The electrostatic potential ψ is the unique solution of the boundary-value problem of Poisson’s equation (2.2), where ε = ε(x) (x ∈ Ω) is the dielectric coefficient assumed to be Lebesgue measurable and bounded above and below by positive constants, and n denotes the unit exterior normal at the boundary Γ. Both ψ∞ : ΓD → R and σ : ΓN → R are given bounded and smooth functions, representing the electrostatic potential at ΓD and the surface charge density at ΓN , respectively. The third term in (2.1) has the form −T S with S being the entropy. The given parameters µi (i = 1, . . . , M ) in the last term in (2.1) are the chemical potentials of the ions. We define ψD : Ω → R and ψN : Ω → R by ∇ · ε∇ψD = 0 in Ω, ∇ · ε∇ψN = 0 in Ω, ψD = ψ∞ on ΓD , ψN = 0 on ΓD , and ε ∂ψD = 0 ε ∂ψN = σ on ΓN , on ΓN , ∂n ∂n respectively. For a given function f : Ω → R that is in the Sobolev space H −1 (Ω) [1, 13], let us denote by Lf : Ω → R the unique weak solution to the following boundary-value
5
problem of Poisson’s equation ∇ · ε∇(Lf ) = −f Lf = 0 ε ∂(Lf ) = 0 ∂n
in Ω, on ΓD , on ΓN .
(The uniqueness is guaranteed if ΓD has a nonzero surface measure.) Notice that L is a linear and self-adjoint operator from the space H −1 (Ω) to that consisting of H 1 (Ω)-functions vanishing on ΓD . With these the solution ψ to the boundary-value problem (2.2) can be decomposed as ψ = ψˆ + ψD + ψN , (2.3) P M where ψˆ = L i=1 qi ci . By integration by parts, we have Z
∂ψN ˆ ε ψ dS ∂n Γ Z Z Z ˆ ˆ ψ∇ · ε∇ψN dV + ε∇ψ · ∇ψN dV = ε∇ψˆ · ∇ψN dV = Ω Ω Ω ! Z Z M X qi ci ψN dV. = −∇ · ε∇ψˆ ψN dV =
σ ψˆ dS =
ΓN
Z
∂ψN ˆ ψ dS = ε ∂n ΓN
Z
Ω
Ω
i=1
Therefore the free-energy functional F [c] can now be rewritten as ! ! Z M M M Z X X 1 X 1 F [c] = qj cj dV + qi c i L ψD + ψN − µi ci dV qi 2 Ω 2 Ω j=1 i=1 i=1 Z Z M X 1 σ (ψN + ψD ) dS. + β −1 ci [ln(vi ci ) − 1] dV + ΓN 2 i=0 Ω By formal and routine calculations using the decomposition (2.3) and the fact that L is self-adjoint, we obtain the first variations [8, 22, 23, 26] 1 vi −1 δci F [c] = qi ψ − ψD − µi + β ln(vi ci ) − ln (v0 c0 ) , i = 1, . . . , M, (2.4) 2 v0 P where ψ is determined by (2.2) with ρ = M i=1 qi ci . The equilibrium conditions δci F [c] = 0 (i = 1, . . . , M ) are then given by vi 1 ln(v0 c0 ) − ln(vi ci ) = β qi ψ − ψD − µi , i = 1, . . . , M. (2.5) v0 2 The following result can be proved by the argument used in [22]: 6
Theorem 2.1. (1) The functional F [c] is convex. Moreover, it has a unique minimizer in PM the set of all concentrations c1 , . . . , cM such that ci ≥ 0 (i = 1, . . . , M ) and i=1 vi ci ≤ 1 in Ω. (2) There exist θ1 , θ2 ∈ R with 0 < θ1 < θ2 < 1 that bound uniformly from below and above, respectively, the free-energy minimizing ionic concentrations and the corresponding concentration of solvent molecules. (3) The minimizer is characterized by the equilibrium conditions (2.5). We now consider a different formulation: minimize the electrostatic free-energy functional Z M Z X 1 1 −1 Fˆ [c] = ρψ dV + σψ dS + β ci [ln(vi ci ) − 1] dV, Ω 2 ΓN 2 Ω i=0 ∇ · ε∇ψ = −ρ in Ω, ψ = ψ∞ on ΓD , with ε ∂ψ = σ on ΓN , Z ∂n subject to ci dV = Ni , i = 1, . . . , M, Z
(2.6)
(2.7)
(2.8)
Ω
P where again ρ = M i=1 qi ci . Note that the functional F [c] defined in (2.1) has the term of the chemical potentials µi . Here the functional Fˆ [c] does not have that term. On the other hand, the minimization of Fˆ [c] is subject to the conservation of total number of ions in each of the M ionic species, defined in (2.8) in which Ni (i = 1, . . . , M ) are given positive numbers. The system (2.7) is the same as (2.2). It again follows from standard arguments that the functional Fˆ [c] with (2.7) has a unique minimizer cˆ = (ˆ c1 , . . . , cˆM ) among all c = (c1 , . . . , cM ) subject to (2.8). By a similar argument used in deriving (2.4), we have for any c = (c1 , . . . , cM ) that vi 1 −1 ˆ ln(vi ci ) − ln (v0 c0 ) , i = 1, . . . , M. (2.9) δci F [c] = qi ψ − ψD + β 2 v0 For the minimizer cˆ, we have Z δci Fˆ [ˆ c] di dV = 0 Ω
∀di :
Z
di dV = 0,
i = 1, . . . , M.
Ω
By Lemma 2.1 below, these lead to the existence of constants, still denoted by µi (i = 1, . . . , M ), such that δci Fˆ [ˆ c] = µi (i = 1, . . . , M ). By (2.9), this is exactly (2.5). We have in fact proved the following: Theorem 2.2. The variational problem of minimizing the functional (2.1) with (2.2) is equivalent to that of minimizing the functional Fˆ [c] with (2.7) and (2.8).
7
We remark that, if one knows the chemical potentials µi (i = 1, . . . , M ), then the first formulation, i.e., minimizing F [c], can be mathematically and computationally simpler than the second one that has the constraint of total numbers of ions. The question is how the chemical potentials can be estimated by experiment or other models. A similar practical issue for the second formulation with the constraint (2.8) is that one needs to have a good estimate of each Ni , the total number of ions of the ith species in the system. Lemma 2.1. Let u ∈ L2 (Ω) be such that Z Z 2 u v dV = 0 ∀v ∈ L (Ω) : v dV = 0. Ω
Ω
Then there exists a constant µ ∈ R such that u = µ almost everywhere in Ω. Proof. Denote for any w ∈ L2 (Ω) Z Z 1 w = − w dV = w dV, |Ω| Ω Ω where |Ω| is the volume (i.e., the Lebesgue measure) of Ω. Note that the integral of w − w over Ω vanishes. For any v ∈ L2 (Ω), we have by the assumption of the lemma that Z Z Z (u − u) v dV = (u − u) (v − v) dV = u (v − v) dV = 0. Ω
Ω
Ω
Therefore u = u almost everywhere in Ω. Set µ = u. The Lemma is proved.
3
Q.E.D.
Equilibrium Ionic Concentrations and Modified Debye Length
We consider the equilibrium conditions (2.5) that determine the equilibrium ionic concentrations. Introducing a shifted potential φ = ψ − ψD /2, we can rewrite (2.5) as vi ln(v0 c0 ) − ln(vi ci ) = β (qi φ − µi ) , v0
i = 1, . . . , M.
(3.1)
As c0 is given in (1.1), this is a system of nonlinear algebraic equations for the unknowns c1 , . . . , cM . In [22], Li proved the following: Theorem 3.1. (1) For each φ ∈ R the system (3.1) has a unique solution ci = Bi (φ) ∈ (0, vi−1 ) (i = 1, . . . , M ). Moreover each Bi : R → (0, vi−1 ) is a smooth function. (2) Define V : R → R by M Z φ X V (φ) = − qi Bi (ξ) dξ + V0 , (3.2) i=1
0
where V0 is a constant. Then the function V : R → R satisfies V ′′ > 0 on R and is thus strictly convex. 8
The relations ci = Bi (φ) (i = 1, . . . , M ) define the generalized Boltzmann distributions for the equilibrium concentrations c1 , . . . , cM . In general, PM explicit formulas of such distribu′ tions seem to be unavailable. Notice that −V (φ) = i=1 qi Bi (φ) with φ = ψ − ψD /2 is the ionic charge density. The proof of V ′′ > 0 in R given in [22] (cf. the proof of Lemma 5.2 in [22]) applies to our present case; cf. the remark in Subsection 3.1. The electrostatic potential ψ is now the unique solution to the boundary-value problem of the generalized and implicit PBE ψD ′ ∇ · ε∇ψ − V ψ − = 0. (3.3) 2 The boundary conditions are the same as in (2.2). The existence and uniqueness of the equilibrium concentrations ci (i = 1, . . . , M ) lead to that of the solution to this boundaryvalue problem.
3.1
Equilibrium ionic concentrations
For convenience, we denote z0 = 0, q0 = 0, and µ0 = 0. We assume that all the M ionic species have different ionic valence-to-volume ratios. We assume there are l ≥ 0 species of anions and m ≥ 0 species of cations. So M = l+m. If l = 0 then all ions are cations. If m = 0 then all ions are anions. We label all the M ionic species and the solvent concentration by −l, . . . , −1, 0, 1, . . . , m, and assume that z−1 z0 z1 zm z−l < ··· < < =0< < ··· < . (3.4) v−l v−1 v0 v1 vm For convenience, we denote J = {−l, . . . , −1, 0, 1, . . . , m}. With our notation, the equilibrium conditions (2.5) then become vi ln(v0 c0 ) − ln(vi ci ) = β (qi φ − µi ) ∀i ∈ J. (3.5) v0 Note that by Theorem 3.1 ci = Bi (φ) for any i ∈ J with i 6= 0. We thus have by (1.1) that " # X c0 = B0 (φ) = v0−1 1 − vi Bi (φ) . (3.6) i∈J,i6=0
Theorem 3.2. Let ci = Bi (φ) (i ∈ J) be the equilibrium concentrations defined by (3.5) and (3.6). Assume (3.4) holds true. (1) Monotonicity of concentrations with the largest and smallest valence-to-volume ratios. We have c′−l (φ) > 0 and c′m (φ) < 0 for all φ ∈ R. (2) Stratification of concentration profiles. If i ∈ {−l, . . . , m − 1} then (vi ci )1/vi < (vi+1 ci+1 )1/vi+1
if and only if
φ
µi−1 /vi−1 − µi /vi . qi−1 /vi−1 − qi /vi
If i ∈ {−l + 1, . . . , m} then (vi ci )1/vi > (vi−1 ci−1 )1/vi−1
9
(3) Asymptotic behavior of concentrations. We have v−l c−l → 1 and ci → 0 (i = −l + 1, . . . , m) as φ → ∞, vm cm → 1 and ci → 0 (i = −l, . . . , m − 1) as φ → −∞.
(3.7) (3.8)
Proof. We first prove that (vi ci )1/vi = e−β(µj /vj −µi /vi ) eβ(qj /vj −qi /vi )φ (vj cj )1/vj
∀i, j ∈ J.
(3.9)
In fact, it follows immediately from (3.5) that (v0 c0 )1/v0 = e−βµi /vi eβ(qi /vi )φ (vi ci )1/vi
∀i ∈ J.
(3.10)
(vi ci )1/vi = eβµi /vi e−β(qi /vi )φ (v0 c0 )1/v0
∀i ∈ J.
(3.11)
This also implies It follows from (3.10) and (3.11) that (3.9) holds true if one of i and j is 0. Now let i, j ∈ J with i 6= j, i 6= 0, and j 6= 0. We have by (3.5) that 1 βqi βµi 1 βqj βµj ln (vi ci ) + φ− = ln (vj cj ) + φ− . vi vi vi vj vj vj Consequently, (vi ci )1/vi = e−β(µj /vj −µi /vi ) eβ(qj /vj −qi /vi )φ (vj cj )1/vj . This, together with (3.10) and (3.11), implies (3.9). We now prove parts (1)–(3). (1) Denote γ = v0 c0 and ηi = eβµi for each i ∈ J. It follows from the equilibrium conditions (3.5) that ηi ci = γ vi /v0 e−βqi φ ∀i ∈ J. (3.12) vi This and the definition of c0 (cf. (1.1)) imply X ηi γ vi /v0 e−βqi φ = 1.
(3.13)
i∈J
Taking the derivative of both sides of (3.13) with respective to φ, and using our notation z0 = 0 and µ0 = 0, we obtain P βv0 γ i∈J ηi qi γ vi /v0 e−βqi φ ′ γ (φ) = P . vi /v0 e−βqi φ i∈J ηi vi γ
Fix k ∈ J. Taking the derivative of both sides of (3.12) with i = k with respect to φ and using the above expression of γ ′ (φ), we get P ηi vi (qi /vi − qk /vk ) γ vi /v0 e−βqi φ c′k (φ) = βηk γ vk /v0 e−βqk φ i∈J P . (3.14) vi /v0 e−βqi φ i∈J ηi vi γ 10
This, the assumption (3.4), and the definition qi = zi e (i ∈ J) imply that c′m (φ) < 0 and c′−l (φ) > 0. (2) This follows directly from (3.9). (3) Since 0 ≤ v−l c−l ≤ 1, we have again by (3.9) with j = −l and (3.4) that for any i ∈ J with i > −l (vi ci )1/vi ≤ e−β(µ−l /v−l −µi /vi ) eβ(q−l /v−l −qi /vi )φ → 0
as φ → ∞.
P
Since i∈J vi ci = 1, we thus also have v−l c−l → 1 as φ → ∞. This proves (3.7). The proof of (3.8) is similar. Q.E.D. Remark. By (3.14) and the Cauchy–Schwarz inequality, we obtain through a series of calculations that h P 2 P i P 2 2 β − X i∈J qi vi ci i∈J qi ci i∈J vi ci P qk c′k (φ) = ≤ 0. 2 i∈J vi ci k∈J Since qi /vi (i ∈ J) are not all the same, this is a strict inequality. Thus the function V : R → R, defined in (3.2) with −V ′ (φ) being the total ionic charge density, is a strictly convex function of φ; cf. Theorem 3.1.
3.2
Equilibrium charge density
We now study the properties of the ionic charge density V ′ (φ) = − the function V is given now by XZ φ V (φ) = − qi Bi (ξ) dξ + V0 . i∈J
P
i∈J
qi Bi (φ). By (3.2),
(3.15)
0
We distinguish three cases: (1) l = 0, m = M , and J = {0, 1, . . . , M }; (2) l = M , m = 0, and J = {−M, . . . , −1, 0}; (3) l ≥ 1, m ≥ 1, l + m = M , and J = {−l, . . . , −1, 0, 1, . . . , m}. We denote f (−∞) = lims→−∞ f (s) and f (∞) = lims→∞ f (s) if the limits exist. Theorem 3.3. Assume (3.4) holds true. (1) If l = 0 and m = M , and J = {0, 1, . . . , M }, then V ′ (−∞) = −qM /vM < 0, V ′ (∞) = 0, and V ′ (φ) < 0 for all φ ∈ R. Moreover, V (−∞) = ∞, and V (∞) = 0 with a suitably chosen constant V0 in (3.2). (2) If l = M , m = 0, and J = {−M, . . . , −1, 0}, then V ′ (−∞) = 0, V ′ (∞) = −q−M /v−M > 0, and V ′ (φ) > 0 for all φ ∈ R. Moreover, V (−∞) = 0 with a suitably chosen constant V0 in (3.2), and V (∞) = ∞. (3) If l ≥ 1, m ≥ 1, l + m = M , and J = {−l, . . . , −1, 0, 1, . . . , m}, then V ′ (−∞) = −qm /vm < 0, V ′ (∞) = −q−l /v−l > 0, and −qm /vm < V ′ (φ) < −q−l /v−l for all φ ∈ R. Moreover, minφ∈R V (φ) = 0 with a suitably chosen constant V0 in (3.2), V (−∞) = ∞, and V (∞) = ∞. 11
In Figure 2, we show the schematic behavior of the convex function V for the three different cases. V
V′
0
φ
0
φ
−qM /vM
Case 1. 0 = z0 /v0 < z1 /v1 < · · · < zM /vM . V
V′ −q−M /v−M
0
φ
0
φ
Case 2. z−M /v−M < · · · < z−1 /v−1 < z0 /v0 = 0. V
V′ −q−l /v−l
0
φ0
φ
0
φ0
φ
−qm /vm
Case 3. z−l /v−l < · · · < z−1 /v−1 < z0 /v0 = 0 < z1 /v1 < · · · < zm /vm . Figure 2: Graphs of V ′ = V ′ (φ) (left) and V = V (φ) (right) for the three different cases. The value −V ′ (φ) is the ionic charge density at a (shifted) electrostatic potential φ.
12
Proof of Theorem 3.3. (1) In this case, we have V ′ (−∞) = −qM /vM < 0 by (3.8) and V ′ (∞) = 0 by (3.7). Since V ′′ > 0 in R by Theorem 3.1, V ′ is a strictly increasing function. Therefore, V ′ (−∞) < 0 and V ′ (∞) = 0 imply that V ′ (φ) < 0 for all φ ∈ R. Fix φ0 ∈ R so that V ′ (φ0 ) < 0. Then for any φ < φ0 we have V (φ) = V (φ0 ) + V ′ (ξ)(φ − φ0 ), where ξ ∈ (φ, φ0 ). Notice that V ′ (ξ) < V ′ (φ0 ) < 0. So, V (φ) → ∞ as φ → −∞, i.e., V (−∞) = ∞. Finally, setting j = 0 in (3.9), we obtain that 1/vi 1/vi βµi /vi (vi Bi (φ)) = (vi ci ) ≤ max e e−β(qi /vi )φ . 1≤i≤M
P Therefore, the integral of M i=1 qi Bi (φ) from 0 to φ > 0 is a bounded function of φ. Thus, by the definition of V (cf. (3.2)), V (φ) is bounded below. Since V ′ < 0 on R, V decreases. Therefore, V (∞) exists. Choosing the constant XZ ∞ V0 = qi Bi (φ) dφ i∈J
0
in (3.15), we then have V (∞) = 0. (2) This part can be proved similarly. (3) The fact that V ′ (−∞) = −qm /vm < 0 and V ′ (∞) = −q−l /v−l > 0 can be proved by the same argument used for proving part (1). Since V ′′ > 0 in R, V ′ (φ) is a strictly increasing function. Therefore, −qm /vm < V ′ (φ) < −ql /vl for all φ ∈ R. Again by the same argument used in part (1), we have V (−∞) = ∞ and V (∞) = ∞. It is clear now that there exists a unique φ0 ∈ R such that V ′ (φ0 ) = 0 and V (φ0 ) = minφ∈R V (φ). Setting in (3.15), V0 =
XZ i∈J
we have minφ∈R V (φ) = 0.
3.3
φ0
qi Bi (φ) dφ, 0
Q.E.D.
Modified Debye length
We now linearize V ′ (φ) at φ = 0 to obtain a modified Debye length. We first define c∞ i = Bi (0) for any i ∈ J. Note that X v0 c∞ = 1 − vi c∞ (3.16) 0 i . i∈J,i6=0
By (3.10) with φ = 0 and ηi = eβµi , we also have that vi /v0 (v0 c∞ = ηi−1 vi c∞ 0 ) i ,
i ∈ J.
We denote again γ = v0 c0 which depends on φ. We have for |φ| ≪ 1 that γ = τ0 + τ1 φ + O φ 2 , 13
(3.17)
(3.18)
where τ0 and τ1 are two constants to be determined. Using our notation, (1.1) becomes X v0 c0 = 1 − vi ci . (3.19) i∈J,i6=0
This and (3.16) imply that τ0 = v0 c0 (0) = 1 −
X
∞ vi c∞ i = v0 c0 .
(3.20)
i∈J,i6=0
Substituting (3.18) into (3.12) for a fixed i ∈ J with i 6= 0, we obtain by (3.20) and (3.17) that v /v ηi ci = τ0 + τ1 φ + O φ2 i 0 1 − βqi φ + O φ2 vi vi /v0 τ1 ηi vi /v0 2 1+ φ+O φ 1 − βqi φ + O φ2 = τ0 vi τ0 v i τ1 ∞ 2 = ci 1 + φ+O φ 1 − βqi φ + O φ2 vτ 0 0 v i τ1 ∞ ∞ = ci + ci − βqi φ + O φ2 . (3.21) v 0 τ0
It now follows from (3.18), (3.19), (3.21), and (3.20) that X X v i τ1 ∞ ∞ 2 − βqi φ + O φ2 . vi ci − vi ci τ0 + τ1 φ + O φ = 1 − 2 ∞ v0 c0 i∈J,i6=0 i∈J,i6=0
Comparing the leading-order terms, we obtain exactly (3.20). Comparing the O(φ)-terms and noting that z0 = 0, we obtain P βv02 c∞ qi v i c ∞ 0 P i∈J2 ∞ i . τ1 = i∈J vi ci
This and (3.21) then lead to ! P ∞ v q v c i j j ∞ P j∈J 2 ∞ j − qi φ + O φ2 ci = c∞ i + βci j∈J vj cj
∀i ∈ J with i 6= 0.
Consequently,
V ′ (φ) = −
X i∈J
qi c i = −
X i∈J
" P # ∞ 2 X q v c i i i i∈J 2 ∞ 2 P − q c qi c ∞ − β . φ + O φ i i i 2 ∞ i∈J vi ci i∈J
In the small |φ| approximation, we then obtain by (3.3) and the shift φ = ψ − ψD /2 the linearized PBE with ionic size effects X 1 ˆ −2 ˆ −2 ψ = − (3.22) ∇ · ε∇ψ − ελ qi c ∞ i − ελD ψ D , D 2 i∈J 14
ˆ D > 0 is defined by where λ ˆ −2 λ D
" # P ∞ 2 q v c β X 2 ∞ i i i . = q c − Pi∈J 2 ∞ ε i∈J i i i∈J vi ci
(3.23)
The fact that the right-hand side of the above equation is nonnegative follows from the Cauchy–Schwarz inequality ! ! !2 X X X qi2 c∞ vi2 c∞ ≤ . qi v i c ∞ i i i i∈J
i∈J
i∈J
In the case that ε is a constant, ψD = 0 on ΓD , and i∈J qi c∞ i = 0 which is the charge neutrality, the linearized equation (3.22) simplifies to the familiar equation P
ˆ −2 ψ = 0. ∆ψ − λ D ˆ D the size-modified Debye length. Recall that the Debye length λD > 0 in the We call λ classical Debye–H¨ uckel theory is defined by λ−2 D =
βX 2 ∞ q c . ε i∈J i i
(3.24)
ˆ D > λD . This indicates that finite ionic sizes lead to a longer Debye It is clear that λ screening length.
4
A Charged Spherical Molecule: Basic Properties
In this section, we consider a charged spherical molecule, centered at the origin and with radius R0 > 0, immersed in a solvent that occupies the region defined by R0 < r < R∞ in the spherical coordinates for some R∞ > R0 , where R∞ can be finite or infinite. Corresponding to the general setting described in the previous section, we have in this case that Ω = {x ∈ R3 : R0 < |x| < R∞ }, ΓN = {x ∈ R3 : |x| = R0 }, ( {x ∈ R3 : |x| = R∞ } if R∞ < ∞, ΓD = {∞} if R∞ = ∞. We consider the boundary-value problem of the PBE, classical and size-modified, for the electrostatic potential ψ (cf. (2.2) and (3.3)): ε∆ψ = V ′ (ψ) in Ω, ∂ψ (4.1) =σ on ΓN , ε ∂n ψ=0 on ΓD . 15
If ψ = ψ(r) (r = |x|) is radially symmetric, then this system is the same as 2 ′′ ′ ε ψ + ψ = V ′ (ψ) if R0 < r < R∞ , r σ ′ ψ (R ) = − , 0 ε ψ(R∞ ) = 0.
(4.2)
Here, we assume both the dielectric coefficient ε > 0 and surface charge density σ are constants. We can consider the inhomogeneous Dirichlet boundary condition ψ = ψ∞ on ΓD for a constant ψ∞ by replacing ψ by ψ − ψ∞ . We assume that V : R → R is smooth, V ′′ > 0 in R, and inf s∈R V (s) = 0. We then have one and only one of the three cases covered in Theorem 3.3: V ′ > 0 in R, V ′ < 0 in R, or there exists a unique φ0 such that V (s) > V (φ0 ) = 0 for any s 6= 0. In the last case, V ′ < 0 in (−∞, φ0 ), V ′ (φ0 ) = 0, and V ′ > 0 in (φ0 , ∞). The term −V ′ (ψ) is the ionic charge density that is related to M ionic species with valences and volumes satisfying P (3.4). For the classical PBE, we have V (φ) = i∈J c∞ e−βqi φ − 1 , where c∞ ∈ J, i 6= 0) are i i (iP bulk ionic concentrations. Under the assumption on the charge neutrality i∈J qi c∞ i = 0, this V satisfies the properties in the third case in Theorem 3.3 with φ0 = 0. So our analysis covers this case. We first consider the case that R∞ < ∞. Theorem 4.1. Let ε > 0, σ ∈ R, and 0 < R0 < R∞ < ∞. Assume that V : R → R is smooth, V ′′ > 0 in R, and inf s∈R V (s) = 0. Then the boundary-value problem (4.1) has a unique weak solution. Moreover it is radially symmetric and smooth on Ω. Proof. We assume that σ 6= 0. The case that σ = 0 can be proved similarly. We denote H = {φ ∈ H 1 (R0 , R∞ ) : φ(R∞ ) = 0} and define I : H → (−∞, ∞] by Z R∞ h i ε ′ 2 I[φ] = 4π φ (r) + V (φ) r2 dr − 4πR02 σφ(R0 ) ∀φ ∈ H. (4.3) 2 R0 If φ ∈ H, then we have with δ = ε/(4σ) > 0 that Z |φ(R0 )| = leading to
R∞
Z p φ (r) dr ≤ R∞ − R0
R∞ ′
′
R0
I[φ] ≥
πεR02
Z
|φ (r) | dr R0
R∞
|φ′ (r) |2 dr − R0
2
1/2
R∞ − R0 ≤ +δ 4δ
4πσ 2 R02 (R∞ − R0 ) ε
Z
R∞
|φ′ (r) |2 dr, R0
∀φ ∈ H.
Notice that by Poincar´e’s inequality, I[φ] controls kφk2H 1 (R0 ,R∞ ) , up to an additive constant. By standard arguments of direct methods in the calculus of variations and by the fact that H 1 (R0 , R∞ ) can be compactly embedded into C([R0 , R∞ ]), we obtain the existence of a minimizer ψ ∈ H of the functional I : H → [0, ∞]. Moreover, the minimum value of I 16
over H is finite. The uniqueness of this minimizer follows from the strict convexity of I. Clearly, the minimizer ψ = ψ(r) is a smooth solution to the corresponding Euler–Lagrange equation, the first equation in (4.2). Moreover, it satisfies the natural boundary condition ψ ′ (R0 ) = −σ/ε. Since ψ ∈ H, ψ(R∞ ) = 0. If we also denote ψ(x) = ψ(|x|), then it is clear that this function is a radially symmetric solution to (4.1), and is smooth on Ω. To prove the uniqueness of solution to (4.1), we first note that the solution ψ minimizes the functional Z h Z i ε 2 Q[φ] = |∇φ| + V (φ) dV − σφ dS Ω 2 ΓN over all φ ∈ H 1 (Ω) such that φ = 0 on ΓD . Indeed, for any η ∈ H 1 (Ω) with η = 0 on ΓD , we have by integration by parts, the convexity of V , and (4.1) that Z Z h i ε 2 |∇η| + ε∇ψ · ∇η + V (ψ + η) − V (ψ) dV − ση dS Q[ψ + η] − Q[ψ] = ΓN Ω 2 Z Z ′ ≥ [ε∇ψ · ∇η + V (ψ)η] dV − ση dS Ω ΓN Z Z Z ∂ψ ′ = [−ε∆ψ + V (ψ)] η dV + ε η dS − ση dS Ω ΓN ∂n ΓN = 0.
(4.4)
Now if ψˆ is another solution to (4.1), which is not assumed to be radially symmetric, then both ψ and ψˆ are minimizers of the functional Q. Since Q is strictly convex, it has at most one minimizer. Therefore, ψˆ = ψ. Q.E.D. We now study the behavior of the equilibrium potential ψ = ψ(|x|) that is the solution to the boundary-value problem (4.1). We denote by T (r) the total amount of ionic charges in the region {x ∈ Ω : R0 < |x| < r}: Z r T (r) = −4π V ′ (ψ(s))s2 ds ∀r ∈ [R0 , R∞ ]. (4.5) R0
Then T∞ := T (R∞ ) represents the total amount of ionic charges in Ω. If we integrate both sides of Poisson’s equation in (4.1), apply integration by parts, and use the boundary condition at r = R0 , then we have 2 T∞ + 4πR02 σ + 4πR∞ εψ ′ (R∞ ) = 0.
(4.6)
This is the charge neutrality of the system, as σ∞ := εψ ′ (R∞ ) can be viewed as the effective surface charge density on the part of the boundary ΓD = {x ∈ R3 : |x| = R∞ } which is set due to the finiteness of the underlying domain. The surface charge density σ in the boundary condition and the effective surface charge density σ∞ = εψ ′ (R∞ ) are the key parameters that determine the monotonicity of the electrostatic potential ψ. Theorem 4.2. Let ε > 0, σ ∈ R, and 0 < R0 < R∞ < ∞. Assume that V : R → R is smooth, V ′′ > 0 in R, and inf s∈R V (s) = 0. Consider V ′ > 0 in R, or V ′ < 0 in R, or V ′ (φ0 ) = 0 for a unique φ0 ∈ R. Let ψ = ψ(|x|) be the solution to the boundary-value problem (4.1). 17
(1) If σ > 0 and σ∞ > 0, then there exists a unique Rc+ ∈ (R0 , R∞ ) such that ψ ′ < 0 in (R0 , Rc+ ) and ψ ′ > 0 in (Rc+ , R∞ ). If σ < 0 and σ∞ < 0, then there exists a unique Rc− ∈ (R0 , R∞ ) such that ψ ′ > 0 in (R0 , Rc− ) and ψ ′ < 0 in (Rc− , R∞ ). (2) If σσ∞ ≤ 0 and at least one of σ and σ∞ is nonzero, then ψ ′ > 0 in (R0 , R∞ ) if σ < 0 or σ∞ > 0, and ψ ′ < 0 in (R0 , R∞ ) if σ > 0 or σ∞ < 0. (3) The case that σ = σ∞ = 0 can only occur when V ′ (0) = 0. In this case ψ = 0 identically in (R0 , R∞ ).
ψ
R0
ψ
R∞
r
R0
(a) σ > 0 and σ∞ < 0.
r
(b) σ > 0 and σ∞ > 0.
ψ R0
R∞
ψ R∞
r
R0
(c) σ < 0 and σ∞ > 0.
R∞
r
(d) σ < 0 and σ∞ < 0.
Figure 3: The graph of electrostatic potential ψ.
The main results of Theorem 4.2 are illustrated in Figure 3. To prove the theorem, we first prove a lemma. Lemma 4.1. Given the assumption as in Theorem 4.2 and A := {r ∈ (R0 , R∞ ) : ψ ′ (r) = 0}. (1) If V ′ > 0 in R or V ′ < 0 in R, then either A = ∅ or A contains exactly one point. (2) If V ′ (φ0 ) = 0 at some unique φ0 ∈ R, then either A = ∅, or A contains exactly one point, or A = (R0 , R∞ ) in which case we have additionally that φ0 = 0 and ψ = 0 identically in (R0 , R∞ ). 18
Proof. Suppose there exist r1 , r2 ∈ R such that R0 ≤ r1 < r2 ≤ R∞ and ψ ′ (r1 ) = ψ ′ (r2 ) = 0. Denote ω = {x ∈ R3 : r1 < |x| < r2 }. Then, by the same argument used in proving Theorem 4.1 (cf. (4.4)) and the strict convexity of V , the restriction of ψ onto ω is the unique minimizer among all φ ∈ H 1 (ω) of the functional Z h i ε G[φ] = |∇φ|2 + V (φ) dV. (4.7) ω 2 (1) Suppose V ′ > 0. Since inf s∈R V (s) = 0, we have V (−∞) = 0. Thus, with φk = −k (k = 1, 2, . . . ), we have 0 ≤ G[ψ] ≤ G[φk ] → 0 as k → ∞. Thus G[ψ] = 0. But this is impossible: If the integral of |∇ψ|2 vanishes, then ψ is a constant. But V (s) > 0 for any s ∈ R. Thus the integral of V (ψ) over ω is positive, a contradiction to G[ψ] = 0. For this case and similarly the case that V ′ < 0, the set A can then contain at most one point. (2) In this case, V (φ0 ) = mins∈R V (s) = 0. Thus, the constant function φ = φ0 also minimizes G over H 1 (ω). By the uniqueness of minimizer which follows from the strict convexity of G, ψ = φ0 in [r1 , r2 ] and [r1 , r2 ] ⊆ A. Let s = inf A ≥ R0 and t = sup A ≤ R∞ . At any point in (s, t), ψ = φ0 and ψ ′ vanishes. Moreover, ψ ′ 6= 0 in (R0 , s) if s > R0 and in (t, R∞ ) if t < R∞ . We prove now s = R0 and t = R∞ . Assume s > R0 and ψ ′ > 0 in (R0 , s). Differentiating both sides of ∆ψ = V ′ (ψ)/ε in R0 < r < t, we obtain 2 ′′ 2 1 ′′ ′′′ ψ (r) + ψ (r) = V (ψ(r)) + 2 ψ ′ (r) ∀r ∈ (R0 , t). (4.8) r ε r Hence, setting u = ψ ′ , we have Lu := ∆u + c(r)u = 0 for R0 < r < t, where c(r) := − [(1/ε)V ′′ (ψ(r)) + 2/r2 ] < 0 and c(r) is bounded in (R0 , s). Clearly u(s) = 0 ≤ u(r) for all r ∈ (R0 , t). Hence u attains its minimum with the minimum value 0 at any interior point x with |x| = s in {x ∈ R3 : R0 < |x| < t}. By the Strong Maximum Principle (cf. Theorem 3.5 of [13]), u must be a constant in R0 < r < t. This contradicts the fact that u > 0 in R0 < r < s and u = 0 in s < r < t. Hence ψ ′ > 0 in (R0 , s) is impossible. Similarly, ψ ′ < 0 in (R0 , s) is impossible. Therefore s = R0 . Similarly, we have t = R∞ . Finally A = (s, t) = (R0 , R∞ ), and since ψ(R∞ ) = 0, ψ = φ0 = 0 in (R0 , R∞ ). Q.E.D. Proof of Theorem 4.2. (1) Assume σ > 0 and σ∞ > 0. Then ψ ′ (R0 ) < 0 and ψ ′ (R∞ ) > 0. Hence by Lemma 4.1 there exists a unique Rc+ ∈ (R0 , R∞ ) such that ψ ′ (Rc+ ) = 0. Clearly, ψ ′ < 0 in (R0 , Rc+ ) and ψ ′ > 0 in (Rc+ , R∞ ). The case that σ < 0 and σ∞ < 0 can be proved similarly. (2) Assume σ > 0 and σ∞ ≤ 0, i.e., ψ ′ (R0 ) < 0 and ψ ′ (R∞ ) ≤ 0. The other cases can be treated similarly. It suffices to show that ψ ′ does not vanish in (R0 , R∞ ). Assume there existed r0 ∈ (R0 , R∞ ) such that ψ ′ (r0 ) = 0. By Lemma 4.1, ψ ′ does not vanish at other points. Thus there are two cases. Case 1: ψ ′ < 0 in (R0 , r0 ) and ψ ′ < 0 in (r0 , R∞ ). Case 2: ψ ′ < 0 in (R0 , r0 ) and ψ ′ > 0 in (r0 , R∞ ). We show that each case is impossible. Consider Case 1. If V ′ > 0 in R, then ∆ψ > 0 on {x ∈ R3 : r0 < |x| < R∞ }. Let x0 ∈ R3 be such that |x0 | = r0 . Since ψ(|x0 |) > ψ(|x|) for any x with |x| ∈ (r0 , R∞ ). Lemma 3.4 in [13] then implies that ψ ′ (r0 ) 6= 0, a contradiction. If V ′ < 0, we obtain the same contradiction ψ ′ (r0 ) 6= 0 by considering (R0 , r0 ). Finally, assume V ′ (φ0 ) = 0. If 19
ψ(r0 ) 6= φ0 , then since ψ ′ (r0 ) = 0, ψ ′′ (r0 ) = V ′ (φ(r0 ))/ε 6= 0. These imply that ψ ′ changes its sign at r0 which is a contradiction in Case 1. Thus Case 1 is impossible. Consider Case 2. In this case, we must have ψ ′ (R∞ ) = 0, since ψ ′ (R∞ ) > 0 would lead to ψ ′ (r1 ) = 0 for some r1 ∈ (r0 , R∞ ) which is impossible. Moreover, ψ reaches its minimum at r0 and in particular ψ ′′ (r0 ) > 0. If V ′ < 0 then ∆ψ < 0 in Ω. The Strong Maximum Principle then implies that ψ must be a constant for r ∈ (R0 , R∞ ). This is impossible as ψ ′ (R0 ) < 0. If V ′ > 0, then ∆ψ > 0 and ψ(R∞ ) = 0 > ψ(r) for r ∈ (r0 , R∞ ). By Lemma 3.4 in [13], we must have ψ ′ (R∞ ) > 0, leading to a contradiction. Assume V ′ (φ0 ) = 0 at some unique φ0 . We cannot have φ0 > ψ(r0 ) as this would imply that V ′ (ψ(r0 )) < 0 and hence ψ ′′ (r0 ) < 0 which contradicts ψ ′′ (r0 ) > 0. We cannot have φ0 ≤ ψ(r0 ) neither, since otherwise we would have that V ′ (ψ) ≥ 0 and hence ∆ψ ≥ 0 on the region r0 < r < R∞ . But ψ(R∞ ) = 0 > ψ(r) for any r ∈ (r0 , R∞ ). Lemma 3.4 in [13] then implies ψ ′ (R∞ ) > 0, a contradiction. Thus Case 2 is also impossible. (3) We have now ψ ′ (R0 ) = 0 and ψ ′ (R∞ ) = 0. As in the proof of Lemma 4.1 with r1 = R0 and r2 = R∞ , ψ is a constant in [R0 , R∞ ]. Since ψ(R∞ ) = 0, ψ = 0 identically in [R0 , R∞ ]. Thus V ′ (0) = V ′ (ψ) = ε∆ψ = 0. Hence φ0 = 0. Q.E.D. We recall that for any r ∈ (R0 , R∞ ), T (r) is the total amount of ionic charges in the region {x ∈ Ω : R0 < |x| < r}; cf. (4.5). In the following corollary, we consider a special and useful case but deliver more results: Corollary 4.1. Let ε > 0, σ ∈ R, and 0 < R0 < R∞ < ∞. Assume that V : R → R is smooth, V ′′ > 0 in R, and V (s) > V (0) = 0 for any s ∈ R with s 6= 0. Let ψ = ψ(|x|) be the solution to the boundary-value problem (4.1). (1) If σ = 0, then ψ = 0 and T = 0 identically. (2) If σ > 0, then ψ ′ < 0, ψ ′′ > 0, and T ′ < 0 in (R0 , R∞ ). (3) If σ < 0, then ψ ′ > 0, ψ ′′ < 0, and T ′ > 0 in (R0 , R∞ ). Proof. (1) It is easy to verify that in this case ψ = 0 is the unique solution to (4.1). Moreover, T ′ (r) = −V ′ (ψ(r)) = −V ′ (0) = 0 and T (R0 ) = 0. Hence T = 0 in (R0 , R∞ ). (2) Assume σ > 0. We then have σ∞ = εψ ′ (R∞ ) ≤ 0. Otherwise, σ∞ > 0. Then by Theorem 4.2 there exists Rc ∈ (R0 , R∞ ) such that ψ ′ < 0 in (R0 , Rc ) and ψ ′ > 0 in (Rc , R∞ ). We have ψ(Rc ) < min(ψ(R0 ), ψ(R∞ )) ≤ 0. Choose τ ∈ R such that ψ(Rc ) < τ < ˆ < I[ψ], where min(ψ(R0 ), ψ(R∞ )). Define ψˆ by ψˆ = ψ if ψ ≥ τ and ψˆ = τ if ψ < τ. Then I[ψ] I[·] is defined in (4.3). This is impossible as ψ minimizes I[·] over H = {φ ∈ H 1 (R0 , R∞ ) : ˆ Now, since σ > 0 and σ∞ ≤ 0, we have ψ ′ < 0 by Theorem 4.2. φ(R∞ ) = 0} that includes ψ. Since ψ(R∞ ) = 0, ψ > 0 in (R0 , R∞ ). Hence V ′ (ψ) > 0 in (R0 , R∞ ). Consequently, 2 1 ψ ′′ (r) = − ψ ′ (r) + V ′ (ψ(r)) > 0 in (R0 , R∞ ). r ε Since V ′ (ψ) > 0 in (R0 , R∞ ), we also have T ′ (r) < 0 for all r ∈ (R0 , R∞ ). (3) This can be proved similarly. Q.E.D. We now consider the case R∞ = ∞. Theorem 4.3. Let ε > 0, σ ∈ R, and R∞ = ∞. 20
(1) Suppose V : R → R is smooth and V ′ (0) 6= 0. Then there is no solution to the boundary-value problem (4.1). (2) Suppose V : R → R is smooth, V ′′ > 0 on R, V (−∞) = ∞ and V (∞) = ∞, and V (0) = inf s∈R V (s) = 0. Then there exists a unique solution to the boundary-value problem (4.1). Moreover, this solution ψ = ψ(|x|) is radially symmetric and smooth on Ω, ψ ′ < 0 and ψ ′′ > 0 in r = |x| if σ > 0, ψ ′ > 0 and ψ ′′ < 0 in r = |x| if σ > 0, and ψ = 0 for all r if σ = 0. We remark that in general all charges in an ionic system are in a finite region. Therefore it is natural to impose the far-field condition ψ(R∞ ) = 0. The properties on V assumed in part (2) of the theorem are exactly those covered in the third case in Theorem 3.3 but with the minimum of V attained exactly at 0. The theorem states that the solution to (4.1) exists only when V ′ (0) = 0. This is precisely the condition of charge neutrality in the bulk. Therefore, for the three cases of V stated in Theorem 3.3, only the third case with V ′ (0) = 0 permits a solution to (4.1). Proof of Theorem 4.3. (1) Assume V ′ (0) > 0. (The case that V ′ (0) < 0 is similar.) Assume ψ is a solution to (4.1), which is not assumed to be radially symmetric. Clearly, ψ is smooth on Ω. Since ψ(∞) = 0 and V ′ (0) > 0, there exist A ≥ R0 and α > 0 such that ε∆ψ = V ′ (ψ) ≥ α for all x ∈ Ω with |x| ≥ A. In the spherical coordinates (r, θ, ϕ) ˆ θ, ϕ), and (r ≥ 0, 0 ≤ θ ≤ π, 0 ≤ ϕ < 2π), ψ = ψ(r, ! ! ˆ ˆ 1 ∂ ∂ 2 ψˆ ∂ ψ 1 ∂ ψ 1 ∂ ∆ψ = 2 ≥α ∀r ≥ A. (4.9) r2 + 2 sin θ + 2 2 r ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂ϕ2 Define
1 u(r) = 4π
Z
2π 0
Z
π
ˆ θ, ϕ) dθdϕ sin θψ(r,
∀r > R0 .
0
Clearly, u(∞) = 0, since ψˆ → 0 as r → ∞. Multiplying both sides of the inequality in (4.9) by sin θ/(4π), and then integrating the resulting sides over θ ∈ [0, π] and ϕ ∈ [0, 2π), we ′ obtain (1/r2 ) (r2 u′ (r)) ≥ α for all r ≥ A. This leads to 1 1 3 1 ′ 2 ′ u (r) ≥ αr + A u (A) − αA ∀r ≥ A. 3 3 r2 With one more time of integration, we get 1 1 3 1 1 2 ′ 2 2 u(r) ≥ u(A) + α r − A − A u (A) − αA − 6 3 r A
∀r ≥ A.
This contradicts u(∞) = 0. Part (1) is proved. (2) We first show the uniqueness. Suppose ψ1 and ψ2 are both solutions to (4.1). Clearly they are smooth on Ω = {x ∈ R3 : |x| ≥ R0 }. Let η = ψ1 − ψ2 . This function attains its maximum on Ω. Suppose P := maxx∈Ω η > 0. If there exists x0 ∈ Ω = {x ∈ R3 : |x| > R0 } such that η(x0 ) = P > 0, then there exists an open ball B(x0 ) ⊂ Ω centered at x0 such that 21
η > 0 on B(x0 ). Since ψ1 and ψ2 are bounded on Ω, ψ1 (x), ψ2 (x) ∈ K for some compact subset of R for all x ∈ Ω. Setting γ0 = mins∈K V ′′ (s) > 0, we obtain ε∆η = V ′ (ψ1 ) − V ′ (ψ2 ) ≥ γ0 η > 0
in B(x0 ).
(4.10)
Since η(x0 ) = maxx∈B(x0 ) η > 0, the Strong Maximum Principle (cf. Theorem 3.5 of [13]) then implies that η is a constant in B(x0 ). Consequently ∆η = 0 in B(x0 ), and (4.10) implies η = 0 in B(x0 ), contradicting to η(x0 ) > 0. Thus there exists y0 ∈ ∂Ω with |y0 | = R0 such that η(y0 ) = P > η(x) for all x ∈ Ω. But this and (4.10) imply ∂n η(y0 ) = 0 by Lemma 3.4 in [13], contradicting ∂n η = 0 on ΓD = {x ∈ R3 : |x| = R0 }. Hence η ≤ 0 and ψ1 ≤ ψ2 in Ω. Similarly ψ2 ≤ ψ1 in Ω. Thus ψ1 = ψ2 in Ω. We now prove the existence of solution to (4.1) that satisfies all the properties. If σ = 0 then ψ = 0 in Ω is the desired solution. Assume σ > 0. (The case σ < 0 is similar.) Choose a sequence Rk ∈ R (k = 1, 2, . . . ) such that R0 < R1 < R2 < · · · < Rk < · · · and Rk → ∞ as k → ∞. For each integer k ≥ 1, let ψk = ψk (|x)|) be the unique solution to the boundary-value problem ε∆ψk = V ′ (ψk ) in Ωk = {x ∈ R3 : R0 < |x| < Rk }, ∂ψk (4.11) ε =σ on ΓkN = {x ∈ R3 : |x| = R0 }, ∂n ψk = 0 on ΓkD = {x ∈ R: |x| = Rk }. By Corollary 4.1, ψk > 0, ψk′ < 0, and ψk′′ > 0 in (R0 , Rk ). Using the spherical coordinates, we have by the first equation in (4.11) that 2 ψk′′ (r) + ψk′ (r) = V ′ (ψk (r)) in (R0 , Rk ). r
(4.12)
Multiply both sides of this equation by ψk′ (r) and integrate the resulting sides over [R0 , Rk ]. Noting that ψk′ ψk′′ = (1/2)(ψk′2 )′ and ψk′ V ′ (ψk ) = (d/dr)V (ψk ), we then obtain by ψk′ (R0 ) = −σ/ε, ψk (Rk ) = 0 and V (0) = 0 that σ 2 Z Rk 2 1 2 ′ 2 (ψk (Rk )) − + (ψk′ (r)) dr = −V (ψk (R0 )). 2 ε R0 r Thus, since V (s) > 0 for any s > 0, we have 0 ≤ V (ψk (r)) ≤ σ 2 /(2ε2 ) for all r ∈ [R0 , Rk ] and k ≥ 1. Since V is positive and strictly increasing in (0, ∞) and V (∞) = ∞, and since ψk (R0 ) ≥ ψk (r) for all r ∈ [R0 , Rk ], we have by setting φc > 0 with V (φc ) = σ 2 /(2ε2 ) that 0 ≤ ψk (r) ≤ φc for all r ∈ [R0 , Rk ] and all k ≥ 1. Moreover, since ψk′′ > 0 in (R0 , Rk ), −σ/ε = ψk′ (R0 ) ≤ ψk′ (r) ≤ ψk′ (R∞ ) ≤ 0 for all r ∈ [R0 , Rk ] and all k ≥ 1. Define ψk = 0 in 1,∞ (Rk , ∞) (k = 1, 2, . . . ). Then the sequence {ψk }∞ (R0 , ∞). k=1 is bounded in W ∞ ∞ There now exists a subsequence {ψk,1 }k=1 of {ψk }k=1 such that ψk,1 → ψ in H 1 (R0 , R1 ) ∞ as k → ∞ for some ψ ∈ H 1 (R0 , R1 ). There exists a subsequence {ψk,2 }∞ k=1 of {ψk,1 }k=1 such 1 1 that ψk,2 → φ in H (R0 , R2 ) as k → ∞ for some φ ∈ H (R0 , R2 ) with φ = ψ on (R0 , R1 ). We define ψ = φ in (R1 , R2 ) so that ψ ∈ H 1 (R0 , R2 ). Repeating, we obtain ψ ∈ C([R0 , ∞)) with 22
∞ ψ ∈ H 1 (R0 , R) for any R > R0 and subsequences {ψk,j }∞ k=1 of {ψk }k=1 (j = 1, 2, . . . ), where ∞ for each j ≥ 1, the sequence {ψk,j+1 }∞ k=1 is a subsequence of {ψk,j }k=1 . Moreover, ψk,j → ψ in 1 ∞ ∞ H (R0 , Rj ) as k → ∞. Now {ψk,k }k=1 is a subsequence of {ψk }k=1 , and {ψk,k }∞ k=1 converges to ψ in H 1 (R0 , R), and hence also in C([R0 , R]), for any R > R0 . The limit ψ = ψ(|x|) is clearly radially symmetric. It is a (weak and hence strong) solution to the first two equations in (4.1). We show now that limr→∞ ψ(r) = 0. For each k ≥ 1, define −1 1 1 ψ (R ) 1 − 1 if R0 ≤ r ≤ Rk , − k 0 R0 Rk r Rk ηk (r) = 0 if r > Rk .
One easily verifies that ηk (|x|) satisfies ε∆ηk = 0 in Ωk (cf. (4.11) for the definition of Ωk ) and ηk = ψk on the boundary of Ωk . Since ε∆(ψk − ηk ) ≥ 0 in Ωk , we thus have by the Maximum Principle that ψk ≤ ηk in Ωk and hence in Ω = {x ∈ R3 : |x| > R0 }. Denote ∞ ∞ {ηk,k }∞ k=1 the subsequence of {ηk }k=1 that corresponds to {ψk,k }k=1 . Then we have 0 ≤ lim sup ψ(r) = lim sup lim ψk,k (r) ≤ lim sup lim ηk,k (r) = lim sup r→∞
r→∞
k→∞
r→∞
k→∞
r→∞
R0 ψ(R0 ) = 0, r
implying limr→∞ ψ(r) = 0. Since ψk′ ≤ 0 on (R0 , Rk ) for all k ≥ 1 ψ ′ ≤ 0 on (R0 , ∞). Let u(r) = ψ ′ (r). Then u(|x|) satisfies ∆u − c(r)u = 0 on {x ∈ R3 : |x| > R0 }, where c(r) = − [(1/ε)V ′′ (ψ(r)) + 2/r2 ] < 0, cf. (4.8). Since u ≤ 0 in (R0 , ∞), the Strong Maximum Principle implies that u 6= 0 in (R0 , ∞). Hence u = ψ ′ < 0 in (R0 , ∞). Since ψ ≥ 0 and V (ψ) ≥ 0, we also have ψ ′′ (r) = −(2/r)ψ ′ + V ′ (ψ) > 0 for all r > R0 . Q.E.D.
5
Numerical Results
We now present numerical results to illustrate our analysis. We consider minimizing numerically the free-energy functional (2.6) with (2.7) and (2.8). We choose Ω = {x ∈ R3 : R0 < |x| < R∞ },
ΓN = {x ∈ R3 : |x| = R0 },
ΓD = {x ∈ R3 : |x| = R∞ },
with R0 , R∞ ∈ R such that R0 < R∞ . We assume the dielectric coefficient ε and surface charge density σ are both constants. We also set ψD = 0 on ΓD . Note that our system is radially symmetric. We use the constrained optimization method developed in our previous work [38] to solve numerically our optimization problem. We set R0 = 10 ˚ A and R∞ = 80 ˚ A. We distribute uniformly some negative charges on the 2 sphere r = R0 with the surface charge density σ = −150/(4π 2 R02 ) e/˚ A . We consider three species of counterions K+ , Ca2+ , Al3+ , and one species of coions Cl− in the solution that occupies Ω. So M = 4. We label these ionic species by 1, 2, 3, and 4, respectively. We label the solvent (water) by 0. The corresponding valences and ionic volumes are (z1 , z2 , z3 , z4 ) = 3 A , respectively. The volume (+1, +2, +3, −1) and (v1 , v2 , v3 , v4 ) = (5.513 , 4.753 , 4.13 , 6.373 ) ˚ 23
˚3 . For these real ionic systems, it happens of a solvent (water) molecule is v0 = 2.753 A that α4 < α0 = 0 < α1 < α2 < α3 , where αk = zk /vk (0 ≤ k ≤ 4) with z0 = 0. We use a3 to approximate the volume of an ion of diameter a. The temperature is T = 300K. The dielectric coefficient is absorbed into the Bjerrum length lB = e2 /(4πεkB T ). We set the Bjerrum length to be lB = 7˚ A. We consider two cases of the numbers of ions in different species. Case I: (N1 , N2 , N3 , N4 ) = (60, 30, 20, 50). By simple calculations, we have that in this case the effective surface charge density σ∞ > 0. Case II: (N1 , N2 , N3 , N4 ) = (50, 50, 50, 50). In this case σ∞ < 0. In Figure 4, we plot our computational results of the electrostatic potential ψ and the equilibrium ionic concentrations. For Case I, the potential ψ is monotonically increasing. This is exactly predicted by Theorem 4.2, since in this case, σ < 0 and σ∞ > 0. Notice that the concentration of +3, which has the largest ratio of valence-to-volume among all the species, monotonically deceases, while that of −1 monotonically increases. These are predicted by Theorem 3.2. For Case II, σ < 0 and σ∞ > 0, and the behavior of the potential ψ is again as predicted by Theorem 4.2. We observe the stratification of the counterion concentrations in both cases.
6
Conclusions
We have analyzed a PB-like mean-filed model for electrostatic interactions of multivalent ions with finite ionic sizes in an ionic solution near a charged surface. Our work continues that presented in [22] but gives more details on how equilibrium ionic concentrations depend on the electrostatic potential and how the important parameters, the valence-to-volume ratios of ions, determine the properties of ionic concentrations near the charged surface. Here we discuss some of our main results. First, our analysis shows that for a strong potential, the layering structure of counterion concentrations does depend on the valenceto-volume ratios. In general, the comparison of concentrations for different ionic species can be subtle. Instead of comparing ci and cj , one may need to compare (vi ci )1/vi and (vj cj )1/vj . The difference of these quantities may only be observed for certain range of the potential. See part (2) of Theorem 3.2. Second, we characterize the ionic charge density through the function V defined in (3.2); cf. also (3.15). We recall for the classical PB theory that the corresponding function V is defined by X −βqi φ V ′ (φ) = − qi c ∞ , i e i∈J,i6=0
c∞ i
where (i ∈ J, i 6= 0) are positive constants. An important difference here is that, when the ionic size effects are included, the function V grows linearly at the −∞ and ∞, rather than exponentially as in the classical PB theory. Third, when the ionic size effects are included, the Debye length is then modified. In ˆ D defined in (3.23) with Table 1 and Table 2, we compare the size-modified Debye length λ the classical Debye length λD defined in (3.24) for a few ionic systems. The difference between these two tables is that the diameters of ions listed in Table 1 do not, while those 24
0.08 25
0
Concentrations (M)
Electrostatic p otential (Volt)
0.1
−0.1 −0.2 −0.3 −0.4 −0.5
0.06
0.04
20
0.02 15 0 65
70
75
80
+1 +2 +3 −1
10
5
−0.6
10
20
30
40
50
60
70
0 10
80
11
12
13
14
15
16
17
18
19
20
Radial distance ( ˚ A)
Radial distance ( ˚ A)
0.08 25
0
Concentrations (M)
Electrostatic p otential (Volt)
0.1
−0.1
−0.2
−0.3
0.04
20
0.02 15 0 65
70
75
80
+1 +2 +3 −1
10
5
−0.4
−0.5 10
0.06
20
30
40
50
60
70
0 10
80
11
12
13
14
15
16
17
18
19
20
Radial distance ( ˚ A)
Radial distance ( ˚ A)
Figure 4: The electrostatic potential ψ (left) and the equilibrium ionic concentrations (right). Top. Case I: (N1 , N2 , N3 , N4 ) = (60, 30, 20, 50), σ < 0, and σ∞ > 0. Bottom. Case II: (N1 , N2 , N3 , N4 ) = (50, 50, 50, 50), σ < 0, and σ∞ < 0. The concentrations profiles are marked by the corresponding ionic valences.
25
in Table 2, include the size of a hydration shell. It is clear that the modification of the Debye length due to the inclusion of the ionic size effect is not significant. This may suggest that a linearized system may not describe well the size effect. Clearly, further tests are needed on real ionic solutions and biomolecular systems. Table 1: A comparison of the Debye length λD and the ionic size-modified Debye length ˆ D for a few systems of ionic solutions. The diameter of an ion does not include that of a λ water molecule (hydration shell). Concentrations mean bulk concentrations. Ionic systems Cl− , Na+ Cl− , K+ Cl− , Ca2+ Cl− , Na+ , K+ Cl− , Na+ , Ca2+ Cl− , Na+ , K+ , Ca2+ Cl− , K+ , Ca2+ , Al3+
Diameters 3.62, 2.04 3.62, 2.76 3.62, 2.00 3.62, 2.04, 3.62, 2.04, 3.62, 2.04, 3.62, 2.76,
of ions (˚ A)
2.76 2.00 2.76, 2.00 2.00, 1.35
Concentrations (M) 0.1, 0.1 0.1, 0.1 0.2, 0.1 0.2, 0.1, 0.1 0.3, 0.1, 0.1 0.4, 0.1, 0.1, 0.1 0.6, 0.1, 0.1, 0.1
λD ( ˚ A) 9.715 9.715 5.609 6.870 4.858 4.345 3.072
ˆ D (˚ λ A) 9.726 9.720 5.618 6.880 4.870 4.358 3.085
Table 2: A comparison of the Debye length λD and the ionic size-modified Debye length ˆ D for a few systems of ionic solutions. The diameter of an ion includes that of a water λ molecule (hydration shell). Concentrations mean bulk concentrations. Ionic systems Cl− , Na+ Cl− , K+ Cl− , Ca2+ Cl− , Na+ , K+ Cl− , Na+ , Ca2+ Cl− , Na+ , K+ , Ca2+ Cl− , K+ , Ca2+ , Al3+
Diameters 6.37, 4.79 6.37, 5.51 6.37, 4.75 6.37, 4.79, 6.37, 4.79, 6.37, 4.79, 6.37, 5.51,
of ions (˚ A)
5.51 4.75 5.51, 4.75 4.75, 4.10
Concentrations (M) 0.1, 0.1 0.1, 0.1 0.2, 0.1 0.2, 0.1, 0.1 0.3, 0.1, 0.1 0.4, 0.1, 0.1, 0.1 0.6, 0.1, 0.1, 0.1
λD ( ˚ A) 9.715 9.715 5.609 6.870 4.858 4.345 3.072
ˆ D (˚ λ A) 9.847 9.763 5.701 6.970 4.974 4.449 3.172
Finally, for a spherical charged molecule immersed in an ionic solution occupying a finite region, the qualitative behavior of the electrostatic potential—either it is monotonic or it changes the monotonicity only once in the entire range of the radial variable—is completely determined by the sign of two surface charge densities. One is the prescribed surface charge density and the other is the effective surface charge density. If the region of ionic solution is the entire complement of the sphere in R3 , then the ionic charge neutrality is necessary for the existence and uniqueness of the electrostatic potential that vanishes at infinity, governed by the PBE (classical or size-modified). Such a potential is always monotonic. 26
We conclude with two remarks on the mathematical model we have studied here. First, for the case of a uniform size for all ions and solvent molecules, the free-energy functional (2.1) can be derived from a lattice-gas model; cf. [20]. For a general and more interesting case where the sizes are different, such a derivation seems not available. It is therefore interesting to give a rigorous derivation of such a free-energy functional using the notion of equilibrium statistical mechanics. Second, we have assumed mainly that the region Ω of ionic solution is bounded. The case that Ω has an infinite volume is a tricky situation to define the functional F [c]; cf. (2.1) and (2.6). Usually, the concentration ci is taken to be the bulk value c∞ i (> 0) at infinity (or far away from the charged surface). Such a value can be experimentally determined. If Ω is unbounded with an infinite volume, then the integral of ci over Ω will be just infinite if ci = c∞ i at infinity. The PBE, however, can be defined on an unbounded domain as the electrostatic potential ψ can exist in the entire space.
Acknowledgment. This work was supported by the National Science Foundation (NSF) through grant DMS0811259 (B.L.), the NSF Center for Theoretical Biological Physics (CTBP) through grant PHY-0822283 (B.L.), the National Institutes of Health through the grant R01GM096188 (B.L.), the Natural Science Foundation of China through grant NSFC-11101276 and NSFC91130012 (P.L. and Z.X.), and the Alexander von Humboldt Foundation (Z.X.). The authors thank the referee for helpful comments.
References [1] R. Adams. Sobolev Spaces. Academic Press, New York, 1975. [2] D. Andelman. Electrostatic properties of membranes: The Poisson–Boltzmann theory. In R. Lipowsky and E. Sackmann, editors, Handbook of Biological Physics, volume 1, pages 603–642. Elsevier, 1995. [3] Y. Bai, K. Travers, V. B. Chu, J. Lipfert, S. Doniach, and D. Herschlag. Quantitative and comprehensive decomposition of the ion atmosphere around nucleic acids. J. Amer. Chem. Soc., 129:14981–14988, 2007. [4] M. L. Bleam, C. F. Anderson, and M. T. J. Record. Relative binding affinities of monovalent cations for double-stranded DNA. Proc. Natl. Acad. Sci. USA., 77:3085– 3089, 1980. [5] I. Borukhov, D. Andelman, and H. Orland. Steric effects in electrolytes: A modified Poisson–Boltzmann equation. Phys. Rev. Lett., 79:435–438, 1997. [6] I. Borukhov, D. Andelman, and H. Orland. Adsorption of large ions from an electrolyte solution: A modified Poisson–Boltzmann equation. Electrochimica Acta, 46:221–229, 2000. 27
[7] A. H. Boschitsch and P. V. Danilov. Formulation of a new and simple nonuniform size-modified Poisson–Boltzmann description. J. Comput. Chem., 33:1152–1164, 2012. [8] 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:3058–3069, 2008. [9] V. B. Chu, Y. Bai, J. Lipfert, D. Herschlag, and S. Doniach. Evaluation of ion binding to DNA duplexes using a size-modified Poisson–Boltzmann theory. Biophys. J, 93:3202– 3209, 2007. [10] M. E. Davis and J. A. McCammon. Electrostatics in biomolecular structure and dynamics. Chem. Rev., 90:509–521, 1990. [11] B. Eisenberg, Y. K. Hyon, and C. Liu. Energy variational analysis of ions in water and channels: Field theory for primitive models of complex ionic fluids. J. Chem. Phys., 133:104104, 2011. [12] F. Fogolari, A. Brigo, and H. Molinari. The Poisson–Boltzmann equation for biomolecular electrostatics: a tool for structural biology. J. Mol. Recognit., 15:377–392, 2002. [13] D. Gilbarg and N. S. Trudinger. Elliptic Partial Differential Equations of Second Order. Springer-Verlag, 2nd edition, 1998. [14] P. Grochowski and J. Trylska. Continuum molecular electrostatics, salt effects and counterion binding—A review of the Poisson–Boltzmann model and its modifications. Biopolymers, 89:93–113, 2008. [15] A. Y. Grosberg, T. T. Nguyen, and B. I. Shklovskii. Colloquium: The physics of charge inversion in chemical and biological systems. Rev. Mod. Phys., 74:329–345, 2002. [16] B. Hille. Ion Channels of Excitable Membranes. Sinauer Associates, 3rd edition, 2001. [17] J. J. Howard, J. S. Perkyns, and B. M. Pettitt. The behavior of ions near a charged wall–dependence on ion size, concentration, and surface charge. J. Phys. Chem. B, 114:6074–6083, 2010. [18] M. S. Kilic, M. Z. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes at large applied voltages. I. Double-layer charging. Phys. Rev. E, 75:021502, 2007. [19] M. S. Kilic, M. Z. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes at large applied voltages. II. Modified Poisson–Nernst–Planck equations. Phys. Rev. E, 75:021503, 2007. [20] V. Kralj-Igliˇc and A. Igliˇc. A simple statistical mechanical approach to the free energy of the electric double layer including the excluded volume effect. J. Phys. II (France), 6:477–491, 1996.
28
[21] D. Lambert, D. Leipply, R. Shiman, and D. E. Draper. The influence of monovalent cation size on the stability of RNA tertiary structures. J. Mol Biol., 390:791–804, 2009. [22] B. Li. Continuum electrostatics for ionic solutions with nonuniform ionic sizes. Nonlinearity, 22:811–833, 2009. [23] B. Li. Minimization of electrostatic free energy and the Poisson–Boltzmann equation for molecular solvation with implicit solvent. SIAM J. Math. Anal., 40:2536–2566, 2009. [24] S. W. Lockless, M. Zhou, and R. MacKinnon. Structural and thermodynamic properties of selective ion binding in a K + channel. PLoS Biology, 5:1079–1088, 2007. [25] J. J. L´opez-Garcia and Jos´e Horno. Poisson–Boltzmann description of the electrical double layer including ion size effects. Langmuir, 27:13970–13980, 2011. [26] B. Z. Lu and Y. C. Zhou. Poisson–Nernst–Planck equations for simulating biomolecular diffusion-reaction processes II: Size effects on ionic distributions and diffusion-reaction rates. Biophys J., 100:2475–2485, 2011. [27] A. C. Maggs. A minimizing principle for the Poisson–Boltzmann equation. EPL, 98:16012, 2012. ´ [28] M. Quesada-P´erez, A. Mart´ın-Molina, and R. Hidalgo-Alvarez. Simulation of electric double layers with multivalent counterions: Ion size effect. J. Chem. Phys., 121:8618– 8626, 2004. [29] R. B. Schoch, J. Han, and P. Renaud. Transport phenomena in nanofluidics. Rev. Mod. Phys., 80:839–883, 2008. [30] K. A. Sharp and B. Honig. Electrostatic interactions in macromolecules: Theory and applications. Annu. Rev. Biophys. Biophys. Chem., 19:301–332, 1990. [31] A. R. J. Silalahi, A. H. Boschitsch, R. C. Harris, and M. O. Fenley. Comparing the predictions of the nonlinear Poisson–Boltzmann equation and the ion size-modified Poisson–Boltzmann equation for a low-dielectric charged spherical cavity in an aqueous salt solution. J. Chem. Theory Comput., 6:3631–3639, 2010. [32] U. P. Strauss and Y. P. Leung. Volume changes as a criterion for site binding of counterions by polyelectrolytes. J. Amer. Chem. Soc., 87:1476–1480, 1965. [33] G. Tresset. Generalized Poisson–Fermi formalism for investigating size correlation effects with multiple ions. Phys. Rev. E, 78:061506, 2008. [34] M. Valisk´o, D. Boda, and D. Gillespie. Selective adsorption of ions with different diameter and valence and highly charged interfaces. J. Phys. Chem. C, 111:15575– 15585, 2007. 29
[35] V. Vlachy. Ionic effects beyond Poisson–Boltzmann theory. Annu. Rev. Phys. Chem., 50:145–165, 1999. [36] H. L. V¨ortler, K. Sch¨afer, and W. R. Smith. Simulation of chemical potentials and phase equilibria in two- and three-dimensional square-well fluids: finite size effects. J. Phys. Chem. B, 112:4656–4661, 2008. [37] J. Wen, S. Zhou, Z. Xu, and B. Li. Competitive adsorption and ordered packing of counterions near highly charged surfaces: From mean-field theory to Monte Carlo simulations. Phys. Rev. E, 85:041406, 2012. [38] S. Zhou, Z. Wang, and B. Li. Mean-field description of ionic size effects with nonuniform ionic sizes: A numerical approach. Phys. Rev. E, 84:021901, 2011.
30
∗
Pei Liu
†
Zhenli Xu
‡
Shenggao Zhou
§
July 18, 2013
Abstract Near a charged surface, counterions of different valences and sizes cluster; and their concentration profiles stratify. At a distance from such a surface larger than the Debye length, the electric field is screened by counterions. Recent studies by a variational mean-field approach that includes ionic size effects and by Monte Carlo simulations both suggest that the counterion stratification is determined by the ionic valence-tovolume ratios. Central in the mean-field approach is a free-energy functional of ionic concentrations in which the ionic size effects are included through the entropic effect of solvent molecules. The corresponding equilibrium conditions define the generalized Boltzmann distributions relating the ionic concentrations to the electrostatic potential. This paper presents a detailed analysis and numerical calculations of such a free-energy functional to understand the dependence of the ionic charge density on the electrostatic potential through the generalized Boltzmann distributions, the role of ionic valenceto-volume ratios in the counterion stratification, and the modification of Debye length due to the effect of ionic sizes.
1
Introduction
Electrostatic interactions between macromolecules and mobile ions in an aqueous solvent generate strong, long-ranged forces that play a key role in biological processes. In such interactions, ionic sizes or excluded volumes can affect many of the detailed chemical and ∗
Department of Mathematics and NSF Center for Theoretical Biological Physics, University of California, San Diego, 9500 Gilman Drive, Mail code: 0112, La Jolla, CA 92093-0112, USA. Email: [email protected]. † Department of Mathematics and Institute of Natural Sciences, Shanghai Jiao Tong University, 800 Dongchuan Rd., Shanghai, 200240, P. R. China. Email: [email protected]. ‡ Department of Mathematics, Institute of Natural Sciences, and Ministry of Education Key Laboratory in Scientific and Engineering Computing, Shanghai Jiao Tong University, 800 Dongchuan Rd., Shanghai, 200240, P. R. China. Email: [email protected]. § Department of Mathematics and NSF Center for Theoretical Biological Physics, University of California, San Diego, 9500 Gilman Drive, Mail code: 0112, La Jolla, CA 92093-0112, USA. Email: [email protected].
1
physical properties of an underlying biological system. For instance, differences in ionic sizes can affect how mobile ions bind to nucleic acids [3, 4, 32]. The size of monovalent cations can influence the stability of RNA tertiary structures [21]. The ionic size effect is more profound in the ion channel selectivity, see, e.g., [16, 24]. Detailed Monte Carlo simulations and integral equations calculations also confirm some of these experimentally observed properties due to the non-uniformity of ionic sizes [17, 28, 34, 36]. The classical Poisson–Boltzmann (PB) equation (PBE) is perhaps the most widely used and efficient model of the electrostatics in ionic solutions, particularly in biomolecular systems [2, 10, 12, 14, 27, 30]. Despite of its many successful applications, this mean-field theory is known to fail in capturing the ionic size effects and ion-ion correlations [15, 17, 28, 35, 36]. Numerical calculations based on the classical PBE often predict unphysically high concentrations of counterions near a charged surface [5, 38]. Recently, there has been a growing interest in studying the ionic size effects within the PB-like mean-field framework [5–7, 9, 11, 18–20, 22, 23, 25, 26, 31, 33, 38]. One of the main approaches is to include the additional entropic effect of solvent molecules in the free-energy functional of all local ionic concentrations, c1 = c1 (x), . . . , cM = cM (x) (M being the total number of ionic species), similar to that in the classical PB model. Given the volume vi of an ion of the ith species (1 ≤ i ≤ M ) and the volume v0 of a solvent molecule, one can define the local solvent concentration c0 = c0 (x) by " # M X vi ci (x) . (1.1) c0 (x) = v0−1 1 − i=1
The size-modified electrostatic free-energy functional is then given by [5, 20, 22, 26] # Z " M M X X 1 −1 F [c1 , · · · , cM ] = ρψ + β ci [ln (vi ci ) − 1] − µi ci dV. 2 i=0 i=1
(1.2)
P Here, ρ = M i=1 qi ci is the ionic charge density, where qi = zi e, and zi is the valence of an ion of the ithe species and e is the elementary charge, β = (kB T )−1 with kB the Boltzmann constant and T the temperature, and µi is the chemical potential of an ion of the ith species. The electrostatic potential ψ is determined by Poisson’s equation −∇ · ε∇ψ = ρ, together with some boundary conditions, where ε is the dielectric coefficient. One easily verifies that the equilibrium conditions δci F = 0 (i = 1, . . . , M ) are given by vi ln (v0 c0 ) − ln (vi ci ) = β (qi ψ − µi ) , v0
i = 1, . . . , M,
(1.3)
where ψ is the electrostatic potential corresponding to the equilibrium ionic concentrations c1 , . . . , cM . It is shown in [22] that this system of nonlinear algebraic equations has a unique solution ci = Bi (ψ) (i = 1, . . . , M ), defining the generalized Boltzmann distributions. If v0 = v1 = · · · = vM , then these distributions are [22, 23] ci =
−βqi ψ c∞ i e
1+
PM
j=1
−βqj ψ − 1) vj c∞ j (e
2
,
i = 1, . . . , M,
where c∞ i =
vi−1 eβµi , P βµj 1+ M j=1 e
i = 1, . . . , M.
Note that, for each i, c∞ i is the value of ci when ψ = 0. In general, if v0 , v1 , . . . , vM are not all the same, then explicit formulas for ci = Bi (ψ) (i = 1, . . . , M ) seem unavailable. In studying the ionic size effects in biomolecular reactions, Lu and Zhou [26] minimized the free-energy functional (1.2) by solving for a steady-state solution of a size-modified Poisson–Nernst–Planck system of partial differential equations. In [38], Zhou et al. developed a constrained optimization method for numerically minimizing the functional (1.2) with Poisson’s equation as a constraint. They found the remarkable stratification of ionic concentrations near a highly charged surface. Moreover, they found that the ionic valenceto-volume ratio is an important parameter that determines the stratification. Specifically, the counterion with the largest valence-to-volume ratio has the highest ionic concentration near the surface, the counterion with the second largest such ratio has also the second highest ionic concentration, and so on. Subsequently, Wen et al. [37] performed Monte Carlo simulations that confirm such a phenomenon. In this work, we further study the size-modified, mean-field free-energy functional (1.2). We begin with a description of two different forms of the electrostatic free-energy functional, one with and the other without the constraint of total number of ions for each ionic species. We show the equivalence of these two forms; cf. Theorem 2.2. We also derive the conditions for the equilibrium concentrations; cf. (1.3). We then focus on the equilibrium conditions (1.3) to analyze the generalized Boltzmann distributions.PIn particular, we characterize the behavior of the equilibrium ionic charge density ρi = M i=1 qi ci as a function of ψ. Our main results of analysis are as follows: (1) The concentration of ions with the largest and that with smallest valence-to-volume ratios monotonically decreases and monotonically increases, respectively, with respect to the electrostatic potential ψ. As ψ → −∞ (or ∞), the concentration of ions with the largest (or smallest) valence-to-volume ratio approaches to the inverse of its corresponding volume, while all other concentrations approach zero. (2) The total ionic charge density ρi = ρi (ψ) is a monotonic function of the potential ψ. If all ions are cations (i.e., all zj > 0), then this density is always positive, ρi (∞) = 0, and ρi (−∞) exists and is finite. If all ions are anions (i.e., all zj < 0), then this density is always negative, ρi (−∞) = 0, and ρi (∞) exists and is finite. If both cations and anions exist, then both ρi (−∞) and ρi (∞) exist and are finite, and the density reaches zero at a unique value of the potential ψ. All these are summarized in Figure 2 in Subsection 3.2. ˆ D , is derived with a weak potential limit. A formula (3) A size-modified Debye length, λ ˆ D is given in (3.23) in Subsection 3.3. of of λ Note that ionizable groups of biomolecules can dissociate in aqueous solution to produce anions or cations. So, studying an ionic system with only cations or only anions is of interest. In general, such a system can serve a model approximation system, as near a charged surface the concentration of coions is usually very low. We finally consider a charged spherical molecule of radius R0 in an ionic solution that 3
occupies a region defined by {R0 < |x| < R∞ } for some R∞ > R0 , and study the PBE, both classical and size-modified. In the case R∞ is finite, the detailed properties of the equilibrium electrostatic potential are described through two parameters: one is the prescribed surface charge density σ at the boundary |x| = R0 and the other is an effective surface charge density σ∞ at the boundary |x| = R∞ . This effective density is set by the total charge neutrality. When σσ∞ < 0 the potential ψ is monotonic. Otherwise, it decreases then increases, or it increases then decreases. These are summarized in Theorem 4.2 and illustrated in Figure 3 in Section 4. With the same setting of a charged spherical molecule, we numerically minimize the functional (1.2) and confirm our analysis; cf. Figure 4 in Section 5. If R∞ = ∞, then only in the case that both cations and anions exist and they neutralize in the bulk does the PBE (classical or size-modified) permit a unique solution, the electrostatic potential, that is radially symmetric and monotonic, and vanishes at infinity. We remark that considering a finite R∞ is of certain interest as it is known that the effect of finite size of charged system can be significant, cf. e.g., [29], and as it is practically useful in numerical computations. The rest of this paper is organized as follows: In Section 2, we describe the electrostatic free-energy functionals and derive the conditions of equilibrium concentrations. In Section 3, we present an analysis on the equilibrium conditions and on the stratification of counterions of multiple valences and sizes near a charged surface. In Section 4, we study a charged spherical molecule immersed in an ionic solution. In Section 5, we report numerical calculations to confirm our analysis. Finally, in Section 6, we draw conclusions.
2
The Free-Energy Functional and Equilibrium Conditions
We consider an ionic solution with M (M ≥ 1) ionic species, occupying a region Ω that is a bounded, open, and connected subset of R3 with a smooth boundary Γ 6= ∅. The boundary Γ is divided into two disjoint and smooth parts ΓD and ΓN (with D for Dirichlet and N for Neumann). We assume that a surface charge density is given on ΓN and the electrostatic potential is specified on ΓD . Our set up covers the case that ΓN = ∅ and ΓD = Γ and that ΓD = ∅ and ΓN = Γ. It also covers the important case of biomolecular solvation with an implicit solvent: The region Ω is the solvent region and the part of the boundary ΓN is the solute-solvent interface; cf. Figure 1. The surface charge density on ΓN effectively replaces the density of partial charges carried by solute atoms. For each i (1 ≤ i ≤ M ), we denote by ci = ci (x) the local ionic concentration of the ith PM ionic species at position x ∈ Ω. As before, ρ = i=1 qi ci is the ionic charge density. We also denote by vi (1 ≤ i ≤ M ) the volume of an ion of the ith species and by v0 the volume of a solvent molecule. The local concentration of the solvent, c0 = c0 (x), is then defined by (1.1). For convenience we denote c = (c1 , . . . , cM ). We consider minimizing the mean-field electrostatic free-energy functional F [c] =
Z
Ω
1 ρψ dV + 2
Z
M
ΓN
X 1 σψ dS + β −1 2 i=0
Z
4
ci [ln(vi ci ) − 1] dV − Ω
M Z X i=1
µi ci dV, Ω
(2.1)
Ω
+ -
-
+
+
+
-
-
+
-
+
+ +
ΓN
-
ΓD
+
+
+ +
-
+
+
-
+ -
+
-
+ +
+
+
+ +
-
+
-
+
-
-
Figure 1: A schematic view of a charged molecule immersed in an ionic solution. The molecule is represented by the inner sphere filled with symbols +.
∇ · ε∇ψ = −ρ ψ = ψ∞ with ε ∂ψ = σ ∂n
in Ω, on ΓD ,
(2.2)
on ΓN .
Here, the first two terms in (2.1) describe the electrostatic potential energy. The electrostatic potential ψ is the unique solution of the boundary-value problem of Poisson’s equation (2.2), where ε = ε(x) (x ∈ Ω) is the dielectric coefficient assumed to be Lebesgue measurable and bounded above and below by positive constants, and n denotes the unit exterior normal at the boundary Γ. Both ψ∞ : ΓD → R and σ : ΓN → R are given bounded and smooth functions, representing the electrostatic potential at ΓD and the surface charge density at ΓN , respectively. The third term in (2.1) has the form −T S with S being the entropy. The given parameters µi (i = 1, . . . , M ) in the last term in (2.1) are the chemical potentials of the ions. We define ψD : Ω → R and ψN : Ω → R by ∇ · ε∇ψD = 0 in Ω, ∇ · ε∇ψN = 0 in Ω, ψD = ψ∞ on ΓD , ψN = 0 on ΓD , and ε ∂ψD = 0 ε ∂ψN = σ on ΓN , on ΓN , ∂n ∂n respectively. For a given function f : Ω → R that is in the Sobolev space H −1 (Ω) [1, 13], let us denote by Lf : Ω → R the unique weak solution to the following boundary-value
5
problem of Poisson’s equation ∇ · ε∇(Lf ) = −f Lf = 0 ε ∂(Lf ) = 0 ∂n
in Ω, on ΓD , on ΓN .
(The uniqueness is guaranteed if ΓD has a nonzero surface measure.) Notice that L is a linear and self-adjoint operator from the space H −1 (Ω) to that consisting of H 1 (Ω)-functions vanishing on ΓD . With these the solution ψ to the boundary-value problem (2.2) can be decomposed as ψ = ψˆ + ψD + ψN , (2.3) P M where ψˆ = L i=1 qi ci . By integration by parts, we have Z
∂ψN ˆ ε ψ dS ∂n Γ Z Z Z ˆ ˆ ψ∇ · ε∇ψN dV + ε∇ψ · ∇ψN dV = ε∇ψˆ · ∇ψN dV = Ω Ω Ω ! Z Z M X qi ci ψN dV. = −∇ · ε∇ψˆ ψN dV =
σ ψˆ dS =
ΓN
Z
∂ψN ˆ ψ dS = ε ∂n ΓN
Z
Ω
Ω
i=1
Therefore the free-energy functional F [c] can now be rewritten as ! ! Z M M M Z X X 1 X 1 F [c] = qj cj dV + qi c i L ψD + ψN − µi ci dV qi 2 Ω 2 Ω j=1 i=1 i=1 Z Z M X 1 σ (ψN + ψD ) dS. + β −1 ci [ln(vi ci ) − 1] dV + ΓN 2 i=0 Ω By formal and routine calculations using the decomposition (2.3) and the fact that L is self-adjoint, we obtain the first variations [8, 22, 23, 26] 1 vi −1 δci F [c] = qi ψ − ψD − µi + β ln(vi ci ) − ln (v0 c0 ) , i = 1, . . . , M, (2.4) 2 v0 P where ψ is determined by (2.2) with ρ = M i=1 qi ci . The equilibrium conditions δci F [c] = 0 (i = 1, . . . , M ) are then given by vi 1 ln(v0 c0 ) − ln(vi ci ) = β qi ψ − ψD − µi , i = 1, . . . , M. (2.5) v0 2 The following result can be proved by the argument used in [22]: 6
Theorem 2.1. (1) The functional F [c] is convex. Moreover, it has a unique minimizer in PM the set of all concentrations c1 , . . . , cM such that ci ≥ 0 (i = 1, . . . , M ) and i=1 vi ci ≤ 1 in Ω. (2) There exist θ1 , θ2 ∈ R with 0 < θ1 < θ2 < 1 that bound uniformly from below and above, respectively, the free-energy minimizing ionic concentrations and the corresponding concentration of solvent molecules. (3) The minimizer is characterized by the equilibrium conditions (2.5). We now consider a different formulation: minimize the electrostatic free-energy functional Z M Z X 1 1 −1 Fˆ [c] = ρψ dV + σψ dS + β ci [ln(vi ci ) − 1] dV, Ω 2 ΓN 2 Ω i=0 ∇ · ε∇ψ = −ρ in Ω, ψ = ψ∞ on ΓD , with ε ∂ψ = σ on ΓN , Z ∂n subject to ci dV = Ni , i = 1, . . . , M, Z
(2.6)
(2.7)
(2.8)
Ω
P where again ρ = M i=1 qi ci . Note that the functional F [c] defined in (2.1) has the term of the chemical potentials µi . Here the functional Fˆ [c] does not have that term. On the other hand, the minimization of Fˆ [c] is subject to the conservation of total number of ions in each of the M ionic species, defined in (2.8) in which Ni (i = 1, . . . , M ) are given positive numbers. The system (2.7) is the same as (2.2). It again follows from standard arguments that the functional Fˆ [c] with (2.7) has a unique minimizer cˆ = (ˆ c1 , . . . , cˆM ) among all c = (c1 , . . . , cM ) subject to (2.8). By a similar argument used in deriving (2.4), we have for any c = (c1 , . . . , cM ) that vi 1 −1 ˆ ln(vi ci ) − ln (v0 c0 ) , i = 1, . . . , M. (2.9) δci F [c] = qi ψ − ψD + β 2 v0 For the minimizer cˆ, we have Z δci Fˆ [ˆ c] di dV = 0 Ω
∀di :
Z
di dV = 0,
i = 1, . . . , M.
Ω
By Lemma 2.1 below, these lead to the existence of constants, still denoted by µi (i = 1, . . . , M ), such that δci Fˆ [ˆ c] = µi (i = 1, . . . , M ). By (2.9), this is exactly (2.5). We have in fact proved the following: Theorem 2.2. The variational problem of minimizing the functional (2.1) with (2.2) is equivalent to that of minimizing the functional Fˆ [c] with (2.7) and (2.8).
7
We remark that, if one knows the chemical potentials µi (i = 1, . . . , M ), then the first formulation, i.e., minimizing F [c], can be mathematically and computationally simpler than the second one that has the constraint of total numbers of ions. The question is how the chemical potentials can be estimated by experiment or other models. A similar practical issue for the second formulation with the constraint (2.8) is that one needs to have a good estimate of each Ni , the total number of ions of the ith species in the system. Lemma 2.1. Let u ∈ L2 (Ω) be such that Z Z 2 u v dV = 0 ∀v ∈ L (Ω) : v dV = 0. Ω
Ω
Then there exists a constant µ ∈ R such that u = µ almost everywhere in Ω. Proof. Denote for any w ∈ L2 (Ω) Z Z 1 w = − w dV = w dV, |Ω| Ω Ω where |Ω| is the volume (i.e., the Lebesgue measure) of Ω. Note that the integral of w − w over Ω vanishes. For any v ∈ L2 (Ω), we have by the assumption of the lemma that Z Z Z (u − u) v dV = (u − u) (v − v) dV = u (v − v) dV = 0. Ω
Ω
Ω
Therefore u = u almost everywhere in Ω. Set µ = u. The Lemma is proved.
3
Q.E.D.
Equilibrium Ionic Concentrations and Modified Debye Length
We consider the equilibrium conditions (2.5) that determine the equilibrium ionic concentrations. Introducing a shifted potential φ = ψ − ψD /2, we can rewrite (2.5) as vi ln(v0 c0 ) − ln(vi ci ) = β (qi φ − µi ) , v0
i = 1, . . . , M.
(3.1)
As c0 is given in (1.1), this is a system of nonlinear algebraic equations for the unknowns c1 , . . . , cM . In [22], Li proved the following: Theorem 3.1. (1) For each φ ∈ R the system (3.1) has a unique solution ci = Bi (φ) ∈ (0, vi−1 ) (i = 1, . . . , M ). Moreover each Bi : R → (0, vi−1 ) is a smooth function. (2) Define V : R → R by M Z φ X V (φ) = − qi Bi (ξ) dξ + V0 , (3.2) i=1
0
where V0 is a constant. Then the function V : R → R satisfies V ′′ > 0 on R and is thus strictly convex. 8
The relations ci = Bi (φ) (i = 1, . . . , M ) define the generalized Boltzmann distributions for the equilibrium concentrations c1 , . . . , cM . In general, PM explicit formulas of such distribu′ tions seem to be unavailable. Notice that −V (φ) = i=1 qi Bi (φ) with φ = ψ − ψD /2 is the ionic charge density. The proof of V ′′ > 0 in R given in [22] (cf. the proof of Lemma 5.2 in [22]) applies to our present case; cf. the remark in Subsection 3.1. The electrostatic potential ψ is now the unique solution to the boundary-value problem of the generalized and implicit PBE ψD ′ ∇ · ε∇ψ − V ψ − = 0. (3.3) 2 The boundary conditions are the same as in (2.2). The existence and uniqueness of the equilibrium concentrations ci (i = 1, . . . , M ) lead to that of the solution to this boundaryvalue problem.
3.1
Equilibrium ionic concentrations
For convenience, we denote z0 = 0, q0 = 0, and µ0 = 0. We assume that all the M ionic species have different ionic valence-to-volume ratios. We assume there are l ≥ 0 species of anions and m ≥ 0 species of cations. So M = l+m. If l = 0 then all ions are cations. If m = 0 then all ions are anions. We label all the M ionic species and the solvent concentration by −l, . . . , −1, 0, 1, . . . , m, and assume that z−1 z0 z1 zm z−l < ··· < < =0< < ··· < . (3.4) v−l v−1 v0 v1 vm For convenience, we denote J = {−l, . . . , −1, 0, 1, . . . , m}. With our notation, the equilibrium conditions (2.5) then become vi ln(v0 c0 ) − ln(vi ci ) = β (qi φ − µi ) ∀i ∈ J. (3.5) v0 Note that by Theorem 3.1 ci = Bi (φ) for any i ∈ J with i 6= 0. We thus have by (1.1) that " # X c0 = B0 (φ) = v0−1 1 − vi Bi (φ) . (3.6) i∈J,i6=0
Theorem 3.2. Let ci = Bi (φ) (i ∈ J) be the equilibrium concentrations defined by (3.5) and (3.6). Assume (3.4) holds true. (1) Monotonicity of concentrations with the largest and smallest valence-to-volume ratios. We have c′−l (φ) > 0 and c′m (φ) < 0 for all φ ∈ R. (2) Stratification of concentration profiles. If i ∈ {−l, . . . , m − 1} then (vi ci )1/vi < (vi+1 ci+1 )1/vi+1
if and only if
φ
µi−1 /vi−1 − µi /vi . qi−1 /vi−1 − qi /vi
If i ∈ {−l + 1, . . . , m} then (vi ci )1/vi > (vi−1 ci−1 )1/vi−1
9
(3) Asymptotic behavior of concentrations. We have v−l c−l → 1 and ci → 0 (i = −l + 1, . . . , m) as φ → ∞, vm cm → 1 and ci → 0 (i = −l, . . . , m − 1) as φ → −∞.
(3.7) (3.8)
Proof. We first prove that (vi ci )1/vi = e−β(µj /vj −µi /vi ) eβ(qj /vj −qi /vi )φ (vj cj )1/vj
∀i, j ∈ J.
(3.9)
In fact, it follows immediately from (3.5) that (v0 c0 )1/v0 = e−βµi /vi eβ(qi /vi )φ (vi ci )1/vi
∀i ∈ J.
(3.10)
(vi ci )1/vi = eβµi /vi e−β(qi /vi )φ (v0 c0 )1/v0
∀i ∈ J.
(3.11)
This also implies It follows from (3.10) and (3.11) that (3.9) holds true if one of i and j is 0. Now let i, j ∈ J with i 6= j, i 6= 0, and j 6= 0. We have by (3.5) that 1 βqi βµi 1 βqj βµj ln (vi ci ) + φ− = ln (vj cj ) + φ− . vi vi vi vj vj vj Consequently, (vi ci )1/vi = e−β(µj /vj −µi /vi ) eβ(qj /vj −qi /vi )φ (vj cj )1/vj . This, together with (3.10) and (3.11), implies (3.9). We now prove parts (1)–(3). (1) Denote γ = v0 c0 and ηi = eβµi for each i ∈ J. It follows from the equilibrium conditions (3.5) that ηi ci = γ vi /v0 e−βqi φ ∀i ∈ J. (3.12) vi This and the definition of c0 (cf. (1.1)) imply X ηi γ vi /v0 e−βqi φ = 1.
(3.13)
i∈J
Taking the derivative of both sides of (3.13) with respective to φ, and using our notation z0 = 0 and µ0 = 0, we obtain P βv0 γ i∈J ηi qi γ vi /v0 e−βqi φ ′ γ (φ) = P . vi /v0 e−βqi φ i∈J ηi vi γ
Fix k ∈ J. Taking the derivative of both sides of (3.12) with i = k with respect to φ and using the above expression of γ ′ (φ), we get P ηi vi (qi /vi − qk /vk ) γ vi /v0 e−βqi φ c′k (φ) = βηk γ vk /v0 e−βqk φ i∈J P . (3.14) vi /v0 e−βqi φ i∈J ηi vi γ 10
This, the assumption (3.4), and the definition qi = zi e (i ∈ J) imply that c′m (φ) < 0 and c′−l (φ) > 0. (2) This follows directly from (3.9). (3) Since 0 ≤ v−l c−l ≤ 1, we have again by (3.9) with j = −l and (3.4) that for any i ∈ J with i > −l (vi ci )1/vi ≤ e−β(µ−l /v−l −µi /vi ) eβ(q−l /v−l −qi /vi )φ → 0
as φ → ∞.
P
Since i∈J vi ci = 1, we thus also have v−l c−l → 1 as φ → ∞. This proves (3.7). The proof of (3.8) is similar. Q.E.D. Remark. By (3.14) and the Cauchy–Schwarz inequality, we obtain through a series of calculations that h P 2 P i P 2 2 β − X i∈J qi vi ci i∈J qi ci i∈J vi ci P qk c′k (φ) = ≤ 0. 2 i∈J vi ci k∈J Since qi /vi (i ∈ J) are not all the same, this is a strict inequality. Thus the function V : R → R, defined in (3.2) with −V ′ (φ) being the total ionic charge density, is a strictly convex function of φ; cf. Theorem 3.1.
3.2
Equilibrium charge density
We now study the properties of the ionic charge density V ′ (φ) = − the function V is given now by XZ φ V (φ) = − qi Bi (ξ) dξ + V0 . i∈J
P
i∈J
qi Bi (φ). By (3.2),
(3.15)
0
We distinguish three cases: (1) l = 0, m = M , and J = {0, 1, . . . , M }; (2) l = M , m = 0, and J = {−M, . . . , −1, 0}; (3) l ≥ 1, m ≥ 1, l + m = M , and J = {−l, . . . , −1, 0, 1, . . . , m}. We denote f (−∞) = lims→−∞ f (s) and f (∞) = lims→∞ f (s) if the limits exist. Theorem 3.3. Assume (3.4) holds true. (1) If l = 0 and m = M , and J = {0, 1, . . . , M }, then V ′ (−∞) = −qM /vM < 0, V ′ (∞) = 0, and V ′ (φ) < 0 for all φ ∈ R. Moreover, V (−∞) = ∞, and V (∞) = 0 with a suitably chosen constant V0 in (3.2). (2) If l = M , m = 0, and J = {−M, . . . , −1, 0}, then V ′ (−∞) = 0, V ′ (∞) = −q−M /v−M > 0, and V ′ (φ) > 0 for all φ ∈ R. Moreover, V (−∞) = 0 with a suitably chosen constant V0 in (3.2), and V (∞) = ∞. (3) If l ≥ 1, m ≥ 1, l + m = M , and J = {−l, . . . , −1, 0, 1, . . . , m}, then V ′ (−∞) = −qm /vm < 0, V ′ (∞) = −q−l /v−l > 0, and −qm /vm < V ′ (φ) < −q−l /v−l for all φ ∈ R. Moreover, minφ∈R V (φ) = 0 with a suitably chosen constant V0 in (3.2), V (−∞) = ∞, and V (∞) = ∞. 11
In Figure 2, we show the schematic behavior of the convex function V for the three different cases. V
V′
0
φ
0
φ
−qM /vM
Case 1. 0 = z0 /v0 < z1 /v1 < · · · < zM /vM . V
V′ −q−M /v−M
0
φ
0
φ
Case 2. z−M /v−M < · · · < z−1 /v−1 < z0 /v0 = 0. V
V′ −q−l /v−l
0
φ0
φ
0
φ0
φ
−qm /vm
Case 3. z−l /v−l < · · · < z−1 /v−1 < z0 /v0 = 0 < z1 /v1 < · · · < zm /vm . Figure 2: Graphs of V ′ = V ′ (φ) (left) and V = V (φ) (right) for the three different cases. The value −V ′ (φ) is the ionic charge density at a (shifted) electrostatic potential φ.
12
Proof of Theorem 3.3. (1) In this case, we have V ′ (−∞) = −qM /vM < 0 by (3.8) and V ′ (∞) = 0 by (3.7). Since V ′′ > 0 in R by Theorem 3.1, V ′ is a strictly increasing function. Therefore, V ′ (−∞) < 0 and V ′ (∞) = 0 imply that V ′ (φ) < 0 for all φ ∈ R. Fix φ0 ∈ R so that V ′ (φ0 ) < 0. Then for any φ < φ0 we have V (φ) = V (φ0 ) + V ′ (ξ)(φ − φ0 ), where ξ ∈ (φ, φ0 ). Notice that V ′ (ξ) < V ′ (φ0 ) < 0. So, V (φ) → ∞ as φ → −∞, i.e., V (−∞) = ∞. Finally, setting j = 0 in (3.9), we obtain that 1/vi 1/vi βµi /vi (vi Bi (φ)) = (vi ci ) ≤ max e e−β(qi /vi )φ . 1≤i≤M
P Therefore, the integral of M i=1 qi Bi (φ) from 0 to φ > 0 is a bounded function of φ. Thus, by the definition of V (cf. (3.2)), V (φ) is bounded below. Since V ′ < 0 on R, V decreases. Therefore, V (∞) exists. Choosing the constant XZ ∞ V0 = qi Bi (φ) dφ i∈J
0
in (3.15), we then have V (∞) = 0. (2) This part can be proved similarly. (3) The fact that V ′ (−∞) = −qm /vm < 0 and V ′ (∞) = −q−l /v−l > 0 can be proved by the same argument used for proving part (1). Since V ′′ > 0 in R, V ′ (φ) is a strictly increasing function. Therefore, −qm /vm < V ′ (φ) < −ql /vl for all φ ∈ R. Again by the same argument used in part (1), we have V (−∞) = ∞ and V (∞) = ∞. It is clear now that there exists a unique φ0 ∈ R such that V ′ (φ0 ) = 0 and V (φ0 ) = minφ∈R V (φ). Setting in (3.15), V0 =
XZ i∈J
we have minφ∈R V (φ) = 0.
3.3
φ0
qi Bi (φ) dφ, 0
Q.E.D.
Modified Debye length
We now linearize V ′ (φ) at φ = 0 to obtain a modified Debye length. We first define c∞ i = Bi (0) for any i ∈ J. Note that X v0 c∞ = 1 − vi c∞ (3.16) 0 i . i∈J,i6=0
By (3.10) with φ = 0 and ηi = eβµi , we also have that vi /v0 (v0 c∞ = ηi−1 vi c∞ 0 ) i ,
i ∈ J.
We denote again γ = v0 c0 which depends on φ. We have for |φ| ≪ 1 that γ = τ0 + τ1 φ + O φ 2 , 13
(3.17)
(3.18)
where τ0 and τ1 are two constants to be determined. Using our notation, (1.1) becomes X v0 c0 = 1 − vi ci . (3.19) i∈J,i6=0
This and (3.16) imply that τ0 = v0 c0 (0) = 1 −
X
∞ vi c∞ i = v0 c0 .
(3.20)
i∈J,i6=0
Substituting (3.18) into (3.12) for a fixed i ∈ J with i 6= 0, we obtain by (3.20) and (3.17) that v /v ηi ci = τ0 + τ1 φ + O φ2 i 0 1 − βqi φ + O φ2 vi vi /v0 τ1 ηi vi /v0 2 1+ φ+O φ 1 − βqi φ + O φ2 = τ0 vi τ0 v i τ1 ∞ 2 = ci 1 + φ+O φ 1 − βqi φ + O φ2 vτ 0 0 v i τ1 ∞ ∞ = ci + ci − βqi φ + O φ2 . (3.21) v 0 τ0
It now follows from (3.18), (3.19), (3.21), and (3.20) that X X v i τ1 ∞ ∞ 2 − βqi φ + O φ2 . vi ci − vi ci τ0 + τ1 φ + O φ = 1 − 2 ∞ v0 c0 i∈J,i6=0 i∈J,i6=0
Comparing the leading-order terms, we obtain exactly (3.20). Comparing the O(φ)-terms and noting that z0 = 0, we obtain P βv02 c∞ qi v i c ∞ 0 P i∈J2 ∞ i . τ1 = i∈J vi ci
This and (3.21) then lead to ! P ∞ v q v c i j j ∞ P j∈J 2 ∞ j − qi φ + O φ2 ci = c∞ i + βci j∈J vj cj
∀i ∈ J with i 6= 0.
Consequently,
V ′ (φ) = −
X i∈J
qi c i = −
X i∈J
" P # ∞ 2 X q v c i i i i∈J 2 ∞ 2 P − q c qi c ∞ − β . φ + O φ i i i 2 ∞ i∈J vi ci i∈J
In the small |φ| approximation, we then obtain by (3.3) and the shift φ = ψ − ψD /2 the linearized PBE with ionic size effects X 1 ˆ −2 ˆ −2 ψ = − (3.22) ∇ · ε∇ψ − ελ qi c ∞ i − ελD ψ D , D 2 i∈J 14
ˆ D > 0 is defined by where λ ˆ −2 λ D
" # P ∞ 2 q v c β X 2 ∞ i i i . = q c − Pi∈J 2 ∞ ε i∈J i i i∈J vi ci
(3.23)
The fact that the right-hand side of the above equation is nonnegative follows from the Cauchy–Schwarz inequality ! ! !2 X X X qi2 c∞ vi2 c∞ ≤ . qi v i c ∞ i i i i∈J
i∈J
i∈J
In the case that ε is a constant, ψD = 0 on ΓD , and i∈J qi c∞ i = 0 which is the charge neutrality, the linearized equation (3.22) simplifies to the familiar equation P
ˆ −2 ψ = 0. ∆ψ − λ D ˆ D the size-modified Debye length. Recall that the Debye length λD > 0 in the We call λ classical Debye–H¨ uckel theory is defined by λ−2 D =
βX 2 ∞ q c . ε i∈J i i
(3.24)
ˆ D > λD . This indicates that finite ionic sizes lead to a longer Debye It is clear that λ screening length.
4
A Charged Spherical Molecule: Basic Properties
In this section, we consider a charged spherical molecule, centered at the origin and with radius R0 > 0, immersed in a solvent that occupies the region defined by R0 < r < R∞ in the spherical coordinates for some R∞ > R0 , where R∞ can be finite or infinite. Corresponding to the general setting described in the previous section, we have in this case that Ω = {x ∈ R3 : R0 < |x| < R∞ }, ΓN = {x ∈ R3 : |x| = R0 }, ( {x ∈ R3 : |x| = R∞ } if R∞ < ∞, ΓD = {∞} if R∞ = ∞. We consider the boundary-value problem of the PBE, classical and size-modified, for the electrostatic potential ψ (cf. (2.2) and (3.3)): ε∆ψ = V ′ (ψ) in Ω, ∂ψ (4.1) =σ on ΓN , ε ∂n ψ=0 on ΓD . 15
If ψ = ψ(r) (r = |x|) is radially symmetric, then this system is the same as 2 ′′ ′ ε ψ + ψ = V ′ (ψ) if R0 < r < R∞ , r σ ′ ψ (R ) = − , 0 ε ψ(R∞ ) = 0.
(4.2)
Here, we assume both the dielectric coefficient ε > 0 and surface charge density σ are constants. We can consider the inhomogeneous Dirichlet boundary condition ψ = ψ∞ on ΓD for a constant ψ∞ by replacing ψ by ψ − ψ∞ . We assume that V : R → R is smooth, V ′′ > 0 in R, and inf s∈R V (s) = 0. We then have one and only one of the three cases covered in Theorem 3.3: V ′ > 0 in R, V ′ < 0 in R, or there exists a unique φ0 such that V (s) > V (φ0 ) = 0 for any s 6= 0. In the last case, V ′ < 0 in (−∞, φ0 ), V ′ (φ0 ) = 0, and V ′ > 0 in (φ0 , ∞). The term −V ′ (ψ) is the ionic charge density that is related to M ionic species with valences and volumes satisfying P (3.4). For the classical PBE, we have V (φ) = i∈J c∞ e−βqi φ − 1 , where c∞ ∈ J, i 6= 0) are i i (iP bulk ionic concentrations. Under the assumption on the charge neutrality i∈J qi c∞ i = 0, this V satisfies the properties in the third case in Theorem 3.3 with φ0 = 0. So our analysis covers this case. We first consider the case that R∞ < ∞. Theorem 4.1. Let ε > 0, σ ∈ R, and 0 < R0 < R∞ < ∞. Assume that V : R → R is smooth, V ′′ > 0 in R, and inf s∈R V (s) = 0. Then the boundary-value problem (4.1) has a unique weak solution. Moreover it is radially symmetric and smooth on Ω. Proof. We assume that σ 6= 0. The case that σ = 0 can be proved similarly. We denote H = {φ ∈ H 1 (R0 , R∞ ) : φ(R∞ ) = 0} and define I : H → (−∞, ∞] by Z R∞ h i ε ′ 2 I[φ] = 4π φ (r) + V (φ) r2 dr − 4πR02 σφ(R0 ) ∀φ ∈ H. (4.3) 2 R0 If φ ∈ H, then we have with δ = ε/(4σ) > 0 that Z |φ(R0 )| = leading to
R∞
Z p φ (r) dr ≤ R∞ − R0
R∞ ′
′
R0
I[φ] ≥
πεR02
Z
|φ (r) | dr R0
R∞
|φ′ (r) |2 dr − R0
2
1/2
R∞ − R0 ≤ +δ 4δ
4πσ 2 R02 (R∞ − R0 ) ε
Z
R∞
|φ′ (r) |2 dr, R0
∀φ ∈ H.
Notice that by Poincar´e’s inequality, I[φ] controls kφk2H 1 (R0 ,R∞ ) , up to an additive constant. By standard arguments of direct methods in the calculus of variations and by the fact that H 1 (R0 , R∞ ) can be compactly embedded into C([R0 , R∞ ]), we obtain the existence of a minimizer ψ ∈ H of the functional I : H → [0, ∞]. Moreover, the minimum value of I 16
over H is finite. The uniqueness of this minimizer follows from the strict convexity of I. Clearly, the minimizer ψ = ψ(r) is a smooth solution to the corresponding Euler–Lagrange equation, the first equation in (4.2). Moreover, it satisfies the natural boundary condition ψ ′ (R0 ) = −σ/ε. Since ψ ∈ H, ψ(R∞ ) = 0. If we also denote ψ(x) = ψ(|x|), then it is clear that this function is a radially symmetric solution to (4.1), and is smooth on Ω. To prove the uniqueness of solution to (4.1), we first note that the solution ψ minimizes the functional Z h Z i ε 2 Q[φ] = |∇φ| + V (φ) dV − σφ dS Ω 2 ΓN over all φ ∈ H 1 (Ω) such that φ = 0 on ΓD . Indeed, for any η ∈ H 1 (Ω) with η = 0 on ΓD , we have by integration by parts, the convexity of V , and (4.1) that Z Z h i ε 2 |∇η| + ε∇ψ · ∇η + V (ψ + η) − V (ψ) dV − ση dS Q[ψ + η] − Q[ψ] = ΓN Ω 2 Z Z ′ ≥ [ε∇ψ · ∇η + V (ψ)η] dV − ση dS Ω ΓN Z Z Z ∂ψ ′ = [−ε∆ψ + V (ψ)] η dV + ε η dS − ση dS Ω ΓN ∂n ΓN = 0.
(4.4)
Now if ψˆ is another solution to (4.1), which is not assumed to be radially symmetric, then both ψ and ψˆ are minimizers of the functional Q. Since Q is strictly convex, it has at most one minimizer. Therefore, ψˆ = ψ. Q.E.D. We now study the behavior of the equilibrium potential ψ = ψ(|x|) that is the solution to the boundary-value problem (4.1). We denote by T (r) the total amount of ionic charges in the region {x ∈ Ω : R0 < |x| < r}: Z r T (r) = −4π V ′ (ψ(s))s2 ds ∀r ∈ [R0 , R∞ ]. (4.5) R0
Then T∞ := T (R∞ ) represents the total amount of ionic charges in Ω. If we integrate both sides of Poisson’s equation in (4.1), apply integration by parts, and use the boundary condition at r = R0 , then we have 2 T∞ + 4πR02 σ + 4πR∞ εψ ′ (R∞ ) = 0.
(4.6)
This is the charge neutrality of the system, as σ∞ := εψ ′ (R∞ ) can be viewed as the effective surface charge density on the part of the boundary ΓD = {x ∈ R3 : |x| = R∞ } which is set due to the finiteness of the underlying domain. The surface charge density σ in the boundary condition and the effective surface charge density σ∞ = εψ ′ (R∞ ) are the key parameters that determine the monotonicity of the electrostatic potential ψ. Theorem 4.2. Let ε > 0, σ ∈ R, and 0 < R0 < R∞ < ∞. Assume that V : R → R is smooth, V ′′ > 0 in R, and inf s∈R V (s) = 0. Consider V ′ > 0 in R, or V ′ < 0 in R, or V ′ (φ0 ) = 0 for a unique φ0 ∈ R. Let ψ = ψ(|x|) be the solution to the boundary-value problem (4.1). 17
(1) If σ > 0 and σ∞ > 0, then there exists a unique Rc+ ∈ (R0 , R∞ ) such that ψ ′ < 0 in (R0 , Rc+ ) and ψ ′ > 0 in (Rc+ , R∞ ). If σ < 0 and σ∞ < 0, then there exists a unique Rc− ∈ (R0 , R∞ ) such that ψ ′ > 0 in (R0 , Rc− ) and ψ ′ < 0 in (Rc− , R∞ ). (2) If σσ∞ ≤ 0 and at least one of σ and σ∞ is nonzero, then ψ ′ > 0 in (R0 , R∞ ) if σ < 0 or σ∞ > 0, and ψ ′ < 0 in (R0 , R∞ ) if σ > 0 or σ∞ < 0. (3) The case that σ = σ∞ = 0 can only occur when V ′ (0) = 0. In this case ψ = 0 identically in (R0 , R∞ ).
ψ
R0
ψ
R∞
r
R0
(a) σ > 0 and σ∞ < 0.
r
(b) σ > 0 and σ∞ > 0.
ψ R0
R∞
ψ R∞
r
R0
(c) σ < 0 and σ∞ > 0.
R∞
r
(d) σ < 0 and σ∞ < 0.
Figure 3: The graph of electrostatic potential ψ.
The main results of Theorem 4.2 are illustrated in Figure 3. To prove the theorem, we first prove a lemma. Lemma 4.1. Given the assumption as in Theorem 4.2 and A := {r ∈ (R0 , R∞ ) : ψ ′ (r) = 0}. (1) If V ′ > 0 in R or V ′ < 0 in R, then either A = ∅ or A contains exactly one point. (2) If V ′ (φ0 ) = 0 at some unique φ0 ∈ R, then either A = ∅, or A contains exactly one point, or A = (R0 , R∞ ) in which case we have additionally that φ0 = 0 and ψ = 0 identically in (R0 , R∞ ). 18
Proof. Suppose there exist r1 , r2 ∈ R such that R0 ≤ r1 < r2 ≤ R∞ and ψ ′ (r1 ) = ψ ′ (r2 ) = 0. Denote ω = {x ∈ R3 : r1 < |x| < r2 }. Then, by the same argument used in proving Theorem 4.1 (cf. (4.4)) and the strict convexity of V , the restriction of ψ onto ω is the unique minimizer among all φ ∈ H 1 (ω) of the functional Z h i ε G[φ] = |∇φ|2 + V (φ) dV. (4.7) ω 2 (1) Suppose V ′ > 0. Since inf s∈R V (s) = 0, we have V (−∞) = 0. Thus, with φk = −k (k = 1, 2, . . . ), we have 0 ≤ G[ψ] ≤ G[φk ] → 0 as k → ∞. Thus G[ψ] = 0. But this is impossible: If the integral of |∇ψ|2 vanishes, then ψ is a constant. But V (s) > 0 for any s ∈ R. Thus the integral of V (ψ) over ω is positive, a contradiction to G[ψ] = 0. For this case and similarly the case that V ′ < 0, the set A can then contain at most one point. (2) In this case, V (φ0 ) = mins∈R V (s) = 0. Thus, the constant function φ = φ0 also minimizes G over H 1 (ω). By the uniqueness of minimizer which follows from the strict convexity of G, ψ = φ0 in [r1 , r2 ] and [r1 , r2 ] ⊆ A. Let s = inf A ≥ R0 and t = sup A ≤ R∞ . At any point in (s, t), ψ = φ0 and ψ ′ vanishes. Moreover, ψ ′ 6= 0 in (R0 , s) if s > R0 and in (t, R∞ ) if t < R∞ . We prove now s = R0 and t = R∞ . Assume s > R0 and ψ ′ > 0 in (R0 , s). Differentiating both sides of ∆ψ = V ′ (ψ)/ε in R0 < r < t, we obtain 2 ′′ 2 1 ′′ ′′′ ψ (r) + ψ (r) = V (ψ(r)) + 2 ψ ′ (r) ∀r ∈ (R0 , t). (4.8) r ε r Hence, setting u = ψ ′ , we have Lu := ∆u + c(r)u = 0 for R0 < r < t, where c(r) := − [(1/ε)V ′′ (ψ(r)) + 2/r2 ] < 0 and c(r) is bounded in (R0 , s). Clearly u(s) = 0 ≤ u(r) for all r ∈ (R0 , t). Hence u attains its minimum with the minimum value 0 at any interior point x with |x| = s in {x ∈ R3 : R0 < |x| < t}. By the Strong Maximum Principle (cf. Theorem 3.5 of [13]), u must be a constant in R0 < r < t. This contradicts the fact that u > 0 in R0 < r < s and u = 0 in s < r < t. Hence ψ ′ > 0 in (R0 , s) is impossible. Similarly, ψ ′ < 0 in (R0 , s) is impossible. Therefore s = R0 . Similarly, we have t = R∞ . Finally A = (s, t) = (R0 , R∞ ), and since ψ(R∞ ) = 0, ψ = φ0 = 0 in (R0 , R∞ ). Q.E.D. Proof of Theorem 4.2. (1) Assume σ > 0 and σ∞ > 0. Then ψ ′ (R0 ) < 0 and ψ ′ (R∞ ) > 0. Hence by Lemma 4.1 there exists a unique Rc+ ∈ (R0 , R∞ ) such that ψ ′ (Rc+ ) = 0. Clearly, ψ ′ < 0 in (R0 , Rc+ ) and ψ ′ > 0 in (Rc+ , R∞ ). The case that σ < 0 and σ∞ < 0 can be proved similarly. (2) Assume σ > 0 and σ∞ ≤ 0, i.e., ψ ′ (R0 ) < 0 and ψ ′ (R∞ ) ≤ 0. The other cases can be treated similarly. It suffices to show that ψ ′ does not vanish in (R0 , R∞ ). Assume there existed r0 ∈ (R0 , R∞ ) such that ψ ′ (r0 ) = 0. By Lemma 4.1, ψ ′ does not vanish at other points. Thus there are two cases. Case 1: ψ ′ < 0 in (R0 , r0 ) and ψ ′ < 0 in (r0 , R∞ ). Case 2: ψ ′ < 0 in (R0 , r0 ) and ψ ′ > 0 in (r0 , R∞ ). We show that each case is impossible. Consider Case 1. If V ′ > 0 in R, then ∆ψ > 0 on {x ∈ R3 : r0 < |x| < R∞ }. Let x0 ∈ R3 be such that |x0 | = r0 . Since ψ(|x0 |) > ψ(|x|) for any x with |x| ∈ (r0 , R∞ ). Lemma 3.4 in [13] then implies that ψ ′ (r0 ) 6= 0, a contradiction. If V ′ < 0, we obtain the same contradiction ψ ′ (r0 ) 6= 0 by considering (R0 , r0 ). Finally, assume V ′ (φ0 ) = 0. If 19
ψ(r0 ) 6= φ0 , then since ψ ′ (r0 ) = 0, ψ ′′ (r0 ) = V ′ (φ(r0 ))/ε 6= 0. These imply that ψ ′ changes its sign at r0 which is a contradiction in Case 1. Thus Case 1 is impossible. Consider Case 2. In this case, we must have ψ ′ (R∞ ) = 0, since ψ ′ (R∞ ) > 0 would lead to ψ ′ (r1 ) = 0 for some r1 ∈ (r0 , R∞ ) which is impossible. Moreover, ψ reaches its minimum at r0 and in particular ψ ′′ (r0 ) > 0. If V ′ < 0 then ∆ψ < 0 in Ω. The Strong Maximum Principle then implies that ψ must be a constant for r ∈ (R0 , R∞ ). This is impossible as ψ ′ (R0 ) < 0. If V ′ > 0, then ∆ψ > 0 and ψ(R∞ ) = 0 > ψ(r) for r ∈ (r0 , R∞ ). By Lemma 3.4 in [13], we must have ψ ′ (R∞ ) > 0, leading to a contradiction. Assume V ′ (φ0 ) = 0 at some unique φ0 . We cannot have φ0 > ψ(r0 ) as this would imply that V ′ (ψ(r0 )) < 0 and hence ψ ′′ (r0 ) < 0 which contradicts ψ ′′ (r0 ) > 0. We cannot have φ0 ≤ ψ(r0 ) neither, since otherwise we would have that V ′ (ψ) ≥ 0 and hence ∆ψ ≥ 0 on the region r0 < r < R∞ . But ψ(R∞ ) = 0 > ψ(r) for any r ∈ (r0 , R∞ ). Lemma 3.4 in [13] then implies ψ ′ (R∞ ) > 0, a contradiction. Thus Case 2 is also impossible. (3) We have now ψ ′ (R0 ) = 0 and ψ ′ (R∞ ) = 0. As in the proof of Lemma 4.1 with r1 = R0 and r2 = R∞ , ψ is a constant in [R0 , R∞ ]. Since ψ(R∞ ) = 0, ψ = 0 identically in [R0 , R∞ ]. Thus V ′ (0) = V ′ (ψ) = ε∆ψ = 0. Hence φ0 = 0. Q.E.D. We recall that for any r ∈ (R0 , R∞ ), T (r) is the total amount of ionic charges in the region {x ∈ Ω : R0 < |x| < r}; cf. (4.5). In the following corollary, we consider a special and useful case but deliver more results: Corollary 4.1. Let ε > 0, σ ∈ R, and 0 < R0 < R∞ < ∞. Assume that V : R → R is smooth, V ′′ > 0 in R, and V (s) > V (0) = 0 for any s ∈ R with s 6= 0. Let ψ = ψ(|x|) be the solution to the boundary-value problem (4.1). (1) If σ = 0, then ψ = 0 and T = 0 identically. (2) If σ > 0, then ψ ′ < 0, ψ ′′ > 0, and T ′ < 0 in (R0 , R∞ ). (3) If σ < 0, then ψ ′ > 0, ψ ′′ < 0, and T ′ > 0 in (R0 , R∞ ). Proof. (1) It is easy to verify that in this case ψ = 0 is the unique solution to (4.1). Moreover, T ′ (r) = −V ′ (ψ(r)) = −V ′ (0) = 0 and T (R0 ) = 0. Hence T = 0 in (R0 , R∞ ). (2) Assume σ > 0. We then have σ∞ = εψ ′ (R∞ ) ≤ 0. Otherwise, σ∞ > 0. Then by Theorem 4.2 there exists Rc ∈ (R0 , R∞ ) such that ψ ′ < 0 in (R0 , Rc ) and ψ ′ > 0 in (Rc , R∞ ). We have ψ(Rc ) < min(ψ(R0 ), ψ(R∞ )) ≤ 0. Choose τ ∈ R such that ψ(Rc ) < τ < ˆ < I[ψ], where min(ψ(R0 ), ψ(R∞ )). Define ψˆ by ψˆ = ψ if ψ ≥ τ and ψˆ = τ if ψ < τ. Then I[ψ] I[·] is defined in (4.3). This is impossible as ψ minimizes I[·] over H = {φ ∈ H 1 (R0 , R∞ ) : ˆ Now, since σ > 0 and σ∞ ≤ 0, we have ψ ′ < 0 by Theorem 4.2. φ(R∞ ) = 0} that includes ψ. Since ψ(R∞ ) = 0, ψ > 0 in (R0 , R∞ ). Hence V ′ (ψ) > 0 in (R0 , R∞ ). Consequently, 2 1 ψ ′′ (r) = − ψ ′ (r) + V ′ (ψ(r)) > 0 in (R0 , R∞ ). r ε Since V ′ (ψ) > 0 in (R0 , R∞ ), we also have T ′ (r) < 0 for all r ∈ (R0 , R∞ ). (3) This can be proved similarly. Q.E.D. We now consider the case R∞ = ∞. Theorem 4.3. Let ε > 0, σ ∈ R, and R∞ = ∞. 20
(1) Suppose V : R → R is smooth and V ′ (0) 6= 0. Then there is no solution to the boundary-value problem (4.1). (2) Suppose V : R → R is smooth, V ′′ > 0 on R, V (−∞) = ∞ and V (∞) = ∞, and V (0) = inf s∈R V (s) = 0. Then there exists a unique solution to the boundary-value problem (4.1). Moreover, this solution ψ = ψ(|x|) is radially symmetric and smooth on Ω, ψ ′ < 0 and ψ ′′ > 0 in r = |x| if σ > 0, ψ ′ > 0 and ψ ′′ < 0 in r = |x| if σ > 0, and ψ = 0 for all r if σ = 0. We remark that in general all charges in an ionic system are in a finite region. Therefore it is natural to impose the far-field condition ψ(R∞ ) = 0. The properties on V assumed in part (2) of the theorem are exactly those covered in the third case in Theorem 3.3 but with the minimum of V attained exactly at 0. The theorem states that the solution to (4.1) exists only when V ′ (0) = 0. This is precisely the condition of charge neutrality in the bulk. Therefore, for the three cases of V stated in Theorem 3.3, only the third case with V ′ (0) = 0 permits a solution to (4.1). Proof of Theorem 4.3. (1) Assume V ′ (0) > 0. (The case that V ′ (0) < 0 is similar.) Assume ψ is a solution to (4.1), which is not assumed to be radially symmetric. Clearly, ψ is smooth on Ω. Since ψ(∞) = 0 and V ′ (0) > 0, there exist A ≥ R0 and α > 0 such that ε∆ψ = V ′ (ψ) ≥ α for all x ∈ Ω with |x| ≥ A. In the spherical coordinates (r, θ, ϕ) ˆ θ, ϕ), and (r ≥ 0, 0 ≤ θ ≤ π, 0 ≤ ϕ < 2π), ψ = ψ(r, ! ! ˆ ˆ 1 ∂ ∂ 2 ψˆ ∂ ψ 1 ∂ ψ 1 ∂ ∆ψ = 2 ≥α ∀r ≥ A. (4.9) r2 + 2 sin θ + 2 2 r ∂r ∂r r sin θ ∂θ ∂θ r sin θ ∂ϕ2 Define
1 u(r) = 4π
Z
2π 0
Z
π
ˆ θ, ϕ) dθdϕ sin θψ(r,
∀r > R0 .
0
Clearly, u(∞) = 0, since ψˆ → 0 as r → ∞. Multiplying both sides of the inequality in (4.9) by sin θ/(4π), and then integrating the resulting sides over θ ∈ [0, π] and ϕ ∈ [0, 2π), we ′ obtain (1/r2 ) (r2 u′ (r)) ≥ α for all r ≥ A. This leads to 1 1 3 1 ′ 2 ′ u (r) ≥ αr + A u (A) − αA ∀r ≥ A. 3 3 r2 With one more time of integration, we get 1 1 3 1 1 2 ′ 2 2 u(r) ≥ u(A) + α r − A − A u (A) − αA − 6 3 r A
∀r ≥ A.
This contradicts u(∞) = 0. Part (1) is proved. (2) We first show the uniqueness. Suppose ψ1 and ψ2 are both solutions to (4.1). Clearly they are smooth on Ω = {x ∈ R3 : |x| ≥ R0 }. Let η = ψ1 − ψ2 . This function attains its maximum on Ω. Suppose P := maxx∈Ω η > 0. If there exists x0 ∈ Ω = {x ∈ R3 : |x| > R0 } such that η(x0 ) = P > 0, then there exists an open ball B(x0 ) ⊂ Ω centered at x0 such that 21
η > 0 on B(x0 ). Since ψ1 and ψ2 are bounded on Ω, ψ1 (x), ψ2 (x) ∈ K for some compact subset of R for all x ∈ Ω. Setting γ0 = mins∈K V ′′ (s) > 0, we obtain ε∆η = V ′ (ψ1 ) − V ′ (ψ2 ) ≥ γ0 η > 0
in B(x0 ).
(4.10)
Since η(x0 ) = maxx∈B(x0 ) η > 0, the Strong Maximum Principle (cf. Theorem 3.5 of [13]) then implies that η is a constant in B(x0 ). Consequently ∆η = 0 in B(x0 ), and (4.10) implies η = 0 in B(x0 ), contradicting to η(x0 ) > 0. Thus there exists y0 ∈ ∂Ω with |y0 | = R0 such that η(y0 ) = P > η(x) for all x ∈ Ω. But this and (4.10) imply ∂n η(y0 ) = 0 by Lemma 3.4 in [13], contradicting ∂n η = 0 on ΓD = {x ∈ R3 : |x| = R0 }. Hence η ≤ 0 and ψ1 ≤ ψ2 in Ω. Similarly ψ2 ≤ ψ1 in Ω. Thus ψ1 = ψ2 in Ω. We now prove the existence of solution to (4.1) that satisfies all the properties. If σ = 0 then ψ = 0 in Ω is the desired solution. Assume σ > 0. (The case σ < 0 is similar.) Choose a sequence Rk ∈ R (k = 1, 2, . . . ) such that R0 < R1 < R2 < · · · < Rk < · · · and Rk → ∞ as k → ∞. For each integer k ≥ 1, let ψk = ψk (|x)|) be the unique solution to the boundary-value problem ε∆ψk = V ′ (ψk ) in Ωk = {x ∈ R3 : R0 < |x| < Rk }, ∂ψk (4.11) ε =σ on ΓkN = {x ∈ R3 : |x| = R0 }, ∂n ψk = 0 on ΓkD = {x ∈ R: |x| = Rk }. By Corollary 4.1, ψk > 0, ψk′ < 0, and ψk′′ > 0 in (R0 , Rk ). Using the spherical coordinates, we have by the first equation in (4.11) that 2 ψk′′ (r) + ψk′ (r) = V ′ (ψk (r)) in (R0 , Rk ). r
(4.12)
Multiply both sides of this equation by ψk′ (r) and integrate the resulting sides over [R0 , Rk ]. Noting that ψk′ ψk′′ = (1/2)(ψk′2 )′ and ψk′ V ′ (ψk ) = (d/dr)V (ψk ), we then obtain by ψk′ (R0 ) = −σ/ε, ψk (Rk ) = 0 and V (0) = 0 that σ 2 Z Rk 2 1 2 ′ 2 (ψk (Rk )) − + (ψk′ (r)) dr = −V (ψk (R0 )). 2 ε R0 r Thus, since V (s) > 0 for any s > 0, we have 0 ≤ V (ψk (r)) ≤ σ 2 /(2ε2 ) for all r ∈ [R0 , Rk ] and k ≥ 1. Since V is positive and strictly increasing in (0, ∞) and V (∞) = ∞, and since ψk (R0 ) ≥ ψk (r) for all r ∈ [R0 , Rk ], we have by setting φc > 0 with V (φc ) = σ 2 /(2ε2 ) that 0 ≤ ψk (r) ≤ φc for all r ∈ [R0 , Rk ] and all k ≥ 1. Moreover, since ψk′′ > 0 in (R0 , Rk ), −σ/ε = ψk′ (R0 ) ≤ ψk′ (r) ≤ ψk′ (R∞ ) ≤ 0 for all r ∈ [R0 , Rk ] and all k ≥ 1. Define ψk = 0 in 1,∞ (Rk , ∞) (k = 1, 2, . . . ). Then the sequence {ψk }∞ (R0 , ∞). k=1 is bounded in W ∞ ∞ There now exists a subsequence {ψk,1 }k=1 of {ψk }k=1 such that ψk,1 → ψ in H 1 (R0 , R1 ) ∞ as k → ∞ for some ψ ∈ H 1 (R0 , R1 ). There exists a subsequence {ψk,2 }∞ k=1 of {ψk,1 }k=1 such 1 1 that ψk,2 → φ in H (R0 , R2 ) as k → ∞ for some φ ∈ H (R0 , R2 ) with φ = ψ on (R0 , R1 ). We define ψ = φ in (R1 , R2 ) so that ψ ∈ H 1 (R0 , R2 ). Repeating, we obtain ψ ∈ C([R0 , ∞)) with 22
∞ ψ ∈ H 1 (R0 , R) for any R > R0 and subsequences {ψk,j }∞ k=1 of {ψk }k=1 (j = 1, 2, . . . ), where ∞ for each j ≥ 1, the sequence {ψk,j+1 }∞ k=1 is a subsequence of {ψk,j }k=1 . Moreover, ψk,j → ψ in 1 ∞ ∞ H (R0 , Rj ) as k → ∞. Now {ψk,k }k=1 is a subsequence of {ψk }k=1 , and {ψk,k }∞ k=1 converges to ψ in H 1 (R0 , R), and hence also in C([R0 , R]), for any R > R0 . The limit ψ = ψ(|x|) is clearly radially symmetric. It is a (weak and hence strong) solution to the first two equations in (4.1). We show now that limr→∞ ψ(r) = 0. For each k ≥ 1, define −1 1 1 ψ (R ) 1 − 1 if R0 ≤ r ≤ Rk , − k 0 R0 Rk r Rk ηk (r) = 0 if r > Rk .
One easily verifies that ηk (|x|) satisfies ε∆ηk = 0 in Ωk (cf. (4.11) for the definition of Ωk ) and ηk = ψk on the boundary of Ωk . Since ε∆(ψk − ηk ) ≥ 0 in Ωk , we thus have by the Maximum Principle that ψk ≤ ηk in Ωk and hence in Ω = {x ∈ R3 : |x| > R0 }. Denote ∞ ∞ {ηk,k }∞ k=1 the subsequence of {ηk }k=1 that corresponds to {ψk,k }k=1 . Then we have 0 ≤ lim sup ψ(r) = lim sup lim ψk,k (r) ≤ lim sup lim ηk,k (r) = lim sup r→∞
r→∞
k→∞
r→∞
k→∞
r→∞
R0 ψ(R0 ) = 0, r
implying limr→∞ ψ(r) = 0. Since ψk′ ≤ 0 on (R0 , Rk ) for all k ≥ 1 ψ ′ ≤ 0 on (R0 , ∞). Let u(r) = ψ ′ (r). Then u(|x|) satisfies ∆u − c(r)u = 0 on {x ∈ R3 : |x| > R0 }, where c(r) = − [(1/ε)V ′′ (ψ(r)) + 2/r2 ] < 0, cf. (4.8). Since u ≤ 0 in (R0 , ∞), the Strong Maximum Principle implies that u 6= 0 in (R0 , ∞). Hence u = ψ ′ < 0 in (R0 , ∞). Since ψ ≥ 0 and V (ψ) ≥ 0, we also have ψ ′′ (r) = −(2/r)ψ ′ + V ′ (ψ) > 0 for all r > R0 . Q.E.D.
5
Numerical Results
We now present numerical results to illustrate our analysis. We consider minimizing numerically the free-energy functional (2.6) with (2.7) and (2.8). We choose Ω = {x ∈ R3 : R0 < |x| < R∞ },
ΓN = {x ∈ R3 : |x| = R0 },
ΓD = {x ∈ R3 : |x| = R∞ },
with R0 , R∞ ∈ R such that R0 < R∞ . We assume the dielectric coefficient ε and surface charge density σ are both constants. We also set ψD = 0 on ΓD . Note that our system is radially symmetric. We use the constrained optimization method developed in our previous work [38] to solve numerically our optimization problem. We set R0 = 10 ˚ A and R∞ = 80 ˚ A. We distribute uniformly some negative charges on the 2 sphere r = R0 with the surface charge density σ = −150/(4π 2 R02 ) e/˚ A . We consider three species of counterions K+ , Ca2+ , Al3+ , and one species of coions Cl− in the solution that occupies Ω. So M = 4. We label these ionic species by 1, 2, 3, and 4, respectively. We label the solvent (water) by 0. The corresponding valences and ionic volumes are (z1 , z2 , z3 , z4 ) = 3 A , respectively. The volume (+1, +2, +3, −1) and (v1 , v2 , v3 , v4 ) = (5.513 , 4.753 , 4.13 , 6.373 ) ˚ 23
˚3 . For these real ionic systems, it happens of a solvent (water) molecule is v0 = 2.753 A that α4 < α0 = 0 < α1 < α2 < α3 , where αk = zk /vk (0 ≤ k ≤ 4) with z0 = 0. We use a3 to approximate the volume of an ion of diameter a. The temperature is T = 300K. The dielectric coefficient is absorbed into the Bjerrum length lB = e2 /(4πεkB T ). We set the Bjerrum length to be lB = 7˚ A. We consider two cases of the numbers of ions in different species. Case I: (N1 , N2 , N3 , N4 ) = (60, 30, 20, 50). By simple calculations, we have that in this case the effective surface charge density σ∞ > 0. Case II: (N1 , N2 , N3 , N4 ) = (50, 50, 50, 50). In this case σ∞ < 0. In Figure 4, we plot our computational results of the electrostatic potential ψ and the equilibrium ionic concentrations. For Case I, the potential ψ is monotonically increasing. This is exactly predicted by Theorem 4.2, since in this case, σ < 0 and σ∞ > 0. Notice that the concentration of +3, which has the largest ratio of valence-to-volume among all the species, monotonically deceases, while that of −1 monotonically increases. These are predicted by Theorem 3.2. For Case II, σ < 0 and σ∞ > 0, and the behavior of the potential ψ is again as predicted by Theorem 4.2. We observe the stratification of the counterion concentrations in both cases.
6
Conclusions
We have analyzed a PB-like mean-filed model for electrostatic interactions of multivalent ions with finite ionic sizes in an ionic solution near a charged surface. Our work continues that presented in [22] but gives more details on how equilibrium ionic concentrations depend on the electrostatic potential and how the important parameters, the valence-to-volume ratios of ions, determine the properties of ionic concentrations near the charged surface. Here we discuss some of our main results. First, our analysis shows that for a strong potential, the layering structure of counterion concentrations does depend on the valenceto-volume ratios. In general, the comparison of concentrations for different ionic species can be subtle. Instead of comparing ci and cj , one may need to compare (vi ci )1/vi and (vj cj )1/vj . The difference of these quantities may only be observed for certain range of the potential. See part (2) of Theorem 3.2. Second, we characterize the ionic charge density through the function V defined in (3.2); cf. also (3.15). We recall for the classical PB theory that the corresponding function V is defined by X −βqi φ V ′ (φ) = − qi c ∞ , i e i∈J,i6=0
c∞ i
where (i ∈ J, i 6= 0) are positive constants. An important difference here is that, when the ionic size effects are included, the function V grows linearly at the −∞ and ∞, rather than exponentially as in the classical PB theory. Third, when the ionic size effects are included, the Debye length is then modified. In ˆ D defined in (3.23) with Table 1 and Table 2, we compare the size-modified Debye length λ the classical Debye length λD defined in (3.24) for a few ionic systems. The difference between these two tables is that the diameters of ions listed in Table 1 do not, while those 24
0.08 25
0
Concentrations (M)
Electrostatic p otential (Volt)
0.1
−0.1 −0.2 −0.3 −0.4 −0.5
0.06
0.04
20
0.02 15 0 65
70
75
80
+1 +2 +3 −1
10
5
−0.6
10
20
30
40
50
60
70
0 10
80
11
12
13
14
15
16
17
18
19
20
Radial distance ( ˚ A)
Radial distance ( ˚ A)
0.08 25
0
Concentrations (M)
Electrostatic p otential (Volt)
0.1
−0.1
−0.2
−0.3
0.04
20
0.02 15 0 65
70
75
80
+1 +2 +3 −1
10
5
−0.4
−0.5 10
0.06
20
30
40
50
60
70
0 10
80
11
12
13
14
15
16
17
18
19
20
Radial distance ( ˚ A)
Radial distance ( ˚ A)
Figure 4: The electrostatic potential ψ (left) and the equilibrium ionic concentrations (right). Top. Case I: (N1 , N2 , N3 , N4 ) = (60, 30, 20, 50), σ < 0, and σ∞ > 0. Bottom. Case II: (N1 , N2 , N3 , N4 ) = (50, 50, 50, 50), σ < 0, and σ∞ < 0. The concentrations profiles are marked by the corresponding ionic valences.
25
in Table 2, include the size of a hydration shell. It is clear that the modification of the Debye length due to the inclusion of the ionic size effect is not significant. This may suggest that a linearized system may not describe well the size effect. Clearly, further tests are needed on real ionic solutions and biomolecular systems. Table 1: A comparison of the Debye length λD and the ionic size-modified Debye length ˆ D for a few systems of ionic solutions. The diameter of an ion does not include that of a λ water molecule (hydration shell). Concentrations mean bulk concentrations. Ionic systems Cl− , Na+ Cl− , K+ Cl− , Ca2+ Cl− , Na+ , K+ Cl− , Na+ , Ca2+ Cl− , Na+ , K+ , Ca2+ Cl− , K+ , Ca2+ , Al3+
Diameters 3.62, 2.04 3.62, 2.76 3.62, 2.00 3.62, 2.04, 3.62, 2.04, 3.62, 2.04, 3.62, 2.76,
of ions (˚ A)
2.76 2.00 2.76, 2.00 2.00, 1.35
Concentrations (M) 0.1, 0.1 0.1, 0.1 0.2, 0.1 0.2, 0.1, 0.1 0.3, 0.1, 0.1 0.4, 0.1, 0.1, 0.1 0.6, 0.1, 0.1, 0.1
λD ( ˚ A) 9.715 9.715 5.609 6.870 4.858 4.345 3.072
ˆ D (˚ λ A) 9.726 9.720 5.618 6.880 4.870 4.358 3.085
Table 2: A comparison of the Debye length λD and the ionic size-modified Debye length ˆ D for a few systems of ionic solutions. The diameter of an ion includes that of a water λ molecule (hydration shell). Concentrations mean bulk concentrations. Ionic systems Cl− , Na+ Cl− , K+ Cl− , Ca2+ Cl− , Na+ , K+ Cl− , Na+ , Ca2+ Cl− , Na+ , K+ , Ca2+ Cl− , K+ , Ca2+ , Al3+
Diameters 6.37, 4.79 6.37, 5.51 6.37, 4.75 6.37, 4.79, 6.37, 4.79, 6.37, 4.79, 6.37, 5.51,
of ions (˚ A)
5.51 4.75 5.51, 4.75 4.75, 4.10
Concentrations (M) 0.1, 0.1 0.1, 0.1 0.2, 0.1 0.2, 0.1, 0.1 0.3, 0.1, 0.1 0.4, 0.1, 0.1, 0.1 0.6, 0.1, 0.1, 0.1
λD ( ˚ A) 9.715 9.715 5.609 6.870 4.858 4.345 3.072
ˆ D (˚ λ A) 9.847 9.763 5.701 6.970 4.974 4.449 3.172
Finally, for a spherical charged molecule immersed in an ionic solution occupying a finite region, the qualitative behavior of the electrostatic potential—either it is monotonic or it changes the monotonicity only once in the entire range of the radial variable—is completely determined by the sign of two surface charge densities. One is the prescribed surface charge density and the other is the effective surface charge density. If the region of ionic solution is the entire complement of the sphere in R3 , then the ionic charge neutrality is necessary for the existence and uniqueness of the electrostatic potential that vanishes at infinity, governed by the PBE (classical or size-modified). Such a potential is always monotonic. 26
We conclude with two remarks on the mathematical model we have studied here. First, for the case of a uniform size for all ions and solvent molecules, the free-energy functional (2.1) can be derived from a lattice-gas model; cf. [20]. For a general and more interesting case where the sizes are different, such a derivation seems not available. It is therefore interesting to give a rigorous derivation of such a free-energy functional using the notion of equilibrium statistical mechanics. Second, we have assumed mainly that the region Ω of ionic solution is bounded. The case that Ω has an infinite volume is a tricky situation to define the functional F [c]; cf. (2.1) and (2.6). Usually, the concentration ci is taken to be the bulk value c∞ i (> 0) at infinity (or far away from the charged surface). Such a value can be experimentally determined. If Ω is unbounded with an infinite volume, then the integral of ci over Ω will be just infinite if ci = c∞ i at infinity. The PBE, however, can be defined on an unbounded domain as the electrostatic potential ψ can exist in the entire space.
Acknowledgment. This work was supported by the National Science Foundation (NSF) through grant DMS0811259 (B.L.), the NSF Center for Theoretical Biological Physics (CTBP) through grant PHY-0822283 (B.L.), the National Institutes of Health through the grant R01GM096188 (B.L.), the Natural Science Foundation of China through grant NSFC-11101276 and NSFC91130012 (P.L. and Z.X.), and the Alexander von Humboldt Foundation (Z.X.). The authors thank the referee for helpful comments.
References [1] R. Adams. Sobolev Spaces. Academic Press, New York, 1975. [2] D. Andelman. Electrostatic properties of membranes: The Poisson–Boltzmann theory. In R. Lipowsky and E. Sackmann, editors, Handbook of Biological Physics, volume 1, pages 603–642. Elsevier, 1995. [3] Y. Bai, K. Travers, V. B. Chu, J. Lipfert, S. Doniach, and D. Herschlag. Quantitative and comprehensive decomposition of the ion atmosphere around nucleic acids. J. Amer. Chem. Soc., 129:14981–14988, 2007. [4] M. L. Bleam, C. F. Anderson, and M. T. J. Record. Relative binding affinities of monovalent cations for double-stranded DNA. Proc. Natl. Acad. Sci. USA., 77:3085– 3089, 1980. [5] I. Borukhov, D. Andelman, and H. Orland. Steric effects in electrolytes: A modified Poisson–Boltzmann equation. Phys. Rev. Lett., 79:435–438, 1997. [6] I. Borukhov, D. Andelman, and H. Orland. Adsorption of large ions from an electrolyte solution: A modified Poisson–Boltzmann equation. Electrochimica Acta, 46:221–229, 2000. 27
[7] A. H. Boschitsch and P. V. Danilov. Formulation of a new and simple nonuniform size-modified Poisson–Boltzmann description. J. Comput. Chem., 33:1152–1164, 2012. [8] 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:3058–3069, 2008. [9] V. B. Chu, Y. Bai, J. Lipfert, D. Herschlag, and S. Doniach. Evaluation of ion binding to DNA duplexes using a size-modified Poisson–Boltzmann theory. Biophys. J, 93:3202– 3209, 2007. [10] M. E. Davis and J. A. McCammon. Electrostatics in biomolecular structure and dynamics. Chem. Rev., 90:509–521, 1990. [11] B. Eisenberg, Y. K. Hyon, and C. Liu. Energy variational analysis of ions in water and channels: Field theory for primitive models of complex ionic fluids. J. Chem. Phys., 133:104104, 2011. [12] F. Fogolari, A. Brigo, and H. Molinari. The Poisson–Boltzmann equation for biomolecular electrostatics: a tool for structural biology. J. Mol. Recognit., 15:377–392, 2002. [13] D. Gilbarg and N. S. Trudinger. Elliptic Partial Differential Equations of Second Order. Springer-Verlag, 2nd edition, 1998. [14] P. Grochowski and J. Trylska. Continuum molecular electrostatics, salt effects and counterion binding—A review of the Poisson–Boltzmann model and its modifications. Biopolymers, 89:93–113, 2008. [15] A. Y. Grosberg, T. T. Nguyen, and B. I. Shklovskii. Colloquium: The physics of charge inversion in chemical and biological systems. Rev. Mod. Phys., 74:329–345, 2002. [16] B. Hille. Ion Channels of Excitable Membranes. Sinauer Associates, 3rd edition, 2001. [17] J. J. Howard, J. S. Perkyns, and B. M. Pettitt. The behavior of ions near a charged wall–dependence on ion size, concentration, and surface charge. J. Phys. Chem. B, 114:6074–6083, 2010. [18] M. S. Kilic, M. Z. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes at large applied voltages. I. Double-layer charging. Phys. Rev. E, 75:021502, 2007. [19] M. S. Kilic, M. Z. Bazant, and A. Ajdari. Steric effects in the dynamics of electrolytes at large applied voltages. II. Modified Poisson–Nernst–Planck equations. Phys. Rev. E, 75:021503, 2007. [20] V. Kralj-Igliˇc and A. Igliˇc. A simple statistical mechanical approach to the free energy of the electric double layer including the excluded volume effect. J. Phys. II (France), 6:477–491, 1996.
28
[21] D. Lambert, D. Leipply, R. Shiman, and D. E. Draper. The influence of monovalent cation size on the stability of RNA tertiary structures. J. Mol Biol., 390:791–804, 2009. [22] B. Li. Continuum electrostatics for ionic solutions with nonuniform ionic sizes. Nonlinearity, 22:811–833, 2009. [23] B. Li. Minimization of electrostatic free energy and the Poisson–Boltzmann equation for molecular solvation with implicit solvent. SIAM J. Math. Anal., 40:2536–2566, 2009. [24] S. W. Lockless, M. Zhou, and R. MacKinnon. Structural and thermodynamic properties of selective ion binding in a K + channel. PLoS Biology, 5:1079–1088, 2007. [25] J. J. L´opez-Garcia and Jos´e Horno. Poisson–Boltzmann description of the electrical double layer including ion size effects. Langmuir, 27:13970–13980, 2011. [26] B. Z. Lu and Y. C. Zhou. Poisson–Nernst–Planck equations for simulating biomolecular diffusion-reaction processes II: Size effects on ionic distributions and diffusion-reaction rates. Biophys J., 100:2475–2485, 2011. [27] A. C. Maggs. A minimizing principle for the Poisson–Boltzmann equation. EPL, 98:16012, 2012. ´ [28] M. Quesada-P´erez, A. Mart´ın-Molina, and R. Hidalgo-Alvarez. Simulation of electric double layers with multivalent counterions: Ion size effect. J. Chem. Phys., 121:8618– 8626, 2004. [29] R. B. Schoch, J. Han, and P. Renaud. Transport phenomena in nanofluidics. Rev. Mod. Phys., 80:839–883, 2008. [30] K. A. Sharp and B. Honig. Electrostatic interactions in macromolecules: Theory and applications. Annu. Rev. Biophys. Biophys. Chem., 19:301–332, 1990. [31] A. R. J. Silalahi, A. H. Boschitsch, R. C. Harris, and M. O. Fenley. Comparing the predictions of the nonlinear Poisson–Boltzmann equation and the ion size-modified Poisson–Boltzmann equation for a low-dielectric charged spherical cavity in an aqueous salt solution. J. Chem. Theory Comput., 6:3631–3639, 2010. [32] U. P. Strauss and Y. P. Leung. Volume changes as a criterion for site binding of counterions by polyelectrolytes. J. Amer. Chem. Soc., 87:1476–1480, 1965. [33] G. Tresset. Generalized Poisson–Fermi formalism for investigating size correlation effects with multiple ions. Phys. Rev. E, 78:061506, 2008. [34] M. Valisk´o, D. Boda, and D. Gillespie. Selective adsorption of ions with different diameter and valence and highly charged interfaces. J. Phys. Chem. C, 111:15575– 15585, 2007. 29
[35] V. Vlachy. Ionic effects beyond Poisson–Boltzmann theory. Annu. Rev. Phys. Chem., 50:145–165, 1999. [36] H. L. V¨ortler, K. Sch¨afer, and W. R. Smith. Simulation of chemical potentials and phase equilibria in two- and three-dimensional square-well fluids: finite size effects. J. Phys. Chem. B, 112:4656–4661, 2008. [37] J. Wen, S. Zhou, Z. Xu, and B. Li. Competitive adsorption and ordered packing of counterions near highly charged surfaces: From mean-field theory to Monte Carlo simulations. Phys. Rev. E, 85:041406, 2012. [38] S. Zhou, Z. Wang, and B. Li. Mean-field description of ionic size effects with nonuniform ionic sizes: A numerical approach. Phys. Rev. E, 84:021901, 2011.
30
