Home

Search

Collections

Journals

About

Contact us

My IOPscience

Ionic size effects: generalized Boltzmann distributions, counterion stratification and modified Debye length

This content has been downloaded from IOPscience. Please scroll down to see the full text. 2013 Nonlinearity 26 2899 (http://iopscience.iop.org/0951-7715/26/10/2899) View the table of contents for this issue, or go to the journal homepage for more Download details: IP Address: 141.161.91.14 This content was downloaded on 23/07/2017 at 13:18 Please note that terms and conditions apply.

You may also be interested in: Continuum electrostatics for ionic solutions with non-uniform ionic sizes Bo Li Asymptotic analysis of the Poisson–Boltzmann equation describing electrokinetics in porous media Grégoire Allaire, Jean-François Dufrêche, Andro Mikeli et al. Non-ideal electrostatic correlations in equilibrium electrolytes Alexandre Ern, Rémi Joubaud and Tony Lelièvre Eternal solutions to a singular diffusion equation with critical gradient absorption Razvan Gabriel Iagar and Philippe Laurençot Yukawa-field approximation of electrostatic free energy and dielectric boundary force Hsiao-Bing Cheng, Li-Tien Cheng and Bo Li High energy rotation type solutions of the forced pendulum equation Patricio Felmer, André de Laire, Salomé Martínez et al. Stationary patterns for an adsorbate-induced phase transition model: II. Shadow system Kousuke Kuto and Tohru Tsujikawa Coupled second order singular perturbations for phase transitions Margarida Baía, Ana Cristina Barroso, Milena Chermisi et al. Universality and critical behaviour in the chiral two-matrix model Steven Delvaux, Dries Geudens and Lun Zhang

IOP PUBLISHING

NONLINEARITY

Nonlinearity 26 (2013) 2899–2922

doi:10.1088/0951-7715/26/10/2899

Ionic size effects: generalized Boltzmann distributions, counterion stratification and modified Debye length Bo Li1 , Pei Liu2 , Zhenli Xu3 and Shenggao Zhou1 1 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 2 Department of Mathematics and Institute of Natural Sciences, Shanghai Jiao Tong University, 800 Dongchuan Rd., Shanghai, 200240, People’s Republic of China 3 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, People’s Republic of China

E-mail: [email protected], [email protected], [email protected] and [email protected]

Received 18 March 2013, in final form 19 July 2013 Published 23 September 2013 Online at stacks.iop.org/Non/26/2899 Recommended by C Le Bris 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. Both recent studies using a variational mean-field approach that includes ionic size effects and Monte Carlo simulations suggest that counterion stratification is determined by the ionic valence-to-volume 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 for 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 valence-to-volume ratios in the counterion stratification and the modification of Debye length due to the effect of ionic sizes. Mathematics Subject Classification: 35J20, 35J60, 35Q99, 82B21, 92C05 (Some figures may appear in colour only in the online journal)

0951-7715/13/102899+24$33.00

© 2013 IOP Publishing Ltd & London Mathematical Society

Printed in the UK & the USA

2899

2900

B Li et al

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 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 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 growing interest in studying ionic size effects within the PBlike 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  −1 vi ci (x) . (1.1) c0 (x) = v0 1 − i=1

The size-modified electrostatic free-energy functional is then given by [5, 20, 22, 26]  
M M   1 −1 F [c1 , . . . , cM ] = (1.2) ci [ln (vi ci ) − 1] − µi ci dV . ρψ + β 2 i=0 i=1  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 i = 1, . . . , M, (1.3) ln (v0 c0 ) − ln (vi ci ) = β (qi ψ − µi ) , v0 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∞ e−βqi ψ ci = i = 1, . . . , M,  , M 1 + j =1 vj cj∞ e−βqj ψ − 1 where vi−1 eβµi , i = 1, . . . , M. ci∞ =  βµj 1+ M j =1 e

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2901

Note that for each i, ci∞ 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 also has 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; see theorem 2.2. We also derive the conditions for the equilibrium concentrations; see (1.3). We then focus on the equilibrium conditions (1.3) to analyse the generalized Boltzmann distributions. In particular, we characterize the behaviour 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 smallest valence-to-volume ratios decreases and increases monotonically, respectively, with respect to the electrostatic potential ψ. As ψ → −∞ (or ∞), the concentration of ions with the largest (or smallest) valenceto-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 section 3.2. (3) A size-modified Debye length, λˆ D , is derived with a weak potential limit. A formula of of λˆ D is given in (3.23) in section 3.3. 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. Finally, we consider a charged spherical molecule of radius R0 in an ionic solution that occupies a region defined by {R0 < |x| < R∞ } for some R∞ > R0 , and study the PBE, both classical and size-modified. In the case where 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; see figure 4 in section 5. If R∞ = ∞, then only

2902

B Li et al

Ω

+ -

-

+

+

+

-

-

+

-

+

+ +

Γ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 +.

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 interest as it is known that the effect of a charged system of finite size can be significant, see 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 our 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 = ∅. 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; see 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 ionic species at position x ∈ . As before, ρ = M 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

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2903

electrostatic free-energy functional   M  M    1 1 −1 F [c] = ρψ dV + σ ψ dS + β ci [ln(vi ci ) − 1] dV − µi ci dV ,  2 N 2 i=0  i=1  (2.1)  ∇ · ε∇ψ = −ρ in ,  on D , (2.2) with ψ = ψ∞  ∂ψ on N . ε ∂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   in , in ,   ∇ · ε∇ψD = 0 ∇ · ε∇ψN = 0   ψD = ψ∞ ψ on D , = 0 on D , N and   ∂ψ ∂ψ D N ε ε   =0 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 problem of Poisson’s equation  ∇ · ε∇(Lf ) = −f in ,    Lf = 0 on D ,  ∂(Lf )  ε =0 on N . ∂n (The uniqueness is guaranteed if D has a nonzero surface measure.) Note 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) M where ψˆ = L( i=1 qi ci ). By integration by parts, we have    ∂ψN ∂ψN σ ψˆ dS = ε ψˆ dS = ε ψˆ dS ∂n ∂n N N    ˆ · ε∇ψN dV + = ε∇ ψˆ · ∇ψN dV = ε∇ ψˆ · ∇ψN dV ψ∇        M

ˆ = qi ci ψN dV . −∇ · ε∇ ψ ψN dV = 



i=1

Therefore the free-energy functional F [c] can now be rewritten as   M M    M      1 1  qi F [c] = ψD + ψN − µi ci dV qi ci L  qj cj  dV + 2  2 i=1 j =1 i=1  +β

−1

M   i=0

 ci [ln(vi ci ) − 1] dV + 

N

1 σ (ψN + ψD ) dS. 2

2904

B Li et al

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  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 i = 1, . . . , M. (2.5) ln(v0 c0 ) − ln(vi ci ) = β qi ψ − ψD − µi , v0 2 The following result can be proved by the argument used in [22]. Theorem 2.1. (1) The functional F [c] is convex. Moreover, it has a unique minimizer in the set of all  concentrations c1 , . . . , cM such that ci  0 (i = 1, . . . , M) and M 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   M   1 1 −1 ˆ F [c] = ci [ln(vi ci ) − 1] dV , (2.6) ρψ dV + σ ψ dS + β  2 N 2 i=0   ∇ · ε∇ψ = −ρ in ,    on D , ψ = ψ ∞ with (2.7)  ∂ψ  ε =σ on N ,  ∂n subject to ci dV = Ni , i = 1, . . . , M, (2.8) M



where again ρ = 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ˆ = (cˆ1 , . . . , 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     1 vi δci Fˆ [c] = qi ψ − ψD + β −1 ln(vi ci ) − i = 1, . . . , M. (2.9) ln (v0 c0 ) , 2 v0 For the minimizer c, ˆ we have 

δci Fˆ [c] ˆ di dV = 0 

 ∀di :

di dV = 0,

i = 1, . . . , M.



By lemma 2.1, 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).

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2905

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). 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   u v dV = 0 ∀v ∈ L2 () : v dV = 0. 



Then there exists a constant µ ∈ R such that u = µ almost everywhere in . Proof. Denote for any w ∈ L2 ()   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    u (v − v) dV = 0. (u − u) v dV = (u − u) (v − v) dV = 





Therefore u = u almost everywhere in . Set µ = u. The lemma is proved.

QED

3. 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 ) , i = 1, . . . , M. (3.1) v0 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 V (φ) = −

M   i=1

φ

qi Bi (ξ ) dξ + V0 ,

(3.2)

0

where V0 is a constant. Then the function V : R → R satisfies V  > 0 on R and is thus strictly convex.

2906

B Li et al

The relations ci = Bi (φ) (i = 1, . . . , M) define the generalized Boltzmann distributions for the equilibrium concentrations c1 , . . . , cM . In general, explicit formulas of such  distributions seem to be unavailable. Note that −V  (φ) = M i=1 qi Bi (φ) with φ = ψ − ψD /2 the ionic charge density. The proof of V  > 0 in R given in [22] (see the proof of lemma 5.2 in [22]) applies to our present case; see the remark in section 3.1. The electrostatic potential ψ is now the unique solution to the boundary-value problem of the generalized and implicit PBE   ψD  = 0. (3.3) ∇ · ε∇ψ − V ψ − 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 boundary-value 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−l z−1 z0 z1 zm < ··· < < =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 = 0. We thus have by (1.1) that    c0 = B0 (φ) = v0−1 1 − vi Bi (φ) . (3.6) i∈J,i=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 cm (φ) < 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

φ
(vi−1 ci−1 )1/vi−1

if and only if

φ>

µi−1 /vi−1 − µi /vi . qi−1 /vi−1 − qi /vi

(3) Asymptotic behaviour 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)

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2907

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)

∀i ∈ J.

(3.11)

This also implies (vi ci )1/vi = eβµi /vi e−β(qi /vi )φ (v0 c0 )1/v0

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 = j , i = 0, and j = 0. We have by (3.5) that  βqj βµj βµi 1  βqi 1 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 ci =

ηi vi /v0 −βqi φ γ e vi

∀i ∈ J.

(3.12)

This and the definition of c0 (see (1.1)) imply  ηi γ vi /v0 e−βqi φ = 1.

(3.13)

i∈J

Taking the derivative of both sides of (3.13) with respect to φ, and using our notation z0 = 0 and µ0 = 0, we obtain  βv0 γ i∈J ηi qi γ vi /v0 e−βqi φ  . γ  (φ) = 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  vi /v0 −βqi φ e i∈J ηi vi (qi /vi − qk /vk ) γ  vk /v0 −βqk φ  e . (3.14) ck (φ) = βηk γ v /v −βq φ i 0 i e i∈J ηi vi γ  (φ) < 0 and This, the assumption (3.4), and the definition qi = zi e (i ∈ J ) imply that cm  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 φ → ∞.

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. QED

2908

B Li et al

Remark. By (3.14) and the Cauchy–Schwarz inequality, we obtain through a series of calculations that  2     2 2 β q v c − q c v c  i i i i i i∈J i∈J i i∈J i qk ck (φ) =
0.  2 v c i∈J i i 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 φ; see 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  φ qi Bi (ξ ) dξ + V0 . V (φ) = − i∈J

 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 (∞) = ∞. In figure 2, we show the schematic behaviour of the convex function V for the three different cases. 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 ). Note 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 e−β(qi /vi )φ . = (vi ci )
max e (vi Bi (φ)) M

1
i
M

Therefore, the integral of i=1 qi Bi (φ) from 0 to φ > 0 is a bounded function of φ. Thus, by the definition of V (see (3.2)), V (φ) is bounded below. Since V  < 0 on R,

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2909

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 φ.

V decreases. Therefore, V (∞) exists. Choosing the constant  ∞ 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),   φ0 qi Bi (φ) dφ, V0 = i∈J

0

we have minφ∈R V (φ) = 0.

QED

2910

B Li et al

3.3. Modified Debye length We now linearize V  (φ) at φ = 0 to obtain a modified Debye length. We first define ci∞ = Bi (0) for any i ∈ J. Note that  v0 c0∞ = 1 − vi ci∞ . (3.16) i∈J,i=0

By (3.10) with φ = 0 and ηi = eβµi , we also have that  ∞ vi /v0 v0 c0 = ηi−1 vi ci∞ , i ∈ J.

(3.17)

We denote again γ = v0 c0 which depends on φ. We have for |φ| 1 that   γ = τ0 + τ1 φ + O φ 2 ,

(3.18)

where τ0 and τ1 are two constants to be determined. Using our notation, (1.1) becomes  v0 c0 = 1 − v i ci . (3.19) i∈J,i=0

This and (3.16) imply that τ0 = v0 c0 (0) = 1 −



vi ci∞ = v0 c0∞ .

(3.20)

i∈J,i=0

Substituting (3.18) into (3.12) for a fixed i ∈ J with i = 0, we obtain by (3.20) and (3.17) that  vi /v0    ηi  ci = τ0 + τ 1 φ + O φ 2 1 − βqi φ + O φ 2 vi     vi /v0    ηi vi /v0 τ1 = τ0 1 + φ + O φ2 1 − βqi φ + O φ 2 vi τ0        τ v i 1 = ci∞ 1 + φ + O φ2 1 − βqi φ + O φ 2 v0 τ0     v i τ1 = ci∞ + ci∞ (3.21) − βqi φ + O φ 2 . v0 τ0 It now follows from (3.18), (3.19), (3.21) and (3.20) that        2 vi τ1 ∞ ∞ τ0 + τ1 φ + O φ = 1 − vi ci − vi ci − βqi φ + O φ 2 . 2 ∞ v c 0 0 i∈J,i=0 i∈J,i=0 Comparing the leading-order terms, we obtain exactly (3.20). Comparing the O(φ)-terms and noting that z0 = 0, we obtain  βv02 c0∞ i∈J qi vi ci∞ τ1 = .  2 ∞ i∈J vi ci This and (3.21) then lead to    vi j ∈J qj vj cj∞   ∞ ∞ ci = ci + βci − qi φ + O φ 2  2 ∞ j ∈J vj cj Consequently, V  (φ) = −

 i∈J

qi ci = −

 i∈J


 qi ci∞ − β

i∈J



qi vi ci∞

2 ∞ i∈J vi ci

∀i ∈ J 2 −

 i∈J

with i = 0. 

  qi2 ci∞ φ + O φ 2 .

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2911

In the small |φ| approximation, we then obtain by (3.3) and the shift φ = ψ − ψD /2 the linearized PBE with ionic size effects  qi ci∞ − 1 ε λˆ −2 (3.22) ∇ · ε∇ψ − ε λˆ −2 D ψ =− D ψD , 2

i∈J

where λˆ D > 0 is defined by
2   qi vi ci∞ β  2 ∞ i∈J −2 . q c −  λˆ D = 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       qi vi ci∞
qi2 ci∞ vi2 ci∞ . i∈J

i∈J

i∈J

 In the case that ε is a constant, ψD = 0 on D , and i∈J qi ci∞ = 0 which is the charge neutrality, the linearized equation (3.22) simplifies to the familiar equation ψ − λˆ −2 D ψ = 0. We call λˆ D the size-modified Debye length. Recall that the Debye length λD > 0 in the classical Debye–H¨uckel theory is defined by β 2 ∞ q c . (3.24) λ−2 D = ε i∈J i i It is clear that λˆ D > λD . This indicates that finite ionic sizes lead to a longer Debye screening length. 4. A charged spherical molecule: basic properties In this section, we consider a charged spherical molecule, centred 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 },  if R∞ < ∞, {x ∈ R3 : |x| = R∞ } D = {∞} if R∞ = ∞. We consider the boundary-value problem of the PBE, classical and size-modified, for the electrostatic potential ψ (see (2.2) and (3.3)):  in , εψ = V  (ψ)    ∂ψ (4.1) ε =σ on N ,   ψ∂n =0 on . D

If ψ = ψ(r) (r = |x|) is radially symmetric, then this system is the same as    2     if R0 < r < R∞ , ε ψ + ψ = V  (ψ)   r  σ ψ  (R0 ) = − ,    ε   ψ(R∞ ) = 0.

(4.2)

2912

B Li et al

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 = 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 (3.4).  For the classical PBE, we have V (φ) = i∈J ci∞ (e−βqi φ − 1), where ci∞ (i ∈ J, i = 0) are bulk ionic concentrations. Under the assumption on the charge neutrality i∈J qi ci∞ = 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 σ = 0. The case that σ = 0 can be proved similarly. We denote H = {φ ∈ H 1 (R0 , R∞ ) : φ(R∞ ) = 0} and define I : H → (−∞, ∞] by  R∞   ε  2 I [φ] = 4π ∀φ ∈ H. (4.3) φ (r) + V (φ) r 2 dr − 4π R02 σ φ(R0 ) 2 R0 If φ ∈ H , then we have with δ = ε/(4σ ) > 0 that  R∞   R∞ 1/2   |φ(R0 )| =  φ  (r) dr 
R∞ − R0 |φ  (r) |2 dr R0

R∞ − R0 +δ
4δ



R0

R∞

|φ  (r) |2 dr,

R0

leading to  I [φ]  πεR02

R∞

R0

|φ  (r) |2 dr −

4π σ 2 R02 (R∞ − R0 ) ε

∀φ ∈ H.

Note that by Poincar´e’s inequality, I [φ] controls φ 2H 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 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     ε Q[φ] = |∇φ|2 + V (φ) dV − σ φ dS  2 N

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2913

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     ε 2 |∇η| + ε∇ψ · ∇η + V (ψ + η) − V (ψ) dV − σ η dS Q[ψ + η] − Q[ψ] =  2 N   " !  σ η dS ε∇ψ · ∇η + V  (ψ)η dV −  N    " ! ∂ψ  = ε σ η dS −εψ + V (ψ) η dV + η 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, ψˆ = ψ. QED We now study the behaviour 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}:  r V  (ψ(s))s 2 ds ∀r ∈ [R0 , R∞ ]. (4.5) T (r) = −4π 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). (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∞ ). 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∞ ).

2914

B Li et al

Figure 3. The graph of electrostatic potential ψ.

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 (see (4.4)) and the strict convexity of V , the restriction of ψ onto ω is the unique minimizer among all φ ∈ H 1 (ω) of the functional    ε (4.7) |∇φ|2 + V (φ) dV . G[φ] = ω 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, ψ  = 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 1  2 V (ψ(r)) + 2 ψ  (r) ∀r ∈ (R0 , t). (4.8) ψ  (r) + ψ  (r) = r ε r Hence, setting u = ψ  , we have Lu := u + c(r)u = 0 for R0 < r < t, where c(r) := −[(1/ε)V  (ψ(r)) + 2/r 2 ] < 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

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2915

value 0 at any interior point x with |x| = s in {x ∈ R3 : R0 < |x| < t}. By the strong maximum principle (see 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∞ ). QED 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 ) = 0, a contradiction. If V  < 0, we obtain the same contradiction ψ  (r0 ) = 0 by considering (R0 , r0 ). Finally, assume V  (φ0 ) = 0. If ψ(r0 ) = φ0 , then since ψ  (r0 ) = 0, ψ  (r0 ) = V  (φ(r0 ))/ε = 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. QED We recall that for any r ∈ (R0 , R∞ ), T (r) is the total amount of ionic charges in the region {x ∈  : R0 < |x| < r}; see (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 = 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∞ ).

2916

B Li et al

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 ) < τ < min(ψ(R0 ), ψ(R∞ )). Define ψˆ by ψˆ = ψ if ψ  τ and ψˆ = τ if ψ < τ. ˆ < I [ψ], where I [·] is defined in (4.3). This is impossible as ψ minimizes Then I [ψ] ˆ Now, since σ > 0 and I [·] over H = {φ ∈ H 1 (R0 , R∞ ) : φ(R∞ ) = 0} that includes ψ. σ∞
0, we have ψ  < 0 by theorem 4.2. Since ψ(R∞ ) = 0, ψ > 0 in (R0 , R∞ ). Hence V  (ψ) > 0 in (R0 , R∞ ). Consequently, 2 1 in (R0 , R∞ ). ψ  (r) = − ψ  (r) + V  (ψ(r)) > 0 r ε Since V  (ψ) > 0 in (R0 , R∞ ), we also have T  (r) < 0 for all r ∈ (R0 , R∞ ). (3) This can be proved similarly.

QED

We now consider the case R∞ = ∞. Theorem 4.3. Let ε > 0, σ ∈ R, and R∞ = ∞. (1) Suppose V : R → R is smooth and V  (0) = 0. Then there is no solution to the boundaryvalue 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 ∂ψ ψ = 2  α ∀r  A. (4.9) r + 2 sin θ + 2 r ∂r ∂r r sin θ ∂θ ∂θ r 2 sin θ ∂ϕ 2 Define 1 u(r) = 4π

 0

2π

 0

π

ˆ θ, ϕ) dθ dϕ sin θ ψ(r,

∀r > R0 .

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2917

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/r 2 )(r 2 u (r))  α for all r  A. This leads to   1 1 1 ∀r  A. u (r)  αr + A2 u (A) − αA3 2 3 3 r With one more time of integration, we get     1 1 1 1  − u(r)  u(A) + α r 2 − A2 − A2 u (A) − αA3 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 ) ⊂  centred at x0 such that η > 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 (see 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 a 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 Nk = {x ∈ R3 : |x| = R0 },  ∂n   ψ = 0 on Dk = {x ∈ R: |x| = Rk }. k 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 in (R0 , Rk ). (4.12) ψk (r) + ψk (r) = V  (ψk (r)) r Multiply both sides of this equation by ψk (r) and integrate the resulting sides over [R0 , Rk ]. Noting that ψk ψk = (1/2)(ψk2 ) and ψk V  (ψk ) = (d/dr)V (ψk ), we then obtain by ψk (R0 ) = −σ/ε, ψk (Rk ) = 0 and V (0) = 0 that  σ 2   R k 2  2 1  2 (ψk (Rk )) − ψk (r) dr = −V (ψk (R0 )). + 2 ε r R0 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,

2918

B Li et al

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 (Rk , ∞) (k = 1, 2, . . .). Then the sequence {ψk }∞ k=1 is bounded in W 1,∞ (R0 , ∞). ∞ There now exists a subsequence {ψk,1 }∞ k=1 of {ψk }k=1 such that ψk,1 → ψ in 1 1 H (R0 , R1 ) as k → ∞ for some ψ ∈ H (R0 , R1 ). There exists a subsequence {ψk,2 }∞ k=1 1 1 of {ψk,1 }∞ k=1 such 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 ψ ∈ 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 ∞ 1 is a subsequence of {ψk,j }k=1 . Moreover, ψk,j → ψ in H (R0 , Rj ) as k → ∞. Now ∞ ∞ 1 {ψk,k }∞ k=1 is a subsequence of {ψk }k=1 , and {ψk,k }k=1 converges to ψ in H (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 1 1  if R0
r
Rk , − − ψk (R0 ) ηk (r) = R0 Rk r Rk  0 if r > Rk . One easily verifies that ηk (|x|) satisfies εηk = 0 in k (see (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 R0 ψ(R0 ) = 0, 0
lim sup ψ(r) = lim sup lim ψk,k (r)
lim sup lim ηk,k (r) = lim sup r r→∞ r→∞ k→∞ r→∞ k→∞ 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/r 2 ] < 0, see (4.8). Since u
0 in (R0 , ∞), the strong maximum principle implies that u = 0 in (R0 , ∞). Hence u = ψ  < 0 in (R0 , ∞). Since ψ  0 QED and V (ψ)  0, we also have ψ  (r) = −(2/r)ψ  + V  (ψ) > 0 for all r > R0 .

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 Å and R∞ = 80 Å. We distribute uniformly some negative charges 2 on the sphere r = R0 with the surface charge density σ = −150/(4π 2 R02 ) e/Å . 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 ) = (+1, +2, +3, −1) and (v1 , v2 , v3 , v4 ) = (5.513 , 4.753 , 4.13 , 6.373 ) Å3 ,

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2919

0.08

25

0

Concentrations (M)

Electrostatic p otential (Volt)

0.1

−0.1 −0.2 −0.3 −0.4 −0.5

0.06

20

0.04

0.02

15 0 65

70

75

80

10

+1 +2 +3 −1

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 (Å)

Radial distance (Å)

0.08

25 0

Concentrations (M)

Electrostatic p otential (Volt)

0.1

−0.1

−0.2

−0.3

20

0.04

0.02

15 0 65

70

75

80

10

+1 +2 +3 −1

5

−0.4

−0.5 10

0.06

20

30

40

50

60

Radial distance (Å)

70

80

0 10

11

12

13

14

15

16

17

18

19

20

Radial distance (Å)

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.

respectively. The volume of a solvent (water) molecule is v0 = 2.753 Å3 . For these real ionic systems, it happens that α4 < α0 = 0 < α1 < α2 < α3 , where αk = zk /vk (0
k
4) with z0 = 0. We use a 3 to approximate the volume of an ion of diameter a. The temperature is T = 300 K. The dielectric coefficient is absorbed into the Bjerrum length lB = e2 /(4πεkB T ). We set the Bjerrum length to be lB = 7 Å. 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. Note 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 behaviour 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 analysed a PB-like mean-field model for electrostatic interactions of multivalent ions with finite ionic sizes in an ionic solution near a charged surface. Our work continues that

2920

B Li et al 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

Diameters of ions (Å)

Concentrations (M)

λD (Å)

λˆ D (Å)

Cl− ,

3.62, 2.04 3.62, 2.76 3.62, 2.00 3.62, 2.04, 2.76 3.62, 2.04, 2.00 3.62, 2.04, 2.76, 2.00 3.62, 2.76, 2.00, 1.35

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

9.715 9.715 5.609 6.870 4.858 4.345 3.072

9.726 9.720 5.618 6.880 4.870 4.358 3.085

Na+

Cl− , K+ Cl− , Ca2+ Cl− , Na+ , K+ Cl− , Na+ , Ca2+ Cl− , Na+ , K+ , Ca2+ Cl− , K+ , Ca2+ , Al3+

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 valence-tovolume 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); see also (3.15). We recall for the classical PB theory that the corresponding function V is defined by  V  (φ) = − qi ci∞ e−βqi φ , i∈J,i=0

ci∞

(i ∈ J, i = 0) are positive constants. An important difference here is that, when the where 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 tables 1 and 2, we compare the size-modified Debye length λˆ D defined in (3.23) with 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 include the size of a hydration shell, while those in table 2 do. 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 the size effect well. Clearly, further tests are needed on real ionic solutions and biomolecular systems. Finally, for a spherical charged molecule immersed in an ionic solution occupying a finite region, the qualitative behaviour 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. 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)

Generalized Boltzmann distributions, counterion stratification and modified Debye length

2921

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

Diameters of ions (Å)

Concentrations (M)

λD (Å)

λˆ D (Å)

Cl− ,

6.37, 4.79 6.37, 5.51 6.37, 4.75 6.37, 4.79, 5.51 6.37, 4.79, 4.75 6.37, 4.79, 5.51, 4.75 6.37, 5.51, 4.75, 4.10

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

9.715 9.715 5.609 6.870 4.858 4.345 3.072

9.847 9.763 5.701 6.970 4.974 4.449 3.172

Na+

Cl− , K+ Cl− , Ca2+ Cl− , Na+ , K+ Cl− , Na+ , Ca2+ Cl− , Na+ , K+ , Ca2+ Cl− , K+ , Ca2+ , Al3+

can be derived from a lattice-gas model; see [20]. For a general and more interesting case where the sizes are different, such a derivation seems unavailable. 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 mainly assumed 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]; see (2.1) and (2.6). Usually, the concentration ci is taken to be the bulk value ci∞ (> 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 just be infinite if ci = ci∞ at infinity. The PBE, however, can be defined on an unbounded domain as the electrostatic potential ψ can exist in the entire space. Acknowledgments This work was supported by the National Science Foundation (NSF) through grant DMS1319731 (BL), the NSF Center for Theoretical Biological Physics (CTBP) through grant PHY-0822283 (BL), the National Institutes of Health through the grant R01GM096188 (BL), the Natural Science Foundation of China through grant NSFC-11101276 and NSFC-91130012 (PL and ZX), and the Alexander von Humboldt Foundation (ZX). The authors thank the referee for helpful comments. References [1] Adams R 1975 Sobolev Spaces (New York: Academic) [2] Andelman D 1995 Electrostatic properties of membranes: the Poisson–Boltzmann theory Handbook of Biological Physics vol 1, ed R Lipowsky and E Sackmann (Amsterdam: Elsevier) pp 603–42 [3] Bai Y, Travers K, Chu V B, Lipfert J, Doniach S and Herschlag D 2007 Quantitative and comprehensive decomposition of the ion atmosphere around nucleic acids J. Am. Chem. Soc. 129 14981–8 [4] Bleam M L, Anderson C F and Record M T J 1980 Relative binding affinities of monovalent cations for doublestranded DNA Proc. Natl Acad. Sci. USA 77 3085–9 [5] Borukhov I, Andelman D and Orland H 1997 Steric effects in electrolytes: a modified Poisson–Boltzmann equation Phys. Rev. Lett. 79 435–8 [6] Borukhov I, Andelman D and Orland H 2000 Adsorption of large ions from an electrolyte solution: a modified Poisson–Boltzmann equation Electrochimica Acta 46 221–9 [7] Boschitsch A H and Danilov P V 2012 Formulation of a new and simple nonuniform size-modified Poisson– Boltzmann description J. Comput. Chem. 33 1152–64 [8] Che J, Dzubiella J, Li B and McCammon J A 2008 Electrostatic free energy and its variations in implicit solvent models J. Phys. Chem. B 112 3058–69

2922

B Li et al

[9] Chu V B, Bai Y, Lipfert J, Herschlag D and Doniach S 2007 Evaluation of ion binding to DNA duplexes using a size-modified Poisson–Boltzmann theory Biophys. J. 93 3202–9 [10] Davis M E and McCammon J A 1990 Electrostatics in biomolecular structure and dynamics Chem. Rev. 90 509–21 [11] Eisenberg B, Hyon Y K and Liu C 2011 Energy variational analysis of ions in water and channels: field theory for primitive models of complex ionic fluids J. Chem. Phys. 133 104104 [12] Fogolari F, Brigo A and Molinari H 2002 The Poisson–Boltzmann equation for biomolecular electrostatics: a tool for structural biology J. Mol. Recognit. 15 377–92 [13] Gilbarg D and Trudinger N S 1998 Elliptic Partial Differential Equations of Second Order 2nd edn (Berlin: Springer) [14] Grochowski P and Trylska J 2008 Continuum molecular electrostatics, salt effects and counterion binding—A review of the Poisson–Boltzmann model and its modifications Biopolymers 89 93–113 [15] Grosberg A Y, Nguyen T T and Shklovskii B I 2002 the physics of charge inversion in chemical and biological systems Rev. Mod. Phys. 74 329–45 (Colloquium) [16] Hille B 2001 Ion Channels of Excitable Membranes 3rd edn (Sunderland, MA: Sinauer) [17] Howard J J, Perkyns J S and Pettitt B M 2010 The behavior of ions near a charged wall-dependence on ion size, concentration, and surface charge J. Phys. Chem. B 114 6074–83 [18] Kilic M S, Bazant M Z and Ajdari A 2007 Steric effects in the dynamics of electrolytes at large applied voltages: I. Double-layer charging Phys. Rev. E 75 021502 [19] Kilic M S, Bazant M Z and Ajdari A 2007 Steric effects in the dynamics of electrolytes at large applied voltages: II. Modified Poisson–Nernst–Planck equations Phys. Rev. E 75 021503 [20] Kralj-Igliˇc V and Igliˇc A 1996 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–91 [21] Lambert D, Leipply D, Shiman R and Draper D E 2009 The influence of monovalent cation size on the stability of RNA tertiary structures J. Mol Biol. 390 791–804 [22] Li B 2009 Continuum electrostatics for ionic solutions with nonuniform ionic sizes Nonlinearity 22 811–33 [23] Li B 2009 Minimization of electrostatic free energy and the Poisson–Boltzmann equation for molecular solvation with implicit solvent SIAM J. Math. Anal. 40 2536–66 [24] Lockless S W, Zhou M and MacKinnon R 2007 Structural and thermodynamic properties of selective ion binding in a K + channel PLoS Biol. 5 1079–88 [25] L´opez-Garcia J J and Jos´e Horno 2011 Poisson–Boltzmann description of the electrical double layer including ion size effects Langmuir 27 13970–80 [26] Lu B Z and Zhou Y C 2011 Poisson–Nernst–Planck equations for simulating biomolecular diffusion-reaction processes: II. Size effects on ionic distributions and diffusion-reaction rates Biophys J. 100 2475–85 [27] Maggs A C 2012 A minimizing principle for the Poisson–Boltzmann equation Europhys. Lett. 98 16012 ´ [28] Quesada-P´erez M, Mart´ın-Molina A and Hidalgo-Alvarez R 2004 Simulation of electric double layers with multivalent counterions: Ion size effect J. Chem. Phys. 121 8618–26 [29] Schoch R B, Han J and Renaud P 2008 Transport phenomena in nanofluidics Rev. Mod. Phys. 80 839–83 [30] Sharp K A and Honig B 1990 Electrostatic interactions in macromolecules: theory and applications Annu. Rev. Biophys. Biophys. Chem. 19 301–32 [31] Silalahi A R J, Boschitsch A H, Harris R C and Fenley M O 2010 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–9 [32] Strauss U P and Leung Y P 1965 Volume changes as a criterion for site binding of counterions by polyelectrolytes J. Am. Chem. Soc. 87 1476–80 [33] Tresset G 2008 Generalized Poisson–Fermi formalism for investigating size correlation effects with multiple ions Phys. Rev. E 78 061506 [34] Valisk´o M, Boda D and Gillespie D 2007 Selective adsorption of ions with different diameter and valence and highly charged interfaces J. Phys. Chem. C 111 15575–85 [35] Vlachy V 1999 Ionic effects beyond Poisson–Boltzmann theory Annu. Rev. Phys. Chem. 50 145–65 [36] V¨ortler H L, Sch¨afer K and Smith W R 2008 Simulation of chemical potentials and phase equilibria in two- and three-dimensional square-well fluids: finite size effects J. Phys. Chem. B 112 4656–61 [37] Wen J, Zhou S, Xu Z and Li B 2012 Competitive adsorption and ordered packing of counterions near highly charged surfaces: From mean-field theory to Monte Carlo simulations Phys. Rev. E 85 041406 [38] Zhou S, Wang Z and Li B 2011 Mean-field description of ionic size effects with non-uniform ionic sizes: a numerical approach Phys. Rev. E 84 021901