Mean-field electrostatics beyond the point-charge description
aa r X i v : . [ c ond - m a t . s o f t ] N ov Mean-field electrostatics beyond the point-charge description
Derek Frydel
School of Chemistry and Chemical Engineering and Institute of Natural Sciences,Shanghai Jiao Tong University, Shanghai 200240, ChinaandLaboratoire de Physico-Chime Theorique, ESPCI,CNRS Gulliver, 10 Rue Vauquelin, 75005 Paris, France (Dated: June 13, 2018)This review explores the number of mean-field constructions for ions whose structure goes beyondthe point-charge description, a representation used in the standard Poisson-Boltzmann equation.The exploration is motivated by a body of experimental work which indicates that ion-specificeffects play a significant role, where ions of the same valence charge but different size, polarizability,or shape yield quite different, and sometimes surprising results. Furthermore, there are many largeions encountered in soft-matter and biophysics that do not fit into a point-charge description, andtheir extension in space and shape must be taken into account of any reasonable representation.
PACS numbers:
I. INTRODUCTION
The present review explores the number of mean-field constructions for ions whose structure goes beyond thepoint-charge description, the representation used in the standard Poisson-Boltzman equation. The structural detailsneglected by a point-charge picture can be related to electrostatic structure of an ion, and they lead to polarizability,asymmetric interactions, or softening of Coulomb interactions if a distribution of an ion charge is extended in space,or they can be linked to the Pauli exclusion principle, and lead to the excluded volume effects or other softer types ofrepulsive interactions.The interest in incorporating these details came with the observation that charge is not the sole parameter thatdescribes an ion [1], as various ion-specific effects emerged in experimental work, beginning with the now famouswork of Hofmeister [2]. In the Hofmeister series ions of the same valance number can be arranged into a sequenceaccording to what effect they have on protein solubility. The larger and more polarizable ions weaken hydrophobicityof a protein, while the smaller ones strengthen it [3]. Fluid interfaces is another place where ion specificity plays animportant role in regulating a surface tension [4]. Dielectric decrement of an electrolyte upon addition of a salt is yetanother example where specificity of an ion can be measured [5, 6] and that is associated with how a given ion affectsa water structure.Another motivation for a model with structured particles is to come up with a more realistic description for water,which within the standard Poisson-Boltzmann equation is represented as a background dielectric constants [7]. Apolar solvent dissociates salts and dissolves ions, but electrostatic interactions between dissolved ions are not the sameas in vacuum. The solvent background screens interactions as dipoles of a solvent molecules align themselves alongfield lines coming from ion centers. This type of screening can be represented as an increased dielectric constant ofa background, which is about eighty times larger than that in vacuum. This description is accurate if alignment ofsolvent molecules is proportional to electrostatic field. But the proportionality relation breaks down when alignmentbecomes complete and solvent no longer responds to an external field, thus the screening effect slowly goes away.Finally, as a number of problems in soft-matter and biophysics increases and involves larger ions with significantsize and definite shape, a point-charge representation no longer serves as an accurate representation. Charge of anion has to be represented as extended in space in order to capture novel results that such ions give rise to [8].The mean-field approximation is a powerful method and yet very simple. It represents pair interactions as aneffective one-body external potential, w ( r ) = Z d r ′ ρ ( r ′ ) u ( r , r ′ ) , (1)which is a mean-potential a particle feels due to other particles in a system. u ( r , r ′ ) is the true pair interaction and ρ ( r ) is the equilibrium density of all particles. To construct a number density we write ρ ( r ) = ρ b e βw b e − β [ V ( r )+ w ( r )] , (2)where ρ b and w b is the value of ρ ( r ) and w ( r ) in a bulk, and V ( r ) is the external potential. The density now needsto be obtained self-consistently.The mean-field approximation includes particle interactions, and this is what sets it apart from the ideal-gas model,but it neglects correlations since it assumes that each particle interacts with a frozen distribution of particles thatdoes not change with the location of a test particle at r . The test particle in this sense is invisible: it measuresbut does not disturb the system. How accurate the mean-field is will depend on how negligible correlations are [9].Among theoretical systems for which the mean-field yields exact density profile are the hard-sphere system in theinfinite dimension limit [10], and the system of particles interacting via potential of the form U ( r ) = λ u ( λr ) in thelimit λ →
0, where u ( r ) is bounded and with finite range [11]. For physically relevant systems λ remains finite,and so the corrections to the mean-field are always present and real. An example of a system that is regarded asweakly-correlated liquid is the Gaussian core model at room temperature and so it is very accurately represented bythe mean-field approximation [12] .Correcting the mean-field with consecutive perturbative terms makes sense up to some value of an expansionparameter. At some point this no longer makes sense and electrostatics enters into the strong-coupling limit. Thissplits electrostatics into two ”worlds” [13]. In the strong-coupling ”wonder world” things stand up on their head: thesame-charged particles attract each other [14], for a counterion only system the distribution of counterions near acharged wall falls off exponentially as in the ideal-gas model [15–18]. The mean-field approximation has nothing tosay in this regime [19–21]. II. THE MEAN-FIELD APPROXIMATION
In this section we go over the mean-field approximation by starting from a complete partition function. To proceed,we postulate the system with scaled particle interactions, u λ ( r , r ′ ) = λu ( r , r ′ ), where λ = 0 corresponds to theideal-gas limit, and λ = 1 recovers the true system. The partition function for this system reads Z λ = 1 N !Λ N Z N Y i =1 d r i e − β P i Within the standard Poisson-Boltzmann equation ions are represented as point-charges. The model is constructedfrom the Poisson equation ǫ ∇ ψ ( r ) = − ρ c ( r ) , (17)which expresses the relation between the electrostatic potential ψ and the charge density ρ c ( r ) = K X i =1 q i ρ i ( r ) , (18)where the subscript i indicates an ion species, K is the number of species, and q i is the charge of an ion species i . Inthe Poisson equation ǫ is a background dielectric constant representing a solvent medium. The relation makes moresense if we transform it into somewhat different format, ψ ( r ) = Z d r ′ ρ c ( r ′ ) C ( r , r ′ ) , (19)where the Green’s function C ( r , r ′ ) denotes the functional form of Coulomb interactions, C ( r − r ′ ) = 14 πǫ | r − r ′ | , (20)and satisfies the fundamental equation ǫ ∇ C ( r − r ′ ) = − δ ( r − r ′ ) . (21)The Poisson equation seen through the format in Eq. (19) is nothing more than a restated definition of an electrostaticpotential, due to some distribution in space of Coulomb charges.Within this description different ionic species are distinguished by their valence number alone. Another approx-imation lies in the treatment of a solvent as a background dielectric constant ǫ . This description assumes that apolarization density of a polar solvent responds linearly to an electrostatic field, P = γ E , where E = − ∇ ψ is thelocal electrostatic field. The total dielectric constant then is ǫ = ǫ + γ , where ǫ is the dielectric constant in vacuum,and the action of a solvent is to screen electrostatic interactions by rising dielectric constant.To construct the Poisson-Boltzmann equation we need an expression for the charge density in terms of an electro-static potential. This can be obtained from the number densities, ρ i , obtained, in turn, from the mean-potential thatan ion of a species i feels due to other ions in a system, w i ( r ) = q i Z d r ′ ρ c ( r ′ ) C ( r , r ′ ) = q i ψ ( r ) , (22)which leads to the mean-field density, ρ i ( r ) = c i e − βq i ψ ( r ) , (23)where c i is the bulk concentration. Inserting this into the Poisson equation in Eq. (17) we arrive at the Poisson-Boltzman equation, ǫ ∇ ψ = − K X i =1 c i q i e − βq i ψ . (24)Later in the work, we refer to the Poisson-Boltzmann equation as the PB equation.For testing different mean-field models in this work we use the wall system where electrolyte is confined by a chargedwall with a surface charge σ c to a half-space x > 0. The PB equation reduces to 1D, ǫψ ′′ = − K X i =1 c i q i e − βq i ψ . (25)The neutrality condition, Z ∞ dx ρ c = − σ c , (26)fixes the boundary conditions at the location of a surface charge, − ǫψ ′ w = σ c , (27)where the subscript w indicates the value of a function at a contact with a wall.Finally, we derive the contact value theorem for the PB equation, which relates the value of a density at a wallcontact to bulk properties. To proceed we multiply the Poisson-Boltzmann equation by ψ ′ , ǫψ ′ ψ ′′ = − ψ ′ K X i =1 c i q i e − βq i ψ , (28)and rewrite it as ∂∂x (cid:18) ǫ ψ ′ (cid:19) = ∂∂x (cid:18) k B T K X i =1 c i e − βq i ψ (cid:19) . (29)The right-hand side term in parentheses is k B T ρ where ρ = P Ki =1 ρ i is the total number density. Integrating Eq. (29)from zero to infinity and using the boundary conditions, we find ρ w = ρ b + βσ c ǫ , (30)where ρ b = P Ki =1 c i . In the exact contact value theorem ρ b → P , where P is a bulk pressure. The present resultreflects the ideal-gas entropy of the Poisson-Boltzmann model where P = P id . B. Dipolar Poisson-Boltzmann equation Another possible point-particle is a point-dipole. The charge distribution of a dipole is comprised of two oppositecharges brought infinitesimally close to each other,lim q →∞ ε → h qδ ( r − r ′ ) − qδ ( r − r ′ + ε n ) i = − ( εq ) (cid:2) n · ∇ δ ( r − r ′ ) (cid:3) , (31)where n is the unit vector, and qε = p is the strength of a dipole moment. The limit q → ∞ is necessary to preventthe two charges from annihilation. Because of the limits, the distribution of a dipole can be represented as a gradientof a delta function.For the system of dipoles an additional stochastic degree of freedom comes out due to dipole orientation. A completeone particle distribution is, therefore, a function of a position and orientation, ̺ ( r , n ), which reduces to the numberdensity ρ i ( r ) = Z d p i ̺ i ( r , n ) . (32)A full distribution, however, is required to obtain a polarization density, P ( r ) = K X i =1 p i Z d n n ̺ i ( r , n ) . (33)To derive the Poisson equation for point-dipoles what is still needed is the formal expression of a charge density.Recalling that a dipole consists of two opposite charges ”glued” together, the charge density can be written as ρ c ( r ) = K X i =1 p i Z d n lim ε → " ̺ i ( r , n ) − ̺ i ( r + ε n , n ) ε = − K X i =1 p i Z d n h n · ∇ ̺ i ( r , n ) i = − ∇ · " K X i =1 p i Z d n n ̺ i ( r , n ) = − ∇ · P ( r ) . (34)The local charge density expressed as divergence of the polarization density can be understood as a charge transferfrom one volume element to another, and the non-zero ( ∇ · P ) implies that the charge that enters a given volumeelement is unbalanced by the charge that leaves it. The Poisson equation for the distribution of dipoles becomes ǫ ∇ ψ = ∇ · P . (35)To obtain the mean-field description of the present system, it is required to have an appropriate expression for P .As before, we start with the expression for a mean-potential. Knowing that an energy of a dipole in an external fieldis ( − p · E ), we write w i ( r , θ ) = p i · ∇ ψ = − p i | ∇ ψ | cos θ, (36)where θ is the angle between a local field E and a dipole p i . The corresponding mean-field distribution is ̺ i ( r , θ ) ∼ c i e βp i | ∇ ψ | cos θ , (37)and it reduces to the mean-field number density ρ i ( r ) ∼ c i Z π dθ sin θ e βp i | ∇ ψ | cos θ . (38)The properly normalized number density is ρ i ( r ) = c i sinh βp i | ∇ ψ | βp i | ∇ ψ | , (39)which reduces to a bulk density c i as a field vanishes. If dipoles were always aligned with a local field, the numberdensity would simply be ρ i → c i e βp i | ∇ ψ | . The different functional form in Eq. (39) reflects the fact that dipolesfluctuate around their preferred orientation.It remains now to obtain an expression for the polarization density, P ∼ K X i =1 p i c i Z d n n e − βp i ( n · ∇ ψ ) . (40)As the polarization vector is aligned with the field, P = P E E ! , (41)we write P as P = K X i =1 c i p i Z π dθ sin θ cos θe βp i | ∇ ψ | cos θ = K X i =1 p i c i sinh βp i | ∇ ψ | βp i | ∇ ψ | !" coth βp i | ∇ ψ | − βp i | ∇ ψ | , (42)or, P = − (cid:18) ∇ ψ | ∇ ψ | (cid:19) K X i =1 p i c i sinh βp i | ∇ ψ | βp i | ∇ ψ | ! L ( βp i | ∇ ψ | ) , (43)where L ( x ) = coth x − x (44)is the Langevin function which describes the degree of alignment of a dipole in a uniform electrostatic field. Anaverage dipole moment of a particle of a species i is given as h p i = p i L ( βp i | ∇ ψ | ) . (45)In the limit βp i | ∇ ψ | → L ≈ βp i | ∇ ψ | / 3. This is reasonable as an average dipole moment should be proportionalto an electrostatic field. But an alignment eventually must reach saturation and L cannot exceed 1, where L = 1signals perfect alignment. Note, however, that the limit L → L ≈ − / ( βp i | ∇ ψ | ).We have now everything that is needed for writing down a modified PB equation for a system of dipoles, ǫ ∇ ψ = − ∇ · " ∇ ψ K X i =1 p i c i sinh( βp i | ∇ ψ | ) βp i | ∇ ψ | L ( βp i | ∇ ψ | ) . (46)The dipolar PB equation was derived in [7], using the field-theory formalism, and was further explored in [22]. Themotivation was to treat water solvent more explicitly than as a background dielectric constant, the way it is donein the standard PB model. The standard PB equation assumes a linear and local relation between the polarizationdensity and the electrostatic field, P = γ E , so that the contributions of a polar solvent are subsumed into a dielectricconstant ǫ → ǫ + γ . The dipolar PB equation allows us to treat a polar solvent explicitly, ∇ · " ǫ + p c d sinh( βp | ∇ ψ | ) βp | ∇ ψ | L ( βp | ∇ ψ | ) ! ∇ ψ = − K X i =1 q i c i e − βq i ψ , (47)where p is the dipole moment of a solvent molecule. The linear relation between the polarization density and thefield are recovered in the limit βp E → P → (cid:18) βc d p (cid:19) | ∇ ψ | , (48)that recovers a space-independent dielectric constant, ǫ eff → ǫ + βc d p ǫ sol , (49)where ǫ sol denotes the dielectric constant of water, ǫ sol /ǫ = 80, and the parameters c d and p are tuned to accuratelyrecover this limit. Within this linear regime the dipolar PB equation behaves like its standard counterpart. For largervalues of βp | ∇ ψ | the linearity breaks down. Within the present model there are two sources of nonlinearity. Thefirst one lies within the Langevin function and captures the saturation of a polarization when a dipole aligns itselfalong a field. The second source of nonlinearity comes from the fact that point-dipoles are incompressible, ρ d ( r ) = c d sinh( βp | ∇ ψ | ) βp | ∇ ψ | , (50)and a local concentration can become arbitrarily large — a description somewhat unrealistic for water which is muchbetter represented as an incompressible fluid.For the wall geometry and symmetric electrolyte 1 : 1 the dipolar PB equation becomes, ǫ ψ ′′ + ∂∂x (cid:20) p c d sinh( βp ψ ′ ) βp ψ ′ L ( βp ψ ′ ) (cid:21) = 2 ec s sinh βeψ. (51)The boundary conditions at a wall are obtained, as they were for the standard PB equation, from the neutralitycondition, − ǫψ ′ = σ c + σ p , (52)where we introduce the polarization surface charge, σ p = p c d sinh( βp ψ ′ w ) βp ψ ′ w L ( βp ψ ′ w ) , (53)a polarization charge that accumulates at a wall due to charge transfer for a nonuniform polarization density. Thepolarization surface charge has always opposite sign to the bare surface charge σ c , and the surface charge can be saidto be screened.In Fig. (1) we plot results for the wall model. The dielectric constant is no longer uniform for the dipolar PBequation as the screening increases in the wall vicinity. Increased electrostatic screening reflects the excess of solvent x[nm] ε e ff / ε PBDPB x[nm] ρ − [ n m - ] PBDPB FIG. 1: The effective dielectric constant, ǫ eff = ǫ + p ρ d L /ψ ′ , and the counterion density, ρ − (the density of coions is denotedas ρ + ), for a wall model with surface charge σ c = 0 . − . The solvent parameters are: c d = 55 M and p = 4 . 78 D, suchthat in the linear polarization regime the dielectric constant of water is recovered, ǫ sol /ǫ = 1 + βc d p / ǫ = 80. The remainingparameters are λ B = βe / πǫ sol = 0 . 72 nm and c s = 0 . molecules near a wall. The saturation effect of the Langevin function is, therefore, not dominant. As a consequence,the counterion density is depleted from the wall region. (A dipole in a uniform electrostatic field adjusts its orientationbut not is position since there can be no gain in energy. In order to transport dipoles, a nonuniform field is required.This type of transport is referred to as dielectrophoresis. The reason for the concentration gradient of dipoles near awall in Fig. (1) reflects the fact that a field is nonuniform due to nonuniform distribution of counterions and coions).The dipolar PB model can be employed for studies of solvent mixtures. Given a mixture of two solvents withdifferent dipole moment, p = p , the solvent with higher polarity will prefer the vicinity of a charged surface, as amore efficient screening medium [1, 6, 22]. The dipolar PB model captures this behavior as seen in Fig. (2), wherethe hydration shell at a charged surface is comprised primarily of solvent of higher polarity. A heterogenous hydrationshell formed around ions dissolved in a solvent mixture can induce additional, ion-hydration interactions, leading stillto other ion-specific effects. x[nm] ρ d [ n m - ] p=4.78Dp=2.39D FIG. 2: Densities of two dipolar species in a solvent mixture. The same parameters as in Fig. (1), except the solvent parametersare: c d = c d = 27 . p = 2 p = 4 . 78 D. C. Langevin Poisson-Boltzmann equation In the dipolar PB equation solvent obeys ideal-gas entropy. But as water is not very compressible, a more realisticrepresentation of a solvent should take into account excluded volume interactions. Such interactions have beenimplemented through a local scheme based on the lattice-gas entropy [23–28]. Here, we make a simple assumptionthat water is incompressible, and the solvent density is uniform everywhere, ρ d ( r ) → c d . This reduction leadsto the Langevin PB equation, where the polarization density is determined solely by the Langevin function, P = c d p L ( βp | ∇ ψ | ) [29, 30], ∇ · " ǫ + c d p L ( βp | ∇ ψ | ) | ∇ ψ | ! ∇ ψ = − K X i =1 q i c i e − βq i ψ . (54)The effective dielectric constant in parentheses has two limiting behaviors. In the limit | ∇ ψ | → 0, it is the same asfor the dipolar PB equation, ǫ → ǫ + βc d p ǫ sol , (55)but now as | ∇ ψ | → ∞ and L → 1, the contributions of a solvent to dielectric response vanish, ǫ → ǫ + c d p | ∇ ψ | (56)and the nonlinear contributions of the Langevin model lead to dielectric decrement.Dielectric decrement has been observed for bulk electrolytes, and it reflects structural rearrangement of water dueto introduction of salt. For salt concentrations ranging between zero and 1 . ǫ eff ( c s ) = ǫ + αc s [5, 6, 31]. The rearrangement of water structure occursaround dissolved ions within the hydration shell. The orientation of these water dipoles is fixed by field lines originatingfrom ion centers, and they cannot respond to an external source of field. This behavior can be quantified with a crudemodel. As dipoles within the hydration shell are excluded from screening an external electrostatic field, the effectivedensity of free water dipoles becomes reduced, c d → c d − ( M + c + + M − c − ), where M ± is the solvation number ofwater molecules in a hydration shell around either a cation or anion. In the linear regime the dielectric constant ofwater is ǫ = ǫ + βc d p / 3. After addition of salt the effective concentration of salt is c d → c d − ( M + c + + M − c − ) andthe dielectric constant becomes ǫ → ǫ − ( c + M + + c − M − ) ( ǫ − ǫ ) c d . (57)In this simple picture M ± is salt specific.For the wall model and 1 : 1 symmetric electrolyte, the Langevin PB equation becomes, ǫ ψ ′′ + c d p L ′ ( βp ψ ′ ) = 2 ec s sinh βeψ. (58)The boundary conditions obtained from the neutrality condition is − ǫ ψ ′ w = σ c + σ p , (59)where σ p = c d p L ( βp ψ ′ w ) , (60)is the polarization surface charge.To obtain the contact value theorem, the Langevin PB equation is multiplied by ψ ′ , ǫ ψ ′ ψ ′′ = 2 ec s sinh( βeψ ) ψ ′ + c d p L ′ ( βp ψ ′ ) ψ ′ , (61)and after some manipulation we get, ∂∂x (cid:0) ǫ ψ ′ (cid:1) = ∂∂x (cid:0) c s β − cosh βeψ (cid:1) + ∂∂x (cid:18) c d p L ( βp ψ ′ ) ψ ′ (cid:19) − ∂∂x (cid:18) c d β − log (cid:20) sinh βp ψ ′ βp ψ ′ (cid:21)(cid:19) . (62)0After integration the contact value relation becomes ρ w = ρ b + β ǫ (cid:18) σ c − σ p (cid:19) − c d log (cid:20) sinh[ βp ( σ c + σ p ) /ǫ ] βp ( σ c + σ p ) /ǫ (cid:21) . (63)The results of the Langevin PB equation are shown in Fig. (3). The dielectric constant near a wall decreases asthe alignment of dipoles with field lines saturates, L → 1. This is an opposite trend to that found in the dipolar PBequation, which shows dielectric increment, see Fig. (1). The dielectric decrement of the Langevin model generatesstronger electrostatic interactions so that counterions stick more tightly to a charged wall. x[nm] ε e ff / ε PBLPB 0 0.1 0.2 0.3 0.4 x[nm] ρ − [ n m - ] PBLPB FIG. 3: The effective dielectric constant, ǫ eff = ǫ + p c d L /ψ ′ , and the counterion density, ρ − (the density of coins is denotedas ρ + ) for a wall model and the Langevin PB equation. The same parameters as in Fig. (1), except now the concentration ofpolar solvent is uniform and fixed at c d = 55 M. D. Point-dipoles with charge For the sake of illustration and as a way of transition to polarizable point-charges, we consider point-charges witha dipole moment. The mean-potential that an ion of a species i feels involves two parts, w i ( r , θ ) = q i ψ − p i | ∇ ψ | cos θ. (64)The corresponding mean-field distribution is ̺ i ( r , θ i ) ∼ c i e − β ( q i ψ − p i | ∇ ψ | cos θ i ) . (65)After integrating out the orientational degrees of freedom we arrive at the usual number density, ρ i ( r ) = c i e − βq i ψ ǫ sol sinh βp i | ∇ ψ | βp i | ∇ ψ | . (66)Before considering the polarization density we note that in the present model polarization density is associated withthe density of ions. In the Langevin model ions and dipoles were separate species. Following Eq. (42) we get P = (cid:18) ∇ ψ | ∇ ψ | (cid:19) K X i =1 p i c i e − βq i ψ sinh βp i | ∇ ψ | βp i | ∇ ψ | L ( βp i | ∇ ψ | ) , (67)and the mean-field Poisson equation becomes ∇ · " ǫ + K X i =1 p i c i e − βq i ψ sinh βp i | ∇ ψ | βp i | ∇ ψ | L ( βp i | ∇ ψ | ) ! ∇ ψ = − K X i =1 c i q i e − βq i ψ sinh βp i | ∇ ψ | βp i | ∇ ψ | . (68)1Here ions themselves contribute to the dielectric response of the medium on account of their inherent dipole moment.This leads to dielectric increment of a solution medium. Ions with permanent dipole moment are not common. Thetypical ions such as Cl − and Na + have spherical distributions. Ions with permanent dipole are encountered in ionicliquids, where ions tend to be larger molecular structures. E. Polarizable Poisson-Boltzmann equation A dipole moment of a polarizable ion is not permanent but is induced by an external electrostatic field accordingto the linear relation p = α E , (69)where α is the ion specific polarizability. Polarizability measures elasticity of an electron cloud of a molecule. Thelarger the electron cloud, the more deformable is the cloud. This behavior is manifest in the sequence for halide ions:F − < Cl − < Br − < I − .Polarizability is in fact a more general concept and measures response of an electron cloud to a time dependent fieldleading to a frequency dependent polarizability. Our strict concern is with static, or zero-frequency polarizability asvariations of an electric field induced by thermal fluctuations of an electrolyte operate at timescales much larger thantimescales of inner dynamics of an electron cloud. Frequency dependent polarizability can lead to other effects, suchas the London forces [32] when fluctuations in electron cloud of two nearby molecules become synchronized. These,however, are less important than the ion-dipole interactions generated by static polarizability [33, 34].Polarizability can be treated by recourse to a harmonic oscillator model, where opposite charges can be displacedrelative to one another under the action of an applied electric field, and the resorting force is proportional to adisplacement and the stiffness parameter k . The mean-potential for a species i is then written as w i ( r ) = q i ψ − p i | ∇ ψ | + k i d i , (70)where the last two terms characterize the energy of an induced dipole. The electrostatic energy of a dipole is alwaysin alignment with a field, and there is no orientational degree of freedom as for the case of a permanent dipole. Thelast term is the energy of a harmonic oscillator. To relate the stiffness parameter k to the polarizability α we startwith the Hooke’s law, F = k d , where the stretching force is electrostatic, F = q E , and a dipole moment is related toa displacement, p = q d ( q is the charge of the two charges being pulled apart). Substituting these definitions into theHooke’s law we get p = (cid:18) q k (cid:19) E , (71)so that from Eq. (69) we get k = q α . (72)The mean-potential in Eq. (70) can now be written as, w i ( r ) = q i ψ − α i | ∇ ψ | , (73)and the mean-field density becomes ρ i ( r ) = c i e − β ( q i ψ − α i | ∇ ψ | / , (74)leading to the following polarization density, P = − ∇ ψ K X i =1 α i c i e − β ( q i ψ − α i | ∇ ψ | / ! . (75)Polarization no longer depends on the Langevin function as it did for ions with a permanent dipole moment. Allnonlinearity of the expression is linked to the local ion density. The polarizable PB equation that results is [35, 36] ∇ · " ǫ + K X i =1 α i c i e − β ( q i ψ − α i | ∇ ψ | / ! ∇ ψ = K X i =1 q i c i e − β ( q i ψ − α i | ∇ ψ | / (76)2For a wall model and a symmetric 1 : 1 electrolyte, where all ions have the same polarizability α , the polarizablePB equation becomes ǫψ ′′ + ∂∂x h αc s ψ ′ cosh( βeψ ) e βαψ ′ / i = 2 ec s sinh( βeψ ) e βαψ ′ / . (77)The boundary conditions at the wall are − ǫψ ′ w = σ c + σ p , (78)where the polarization surface charge is σ p = 2 αc s ψ ′ w cosh( βeψ w ) e βαψ ′ w / . (79)Finally, the contact value theorem for the present model is ρ w = ρ b + β ǫ (cid:16) σ c − σ p (cid:17) , (80)obtained according to the procedure in Eq. (61). Equations (78), (79), and (80) can be combined to yield a singleequation for either ρ w or σ p . Below we write down the equation for the ratio σ p /σ c , which can be considered as ameasure of polarizability, (cid:18) σ p σ c (cid:19) + (cid:18) σ p σ c (cid:19) − ǫρ b βσ c + 2 ǫ βασ c !(cid:18) σ p σ c (cid:19) − ǫρ b βσ c ! = 0 . (81) σ p /σ c spans the range [0 , − 1) as α increases. σ p /σ c = − σ p cancels out the bare surface charge. Fromthe cubic equation above we obtain the dimensionless polarizability parameter, α ∗ = (cid:18) βσ c α ǫ (cid:19) , that controls the polarizability contributions at a charged surface. The other parameter, ǫρ b βσ c , depends on a saltconcentration.If we take an electrolyte at room temperature, the surface charge σ c = 0 . − , and the polarizability α/ (4 πǫ ) =10 ˚A , then the dimensionless polarizability parameter is α ∗ ≈ . 04. Polarizability α/ (4 πǫ ) = 10 ˚A correspondsroughly with the polarizability of iodide ion I − and is already rather high. We conclude then that polarizability oftypical salts has small effect on electrolytes.Polarizability contributions can be increased for dielectric media with low dielectric constant. Such a situation isrealized in ionic liquids, where the absence of a polar solvent permits unscreened electrostatic interactions, as ionicliquids are melted salts. In Fig. (4) we consider an electrolyte with reduced dielectric constant, ǫ/ǫ = 10 (whichyields a larger Bjerrum length, λ B = 5 . 76 nm). The increased dielectric constant near a wall region, which reflectsa counterion profile, ǫ eff = ǫ + 2 αc s cosh βeψ , generates a weaker attraction to a surface charge, so that counterionsbecome more spread out.The present model can by applied to study of ion specificity by assigning different polarizabilities to ion species. InFig. (5) we show density profiles for counterions with the same charge but different polarizability, α / (4 πǫ ) = 0˚A and α / (4 πǫ ) = 10˚A . Polarizable counterions being a better screening agent are preferred near a wall. 1. Negative Excess polarizability The present mean-field framework developed for polarizable ions has been used to capture the physics of dielectricdecrement caused by the restructuring of water as hydration shells form around dissolved ions [37, 38]. Since polariz-able ions in general cause dielectric increment, dielectric decrement can easily be realized when using negative valuesof polarizability. Negative polarizability occurs in quantum mechanics for molecules in excited state or for non-staticpolarizabilities, but in soft-matter it is an effective phenomena. Negative polarizability means that induced dipoleacts opposite to the local field. This conceptually captures the fact that the water dipoles in a hydration shell do notrespond to electrostatic field.3 x[nm] ε e ff / ε PBPPB x[nm] ρ − [ n m - ] PBPPB FIG. 4: Effective dielectric constant ǫ eff = ǫ + 2 αc s cosh βeψ , and the counterion density profile for reduced dielectric constant, ǫ/ǫ = 10 (in water ǫ/ǫ = 80). The relevant system parameters are: σ c = 0 . − , λ B = 5 . 76 nm, c s = 0 . α/ (4 πǫ ) = 10˚A . x[nm] ρ − [ n m - ] α α FIG. 5: The counterion profiles for polarizable and non-polarizable ions. The same parameters as in Fig. (4) but now only halfof the ions are polarizable: α = 0˚A and α = 10˚A . According to Eq. (57), the dielectric increment/decrement depends linearly on the salt concentration. The samelinear dependence is found for the polarizable PB model, ǫ eff − ǫ = α + ρ + + α − ρ − , (82)where, after comparing with Eq. (57), the negative excess polarizability can be approximated as α ± = − M ± (cid:18) ǫ − ǫ c d (cid:19) . (83)Even for a modest value of a solvation number, M ± = 4, the excess polarizability is already significant, α/ (4 πǫ ) ≈− , where we assume M + = M − = M and α + = α − = α .In Fig. (6) we show the results for α/ (4 πǫ ) = − . We compare the plots with positive polarizability of thesame magnitude, α/ (4 πǫ ) = 300˚A . The negative polarizability, as expected, lowers the dielectric constant near a4wall. Less intuitive is the fact that this leads to depletion of counterions from the wall region. The depletion is,furthermore, more significant than that for positive polarizability. If the electrostatic screening is reduced in the wallregion than counterions should hug to the wall more tightly. This is, at least, what we see in Fig. (3) for the LangevinPB model. So why does the Langevin PB model satisfies our intuitions and the polarizable PB model for negativepolarizabilities does not? Formal answer to this puzzle can be found by examining the contact value theorem inEq. (80). The sign of the polarization surface charge, σ p , does not matter, and any polarizability lowers the contactdensity, ρ w . x[nm] ε e ff / ε PBPPB, α>0 PPB, α<0 x[nm] ρ − [ n m - ] PBPPB, α>0 PPB, α<0 FIG. 6: Effective dielectric constant ǫ eff = ǫ +2 αc s cosh βeψ , and the counterion density profile for negative excess polarizability α/ (4 πǫ ) = − . The dielectric constant of a solvent background is that of water, ǫ/ǫ = 80. The remaining parametersare: σ c = 0 . − , λ B = 0 . 72 nm, c s = 0 . The two models, the Langevin and the polarizable PB equation with α < et al. [39] the Langevin PB equation with correlations hasbeen solved and it yields a non-monotonic density profile with a valey at a wall followed by a peak further away froma wall. The depletion is, therefore, captured but immediately at a wall and not for the entire profile. A possibleweak point of the negative polarizability model is the linearity assumption, ǫ eff − ǫ ∼ − c s . For homogenous solutionslinearity breaks down for higher concentrations [22], as hydration shells begin to overlap. This suggests that near awall, where concentrations are high, this too could have its effect. IV. FINITE-SPREAD POISSON-BOLTZMANN EQUATION An alternative approach to introduce structure of a charged particle is to smear its net charge within a finite volumeaccording to a desired distribution ω ( r − r ) such that q = Z d r ω ( r − r ) . (84)An arbitrary distribution is expected to depend on, in addition to the position r , the orientation characterized bythree angles. An arbitrary distribution embodies dipole, p = Z d r ( r − r ) ω ( r − r ) , (85)and higher order multipoles. If the two distributions at r i and r j , do not overlap, there is no difference betweenthe finite-spread and point-ion representation. The difference occurs for overlapping separations and the resulting5potential, U ( r i − r j ) = Z d r Z d r ′ ω ( r − r i ) ω ( r ′ − r j )4 πǫ | r − r ′ | , (86)is no longer described as a truncated series of multipoles.The finite-spread model does not want to provide a detailed electronic structure of an ion. This is beyond thescope of classical physics. But there are particles whose charge distribution is better described as extended in space,rather than as a sequence of multipoles. Among examples are charged rods, dumbbell shaped particles [8, 41–43], ormacromolecules whose non-electrostatic interactions are ”ultrasoft”, allowing interpenetration, and the distributionof charge in space is a sensible representation [44]. A perfect example is a polyelectrolyte in a good solvent whosecharges along a polymer chain appear on average as a smeared-out cloud due to quickly alternating configurations.Uncharged, two chain polymers interact via a Gaussian potential representing steric interactions of two self-avoidingpolymer chains [45]. Dendrimers offer another example of a soft, flexible macromolecule [46].There is also a more fundamental aspect of smeared-out charges: a smear-out point-charge is rid of divergence. Forsame-charged ions this eliminates effective excluded volume effects of a Coulomb potential and permits interpenetrationof two or more charges. For opposite-charged ions it leads to a new type of a Bjerrum pair where two ions collapse intoa neutral but polarizable entity [47–50]. The usual Bjerrum pair, formed between ions with hard-core interactions, isrepresented as a permanent dipole [51].Ultrasoft repulsive interactions (without the long-range Coulomb part) have been extensively studied, both for itstheoretical aspects and as a description of a soft matter system. Studies reveal two distinct behaviors. Some ultrasoftpotentials supports ”stacked” configurations, where two or more particles collapse, even though no true attractiveinteractions come into play [52, 53]. This behavior leads to a peak in a correlation function around r = 0. To this classof potentials belongs the penetrable sphere model [54]. The Gaussian core model [55], on the other hand, representsthe class of sort particles unable to support stacked configurations.As a note of interest, we mention that the interest in ultrasoft interactions is not confined to soft-matter. Thesoft-core boson model with interactions U ( r ) ∼ ( R + r ) − , where R is the soft-core radius, has been studied in[56] in connection to superfluidity. By removing the singularity from the potential, boson particles cluster and formcrystal with multiple particles occupying the same lattice sites. A. Spherical Distribution ω ( | r − r | ) In this section we consider ion species with spherically symmetric distribution ω i ( | r − r | ). As in previous mean-fieldconstructions, we start with the mean-potential that an ion of the species i feels, w i ( r ) = Z d r ′ ω i ( | r − r ′ | ) ψ ( r ′ ) . (87)The non-locality of the expression reflects the finite distribution of an ion charge in space and the fact that every partof this distribution interacts with an electrostatic field. The number density that follows is ρ i ( r ) = c i e − β R d r ′ ω i ( | r − r ′ | ) ψ ( r ′ ) . (88)To obtain the appropriate PB equation we need an expression for the charge density, which is given as the convolutionof the number density, ρ c ( r ) = K X i =1 Z d r ′ ω i ( | r − r ′ | ) ρ i ( r ′ ) . (89)Convolution is, again, a result of the finite extension of a charge. Using an explicit expression for ρ i we get ρ c ( r ) = K X i =1 c i Z d r ′ ω i ( r − r ′ ) e − β R d r ′′ ω i ( r ′ − r ′′ ) ψ ( r ′′ ) , (90)and the finite-spread PB equation for smeared-out ions is [57] ǫ ∇ ψ = − K X i =1 c i Z d r ′ ω i ( | r − r ′ | ) e − β R d r ′′ ω i ( | r ′ − r ′′ | ) ψ ( r ′′ ) . (91)6To complete the model, we still need to choose a specific spherical distribution. We model ions as uniformlydistributed charges within a spherical volume, ω i ( | r − r ′ | ) = 3 q i πR i θ ( R i − | r − r ′ | ) , (92)where q i and R i are the charge and the radius of an ion species i , respectively. The pair interaction between two ionswith charge q and size R , when the two ions overlap, is U ( r ≤ R ) = q πǫR " − (cid:16) rR (cid:17) + 316 (cid:16) rR (cid:17) − (cid:16) rR (cid:17) , (93)when not overlapping the usual Coulomb potential is recovered, U ( r > R ) = q πǫR (cid:16) rR (cid:17) − . (94)If overlap is complete the pair interaction remains finite, U (0) = 65 q πǫR , (95)and is said to be bounded. In Fig. (7) we plot various realizations of the pair potential U ( r ) for different ion size R .The degree of penetration clearly increases with increasing R . r/ λ B β U R/ λ B =0R/ λ B =0.1R/ λ B =0.2R/ λ B =0.4 FIG. 7: Pair potential between two charge distributions in Eq. (92) for different R . At overlapping separations the functionalform of a pair potential is that in Eq. (93). The mean-field Poisson equation for a symmetric 1 : 1 electrolyte, with ion distributions in Eq. (92), is ǫ ∇ ψ = 6 ec s πR Z d r ′ θ ( R − | r − r ′ | ) × sinh (cid:20) − βe πR Z d r ′′ θ ( R − | r ′ − r ′′ | ) ψ ( r ′′ ) (cid:21) . (96)For the wall model we the integral terms simplify, Z d r ′ θ ( R − | r − r ′ | ) f ( z ) = π Z R − R dz ′ f ( z + z ′ )( R − z ′ ) , (97)7where for f ( z ) = 1 we recover 4 πR / 3, a volume of a sphere. And the finite-spread PB equation becomes ǫψ ′′ = 6 ec s R Z R − R dz ′ ( R − z ′ ) sinh (cid:20) − βe R Z R − R dz ′′ ψ ( z + z ′ + z ′′ )( R − z ′′ ) (cid:21) . (98)For the wall model particle centers are confined to the half-space x > x > − R dueto finite size of an ionic charge as half of a sphere sticks out. The boundary conditions, therefore, are not determinedat the wall, x = 0, but at x = − R , − ǫψ ′ ( − R ) = σ c . (99)This implies that the surface charge is at x = − R . The contact value theorem, however, is not effected and is thesame as for the standard PB equation, ρ w = ρ b + βσ c ǫ , (100)where ρ w = ρ (0).In Fig. (8) we plot electrostatic quantities of penetrable ions: a charge density and an electrostatic potential. Unlikethe number density, these quantities are not confined to the region x > 0, and extend to x = − R as the charge ofan ion sticks out. Note how the sharp peak in the charge density for the standard PB model is smoothed-out in thefinite-spread model. x[nm] -3-2-10 ρ c [ n m - ] R=0R=0.8nm x[nm] e β ψ R=0R=0.8nm FIG. 8: The charge density and electrostatic potential for penetrable ions with charge distribution in Eq. (92) with R = 0 . x > x = 0 marks the half-space available to ioncenters. The results for R = 0 correspond to those for the standard PB equation. The system parameters are σ c = 0 . / m , λ B = 0 . 72 nm, and c s = 1 M. Fig. (9) shows number density profiles for penetrable ions. The first striking feature is that profiles are non-monotonic. More surprising still is the fact that a surface charged is overcharged: more counterions accumulate at awall than needed for neutralizing it. Overcharging is signaled by a peak in coion density and indicates attraction ofcoions to a same-charged surface. Attraction between same-charged plates was found for dumbbell shaped counterionsin [8], suggesting that charge inversion is a common feature of charges extended in space.A closer look into plots reveals that overcharging and consequent charge inversion is a more complex phenomenon.To magnify these features we plot in Fig. (10) the number densities for a symmetric 3 : 3 electrolyte. What we seeis not a simple charge inversion but rather an alternating layers of counterions and coions leading to oscillations indensity profiles. This behavior is reminiscent of polyelectrolyte layer-by-layer adsorption onto a charged substrate[58–60].In Fig. (11) we show Monte Carlo snapshots for counterions adsorbed onto a charged wall. The figure compares twosystems: counterions that approximate a point-charges with size R = 0 . R = 0 . x[nm] ρ − [ n m - ] R=0R=0.8nm x[nm] ρ + [ n m - ] R=0R=0.8nm FIG. 9: The number density profiles for counter- and co-ions for penetrable ions near a charged wall. The same parameters asin Fig. (8). x[nm] ρ [ n m - ] counterionscoions FIG. 10: Number density profiles for counterions and coions with valance number 3. The increased electrostatic interactionsmagnify the features in Fig. (9). Otherwise the same parameters as those in Fig. (8). observed for this system, although deviations from the mean-field are sufficiently significant to yield counterion densityprofiles different than those of the standard PB equation by being shifted closer to a wall. On the other hand, theconfiguration for R = 0 . R , see Eq. (93). By taking the limit R → ∞ we recoverthe ideal gas limit. A divergence in the pair interactions for point-ions acts as an effective hard-shpere interactionwhose radius depends on the Bjerrum length as well as a surface charge. This leads to the effective excluded volume9 FIG. 11: Monte Carlo configuration snapshots for counterions adsorbed onto a charged wall (within the slice 0 < x < . 35 nm).The conditions are the same as for results in Fig. (12). The circles representing particles have diameter σ = 0 . R = 0 . R = 0 . ρ d = 2 . − and ρ d =2 . − , respectively. For comparison, the surface charge density is σ c /e = 2 . 50 nm − , indicating overcharging for R = 0 . interactions. By removing a divergence and permitting interpenetration, the effective excluded volume interactionsare eliminated. This completely changes the lateral structure of adsorbed counterions. 1. stacked configurations A configuration snapshot for penetrable ions in Fig. (11) gives impression that counterions form stacked config-urations, detected from the correlation function by the presence of a peak around r = 0. The presence of stackedformations is, furthermore, linked to instability of the Kirkwood analysis [11, 61–64]. The Kirkwood analysis is themean-field type of an analysis. As already said in section II, even though the mean-field approximation does not incor-porate explicit correlations in its free energy formulation, the correlations can be extracted using exact thermodynamicrelations. This is possible because the mean-field is not a self-consistent approximation. Thus, the exact relation inEq. (12) yields Eq. (14) from the mean-field density, which then can be put in the form of the Ornstein-Zernikeequation, as done in Eq. (15), which, when Fourier transformed, becomes, h ( k ) = − βu ( k )1 + ρβu ( k ) , (101)and the corresponding structure factor, defined as S ( k ) = 1 + ρh ( k ), is S ( k ) = 11 + ρβu ( k ) . (102)All is good as long as u ( k ) is a non-negative function. But if for some modes k the Fourier transformed pair potentialis negative, S ( k ) becomes divergent for some wave number k , indicating divergent fluctuations. In addition, theonset of this so-called Kirkwood instability coincides with an onset of the long-range correlation order and with abifurcation point where a constant density no longer yields minimum free energy and another periodic solution takesprecedence [11]. The instability was later linked to the spinodal of the supercooled liquid.0 x[nm] ρ + [ n m - ] mean-fieldsimulation FIG. 12: The coion density profile near a charged wall. The mean-field theory very accurately reproduces the exact resultsof the Monte Carlo simulation. The number of particles in the simulation box is N + = N − = 1200, and the box size is L y = L z = 16nm and L x = 12nm. The periodic boundary conditions are in the lateral ( y, z )-directions. The other parametersare: λ B = 0 . 72 nm, R = 0 . σ c = 0 . − . Simulations of the penetrable sphere model (exhibiting the Kirkwood instability) showed the existence of stackedconfigurations (referred to as ”clumps” in that work) at temperatures below instability [65]. Individual stacks arrangedinto crystal structure (corresponding to the global minimum) or amorphous glassy structures (corresponding to a localminimum). Stacked structures gave rise to a peak around r = 0 in the correlation function. A more thorough analysissupported by simulations linked the instability to crystals with multiply occupied lattice sites [52, 66, 67]. Potentialswhose Fourier transformed potential u ( k ) is positive remain always stable and are not found to form stacked crystalformations. Instead their solid phase exhibits reentrant melting upon squeezing, a behavior seen in water [68].In the present work we are interested in the liquid structure before the onset of instability. In particular, we wantto know if the penetrable ions exhibit Kirkwood instability linked to stacked formations, and if yes, what role theyplay in a charge inversion mechanism.It is enough to consider the one component plasma of penetrable ions. The Fourier transformed pair potentialdepends on the distribution ω and is obtained from Eq. (86), U ( k ) = ω ( k ) ǫk . (103)For point-ions ω ( k ) = q , and for penetrable ions lim k → ω ( k ) = q , since at large separations the usual Coulombinteractions are recovered. The difference between point and spread-out ions is seen for large k , where S ( k ) reflectsbehavior for small separations. Regardless of the distribution ω ( k ), U ( k ) ≥ k and the Kirkwood instabilitydoes not occur. Stacked formations, therefore, can be excluded as playing any part in a mechanism for charge inversion.It is not clear, however, whether the conclusion holds for all Coulomb potentials with soft-core, or if it is specificto smeared-out ions. To address this concern, we consider the following soft-core Coulomb potential, u ( r ) = ( q πǫr if r ≥ σ q πǫσ if r < σ, (104)whose Fourier transform is u ( k ) = q ǫk sin kσkσ . (105) U ( k ) is no longer a non-negative function and yields instability in the mean-field structure factor. The mean-fieldcorrelation function, furthermore, shows a peak at r = 0, h (0) = − π Z ∞ dk βU ( k ) k ρβu ( k ) , (106)1under certain conditions. We conclude that while the smearing-out procedure cannot lead to Kirkwood instability,this behavior is not general to all Coulomb potentials with soft-core.. B. Needle-ions: the case for non-spherical ω ( r − r ′ ) The modified PB equation for non-spherical distributions is more complicated as there are three additional degreesof freedom for a particle orientation to cope with. Levy et al. [22] derived the modified PB equation for a generaldistribution using the field-theory methodology.In this section we consider, as an example of non-spherical distribution, needle-ions. We perform our constructionas before, by writing down a mean-potential from which we obtain a number and charge densities. A needle-ionconsists of a charge q uniformly distributed along a line of length d . The charge distribution of a needle-ion is ω ( r − r , n ) = qd Z d/ − d/ ds δ ( r + s n − r ) , (107)where r is a midpoint and n is a unit vector that designates orientation. The mean-potential that a needle-ion of aspecies i feels when its center is at r is w i ( r , n ) = q i d i Z d r ′ ψ ( r ′ ) Z d i / − d i / ds δ ( r + s n − r ′ ) , (108)or, suppressing the delta function, we may alternatively write w i ( r , n ) = q i d i Z d i / − d i / ds ψ ( r + s n ) . (109)A nonlocal contribution comes from particle’s finite extension in space. The mean-field expression for the numberdensity in space and orientation is ̺ i ( r , n ) ∼ c i exp " − βq i d i Z d i / − d i / ds ψ ( r + s n ) . (110)We still need an expression for a charge density to complete the construction. A charge density at location r hasnonlocal contributions from neighboring ions that lie within a spherical region of radius d/ 2. This region is describedby the Heaviside step function θ ( d/ − | r − r ′ | ). However, being located within this region is not sufficient conditionfor contributing to the charge density at r . There is additional condition of orientation: only ions with orientation n = r ′ − r | r ′ − r | , (111)contribute to the charge density at r . Each species’ contribution to the charge density is ρ ic ( r ) ∼ q i Z d r ′ θ (cid:18) d i − | r ′ − r | (cid:19) ̺ i (cid:18) r ′ , r ′ − r | r ′ − r | (cid:19) . (112)Using Eq. (110) to substitute for ̺ i , the total charge density becomes ρ c ( r ) = N X i =1 q i c i πd i ! Z d r ′ θ (cid:18) d i − | r ′ − r | (cid:19) × exp " − βq i d i Z d i / − d i / ds ψ (cid:18) r ′ + s r ′ − r | r ′ − r | (cid:19) , (113)where the coefficient 6 / ( πd i ) comes from the limit ψ → 0, where all orientations are equally probable and the chargedensity properly recovers its bulk value, ρ c → P Ki =1 c i q i . For a symmetric 1 : 1 electrolyte and ions of the same length2we get ρ c ( r ) = 12 ec s πd Z d r ′ θ (cid:18) d − | r ′ − r | (cid:19) × sinh " − βed Z d/ − d/ ds ψ (cid:18) r ′ + s r ′ − r | r ′ − r | (cid:19) . (114)By inserting this result into the Poisson equation, ǫ ∇ = − ρ c , we obtain the desired modified PB equation forneedle-ions, − ǫ ∇ ψ = 12 ec s πd Z d r ′ θ (cid:18) d − | r ′ − r | (cid:19) × sinh " − βed Z d/ − d/ ds ψ (cid:18) r ′ + s r ′ − r | r ′ − r | (cid:19) . (115) C. Dumbbell ions There are cases when multivalent organic ions, such as certain DNA condensing agents or short stiff polyeletrolyes,have a rod-like structure wherein charges are spatially separated from each other [8, 40–43]. These separations arenot small and are comparable to typical screening lengths, ∼ d . The distribution of a dumbbell ion located at r is ω ( r − r , n ) = qδ ( r − r ) + qδ ( r − r − d n ) . (116)These dumbbell counterions were found to give rise to attraction between two same-charged plates by bridging twosurfaces, thereby providing a finite equilibrium distance between surfaces.The mean-field construction for dumbbell ions follows the usual route. The mean-potential for an ion species i is w i ( r ) = q i ψ ( r ) + q i ψ ( r + d i n ) , (117)which then leads to the following distribution, ̺ i ( r , n ) ∼ c i e − βq i ψ ( r ) − βq i ψ ( r + d i n ) . (118)A properly normalized charge density then is written as ρ c ( r ) = K X i =1 q i c i e − βq i ψ ( r ) Z d r ′ δ ( d i − | r − r ′ | )4 πd i e − βq i ψ ( r ′ ) , (119)and it remains now to put this into the Poisson equation. The mean-field Poisson equation for a 1 : 1 dumbbellelectrolyte with equal sized particles becomes ǫ ∇ ψ ( r ) = 4 ec s Z d r ′ δ ( d − | r − r ′ | )4 πd sinh (cid:20) βeψ ( r ) + βeψ ( r ′ ) (cid:21) . (120)For the wall model this becomes ǫψ ′′ = 2 ec s d Z d − d ds sinh (cid:2) βeψ ( x ) + βeψ ( x + s ) (cid:3) θ ( x ) θ ( x + s ) , (121)where the Heaviside step functions ensure that dumbbell ions do not go through the wall.3 V. SHORT-RANGE NON-ELECTROSTATIC INTERACTIONS So far we have considered only interactions due to electrostatic structure of an ion. In addition to these there are alsonon-electrostatic interactions, generally short-ranged and repulsive, the most obvious of which are the excluded volumeinteractions whose source is traced to the Pauli exclusion principle which prohibits two electrons from occupying thesame quantum state [72].For ions in aqueous solution the excluded volume interactions are enhanced due to formation of a hydration shell.Excluded interactions can, furthermore, lead to effective, softer type of interactions. For example, the effectiveinteractions between two linear polymers in a good solvent can be represented with the Gaussian functional formand result from a self-avoiding walk between dissolved polymer chains [45]. There are many other types of exoticinteractions in soft-matter systems, and some further examples include effective interactions between star polymers,dendrimers, etc. [44]. In this section we consider different schemes for incorporating short-range interactions ofnon-electrostatic origin. A. the mean-field implementation The simplest way to implement short-range interactions is to use the mean-field framework. Considering point-charges, the mean-potential for such an implementation is βw i ( r ) = βq i ψ ( r ) + β K X j =1 Z d r ′ ρ j ( r ′ ) u ij ( r − r ′ ) , (122)where u ij designates the non-electrostatic interactions between particles of the species i and j . The correspondingmean-field density is ρ i = c i e − βq i ψ e − β P Kj =1 R d r ′ ( ρ j ( r ′ ) − c j ) u ij ( r − r ′ ) , (123)which recovers bulk density in the limit ψ → 0. The charge density is ρ c = P Ki =1 q i ρ i and the resulting mean-fieldPoisson equation is − ε ∇ ψ = K X i =1 c i q i e − βq i ψ e − β P j R d r ′ ( ρ j ( r ′ ) − c j ) u ij ( r − r ′ ) . (124)The approximation consists of two coupled equations, Eq. (123) and (124). Note that the implementation of non-electrostatic interactions leads to nonlocal approximation where the density is convoluted with the pair interaction u ij ( r − r ′ ).If the electrolyte is symmetric, 1 : 1, and there is only one type of short-range interactions for each particles, Eq.(124) reduces to a more familiar form, ε ∇ ψ ( r ) = 2 c s sinh( βeψ ) e − β R d r ′ ( ρ ( r ′ ) − c s ) u ( r − r ′ ) , (125)where the total number density, ρ = ρ + + ρ − , is given by ρ = 2 c s cosh( βeψ ) e − β R d r ′ ( ρ ( r ′ ) − c s ) u ( r − r ′ ) . (126)As a specific example, we consider penetrable sphere ions (PSM) whose short-range repulsive interaction is βu ( | r − r ′ | ) = ε θ ( σ − | r − r ′ | ) , (127)where σ is the diameter of a penetrable sphere, and ε is the strength. In the limit ε → ∞ the hard-core interactionsare recovered.At first we consider uncharged penetrable spheres and compare results with those from simulation. Results fora wall model are plotted in Fig. (13) which shows density profiles of penetrable spheres near a planar wall. Themean-field becomes less accurate as ε increases where it overestimates the contact values, ρ w , which are related tothe bulk pressure via the contact value theorem, ρ w = βP. (128)4 x/ σ ρπ σ / ε=1ε=0.5ε=2 PSM simulationmean-field FIG. 13: Density profiles of uncharged penetrable spheres near a planar wall. Particle centers are confined in the x -axis, x ∈ [0 , . σ ]. The number of particles is fixed, R dx ρ ( x ) = N = 1000. For Monte Carlo simulation the dimensions of thesimulation box are 12 . × . × . σ . The box encloses N = 1000 particles. In the y and z directions periodic boundaryconditions are used. The influence of the second wall is minor, and we refer to this system as the wall model. Overestimated contact density values imply that the mean-field pressure is larger than the true one, P mf > P . Toobtain the mean-field pressure we use the virial equation [73], βP = ρ b − π ρ b Z ∞ dr r g ( r ) ∂βu ( r ) ∂r , (129)from which we discard correlations, g ( r ) = 1, according to the mean-field procedure, and we get βP mf = ρ b + ερ b (cid:16) πσ (cid:17) , (130)where we used ∂βu ( r ) ∂r = − εδ ( r − σ ). The mean-field contact value theorem for penetrable spheres, therefore, is ρ w = ρ b + ερ b (cid:16) πσ (cid:17) . (131)The lower contact density for a true system implies a neglect of correlations in the mean-field approximation.Having in mind hard-spheres as a model system for excluded volume interactions, we can make contact with it bysetting ε = 1, where the resulting mean-field pressure, βP mf = ρ b (1 + 4 η ) , (132)agrees to the second virial term with the the pressure for hard-spheres, where η = πρσ / ε = 1 and solve Eq. (125) and Eq. (126) for symmetric 1 : 1electrolyte. For the wall model the boundary conditions are the same as for the standard PB equation, and thecontact value theorem is ρ w = ρ b (1 + 4 η ) + βσ c ǫ . (133)In Fig. (14) we plot density profiles for counterions near a charged wall. In comparison with the standard PB equation,the penetrable sphere ions generate a non-monotonic structure of a double-layer, where we see the emergence of asecondary peak. The structure is a result of overcrowding, where counterions coming to neutralize the surface chargecannot be packed too closely together. The simulation results for hard-sphere particles with the same diameter yielda profile with stronger structure. The penetrable sphere model captures only qualitatively these features.The present model can be used to study ion specific effects. For neutral surfaces, size asymmetry can lead todifferent density profiles of ions with the same valance number. This leads to charge build-up across an interface.5 x[nm] ρ − [ n m - ] PBPB+MFsimulation FIG. 14: The counterion density near a charged wall. The system is confined between two charged walls at x = 0 and x = 6 nm and the surface charge is σ c = − . − . The other parameters are: the Bjerrum length λ B = 0 . 72 nm, and thediameter of penetrable spheres σ = 0 . R dx ρ i ( x ) = N i = 300 and N + = N − . The penetrability parameter is set to ε = 1. For simulation we used the hard-spherelimit, ε → ∞ . In Fig. (15) we show density profiles near a neutral wall for a 1 : 1 electrolyte with size asymmetry. The largercations exhibit greater structure and are squeezed against the wall which, in turn, leads to a charge build-up thatpulls anions. Simulation results for hard-spheres with the same diameters show similar profiles, however, anions havegreater structure as they are depleted from the immediate wall vicinity. x[nm] ρ [ n m - ] cationsanionssimulation x[nm] ρ c [ n m - ] FIG. 15: Density profiles of cations, anions, and of total charge near an uncharged wall. Two parallel uncharged plates at x = 0 and x = 6 nm confine all particle centers. The dielectric constant is the same across an interface. The Bjerrum lenght is λ B = 0 . 72 nm, the ion sizes are σ ++ = 0 . σ −− = 0 nm, and σ + − = 0 . R dxρ i ( x ) = N i where N + = N − = 300.The simulation was done for the same system but for hard-spheres, ε → ∞ , while in the numerical model we used ε = 1. B. short-range interactions beyond the mean-field In this section we develop a more accurate implementation of short-range interactions, while keeping electrostaticsat the mean-field level. Such a procedure introduces asymmetric treatment of different parts of a pair potential. To6formally set up and justify this asymmetry of methods, we consider the scaled pair potential, u λij = u hs ij + λ q i q j πǫ | r − r ′ | ! , (134)where u hs ij is the hard-sphere potential, and λ is the scaling parameter. For λ = 0 a hard-sphere system is recovered.The density of all ions is independent of λ and is kept fixed by the external electrostatic potential, ψ λ ext . The partitionfunction for this system is Z λ = 1 Q Kj =1 ( N j !Λ N j ) Z N Y i =1 d r i e − β P Ni Q i ψ λ ext ( r i ) e − β P Ni,j (cid:2) λ QiQj πǫ | r i − r j | + u hs ij ( r i , r j ) (cid:3) , (135)where N = P Kj =1 N j is the total number of particles, N j is the number of particles of a species i , and the charges { Q i } have the following values Q N j − +1 = Q N j − +2 = · · · = · · · = Q N j − + N j = q j , where q i is the charge of a species i . The exact functional form of ψ λ ext ( r ) is not needed and it will not appear in the final result. It is sufficient to knowthat it keeps densities fixed at their physical shape for any value λ , and ψ λ =1ext recovers the true external potential, ψ ext .The free energy, βF = − log Z , is obtained from thermodynamic integration, F [ { ρ i } ] = F λ =0 [ { ρ i } ] + Z dλ ∂F λ ∂λ = F id [ { ρ i } ] + F hsex [ { ρ i } ] + Z d r ρ c ( r ) ψ λ =0ext ( r ) + Z dλ ∂F λ ∂λ , (136)where F hsex [ ρ ] is the excess free energy due to hard-sphere interactions and is a functional of density only. The integrandof the last term after evaluation is ∂F λ ∂λ = Z d r ρ c ( r ) ∂ψ λ ext ( r ) ∂λ + 12 K X i,j Z d r Z d r ′ q i ρ i ( r ) q j ρ j ( r ′ )4 πǫ | r − r ′ | g λij ( r , r ′ ) , (137)where g ij = 1 + h ij . Inserting this into Eq. (136) we get F [ { ρ i } ] = F id [ { ρ i } ] + F hsex [ { ρ i } ] + Z d r ρ c ( r ) ψ ext ( r ) + 12 Z d r Z d r ′ ρ c ( r ) ρ c ( r ′ )4 πǫ | r − r ′ | + 12 K X i,j Z d r Z d r ′ q i ρ i ( r ) q j ρ j ( r ′ )4 πǫ | r − r ′ | Z dλ h λij ( r , r ′ ) . (138)By setting correlations to zero, h λij = 0, F [ { ρ i } ] ≈ F id [ { ρ i } ] + F hsex [ { ρ i } ] + Z d r ρ c ( r ) ψ ext ( r ) + 12 Z d r Z d r ′ ρ c ( r ) ρ c ( r ′ )4 πǫ | r − r ′ | , (139)we have an approximation that completely neglects terms coupling the electrostatic and hard-core interactions, andthat consists of the free energy for hard-spheres plus the mean-field electrostatic correction.Densities are obtained from the minimum condition, δFδρ i ( r ) = 0, and ρ i ( r ) = c i exp h − βq i ψ ( r ) − δβF hsex δρ i ( r ) + βµ ex i , (140)where µ ex = ∂F hsex ∂ρ | ρ = ρ b is the excess chemical potential over the ideal contribution, µ id = log ρ b Λ , and ψ is the totalelectrostatic potential. Implementing this into the Poisson equation we get − ǫ ∇ ψ = K X i =1 c i q i exp h − βez i ψ − δβF hsex δρ i ( r ) + βµ ex i . (141)7Then for 1 : 1 electrolyte with all ions having the same size, ǫ ∇ ψ = 2 c s e sinh( βeψ ) exp h − βδF hsex δρ + βµ ex i (142)where ρ = ρ + + ρ − is ρ = 2 c s cosh( βeψ ) exp h − βδF hsex δρ + βµ ex i . (143) 1. perturbative expansion and the dilute limit To complete the approximation it remains to find an expression for the excess free energy. The first two terms ofthe virial expansion for F hsex are [74] F ex = 12 Z d r Z d r ρ ( r ) ρ ( r ) ¯ f ( r )+ 16 Z d r Z d r Z d r ρ ( r ) ρ ( r ) ρ ( r ) ¯ f ( r ) ¯ f ( r ) ¯ f ( r )+ . . . (144)where ¯ f ( r ) = 1 − e − βu ( r ) is the negative Mayer f -function, and for hard-spheres is given by the Heaviside function,¯ f ( r ) = θ ( σ − r ). In the dilute limit the first term dominates and constitutes an accurate approximation,lim ρ → F hsex = 12 Z d r Z d r ′ ρ ( r ′ ) ρ ( r ′ ) θ ( σ − | r − r ′ | ) , (145)which yields δF hsex δρ ( r ) = Z d r ′ ρ ( r ′ ) θ ( σ − | r − r ′ | ) , (146)and the number density becomes ρ ± ( r ) = c b e ∓ βeψ ( r ) − R d r ′ ( ρ ( r ′ ) − c b ) θ ( σ −| r − r ′ | ) . (147)Incidentally, the dilute limit approximation is the same as the mean-field implementation of the penetrable sphereinteractions with ε = 1, as both approximations are designed to give the lowest order term of the virial expansion for F ex (see Fig. (14) and Fig. (15) for performance of the dilute limit approximation). 2. nonperturbative approach Further expansion of the excess free energy does not constitute an efficient scheme. Already the second lowest terminvolves the three-body overlap contributions that numerically is difficult to deal with. A more powerful approach is anonperturbative scheme. A nonperturbative construction keeps numerical complexity of the dilute limit approximationbut incorporates additional terms (generally an infinite set of terms) that lead to accurate behavior for some limitingcondition.One example is the weighted density approximation of Ref.[75, 76]. This approximation is constructed in terms ofweighted density which constitutes a building block of the theory and is suggested from the lowest order term of thevirial series for F hsex , ¯ ρ ( r ) = 18 Z d r ′ ρ ( r ′ ) θ ( σ − | r − r ′ | ) . (148)¯ ρ is dimensionless and normalized to recover the packing fraction in a bulk, ¯ ρ b = η = πσ ρ b / 6. The approximationassumes that F hsex has a general form F hsex = Z d r ρ ( r ) φ ex (¯ ρ ( r )) , (149)8where φ ex denotes an excess free energy per particle and is a function of ¯ ρ ( r ). Tarazona suggested a generalizedCarnahan-Starling approach, where for φ ex he used the quasi-exact Carnahan-Starling equation [75], φ csex (¯ ρ ( r )) = ¯ ρ ( r )(4 − ρ ( r ))(1 − ¯ ρ ( r )) , (150)but defined as a function of a weighted density. Now, in addition to recovering the dilute limit exactly, the approxi-mation recovers the homogenous limit. Furthermore, the construction satisfies the contact value theorem, ρ w = P cs ,where P cs is the Carnahan-Starling expression for hard-sphere pressure. The excess chemical potential of this gener-alized Carnahan-Starling approach is δF hsex δρ ( r ) = f csex (¯ ρ ( r )) + 18 Z d r ′ ρ ( r ′ ) ∂f csex ∂ ¯ ρ ( r ′ ) θ ( σ − | r − r ′ | ) , (151)and the densities of ionic species are obtained from Eq. (140).A nonperturbative construction can be further improved by increasing the number of weighted densities as buildingblocks of the theory. Some improvements where implemented as a result of careful studies of the direct correlationfunction, which suggested a density dependent weight function [77]. The breakthrough approach, however, came withthe Rosenfeld’s fundamental measure theory [78]. Motivated (at least in part) by desire to construct a theory thatrecovers the 1D limit behavior (a property later referred to as the dimensional crossover), Rosenfeld obtained a newset of weight functions by decomposing the Heaviside step function, θ ( σ ij − | r i − r j | ) = ω i ⊗ ω j + ω i ⊗ ω j + ω i ⊗ ω j + ω i ⊗ ω j − ω i ⊗ ω j − ω i ⊗ ω j , (152)where ω α ⊗ ω β = Z d r ′ ω iα ( r ′ − r i ) ω jβ ( r ′ − r j ) , (153)and the relevant weight functions are ω i ( r ) = θ ( R i − r ) ,ω i ( r ) = δ ( R i − r ) , ω i ( r ) = r r δ ( R i − r ) , and ω i ( r ) = ω i ( r ) / (4 πR i ), ω i ( r ) = ω i ( r ) / (4 πR i ), and ω i ( r ) = ω i ( r ) / (4 πR i ). The six weighted densities that resultare, n α = K X i =1 Z d r ′ ρ i ( r ′ ) ω iα ( r − r ′ ) , (154)and a general formula for the excess free energy is βF hsex = Z d r Φ RF ( { n α ( r ) } ) . (155)Based on the scaled particle theory results [79–81], Rosenfeld came up with the following functional form [78, 82–84],Φ RF = − n log(1 − n ) + n n − n · n − n + n − n ( n · n )24 π (1 − n ) . (156)The construction also recovers the PY direct correlation function for homogenous liquids.In Fig. (16) we plot the counterion density profiles for a symmetrical electrolyte 1 : 1 confined between twoparallel hard walls. The conditions are the same as in Fig. (14). The WDA gives improvements over the dilute limitapproximation, but the DFT fundamental measure theory results agree most closely with the simulation. In Fig.(17) we plot density profiles for a 1 : 1 electrolyte with size asymmetry confined by uncharged walls. The conditionsare the same as in Fig. (15). The DFT here is less accurate, although the charge density profile quite well agrees withsimulation. The disparity can be traced to the lack of correlations in the mean-field treatment of electrostatics, whichbecome important for an uncharged wall system. The presence of correlations is best seen in the contact density,related to the pressure via the contact value theorem, which is lower in the simulation results, and which indicatesnegative correlational contributions due to formation of Bjerrum pairs [51].9 x[nm] ρ − [ n m - ] DFTWDAdilute limitsimulation FIG. 16: The counterion density near a charged wall for the DFT scheme. Conditions as in Fig. (14). DFT denotes the densityfunctional theory based on the fundamental measure theory [78], and WDA denotes the weighted density approximation basedon the generalized Carnahan-Starling equation [75]. x[nm] ρ [ n m - ] DFTdilute limitsimulation x[nm] ρ c [ n m - ] FIG. 17: Density profiles near an uncharged wall for various approximations. Conditions as in Fig. (15). 3. Correlations The results in Fig. (17) for neutral confinement indicate that despite of highly accurate expression for hard-coreinteractions, F hsex , the mean-field treatment of electrostatics is not sufficient and the correlations play a dominant role.The neglected correlational contribution to the free energy, taken out of the complete expression in Eq. (138), is F c = 12 K X i,j Z d r Z d r ′ q i ρ i ( r ) q j ρ j ( r ′ )4 πǫ | r − r ′ | Z dλ h λij ( r , r ′ ) . (157)As λ → h λij does not vanish, but instead h λij → h hs ij .Applied to homogenous electrolytes, the formula in Eq. (157) yields the charging process formula, f c = 12 K X i =1 q i c i Z dλ " π Z ∞ dr r K X j =1 q j c j h λij ( r )4 πǫr , (158)0where f c = F c /V . By identifying the term ρ λc,i ( r ) = K X j =1 q j c j h λij ( r ) (159)as the charge distribution around an ion of the species i fixed at the origin (constituting the charge correlation hole),the term in brackets becomes an electrostatic potential that a test ion of the species i feels due to surrounding ionsin the system, f c = 12 K X i =1 q i c i Z dλ ψ λi , (160)where the superscript λ indicates that the interactions are scaled. The formula is the expression of the ”chargingprocess”, a common route for obtaining the correlational free energy. Substituting for ψ λ the linear Debye-H¨uckelsolution leads to the Debye charging process, which for 1 : 1 electrolyte of ions of the same size becomes [20], q i ψ λi = − e πǫ " κ √ λ κσ √ λ , (161)and βf c = − πσ " log( κσ + 1) − κσ + ( κσ ) . (162)This correlation term constitutes a weak-coupling correction and it does not capture the formation of Bjerrum pairs[51]. It gives, however, some estimate of what the contributions of neglected correlations are.Application of the charging process to inhomogeneous electrolytes is not easy as it requires the precise functionalform for ψ λ ext ( r ) which ensures that densities remain constant throughout charging. The implementation of coupledcontributions of hard-core and electrostatic contributions remains a challenge. There are some perturbative extensionsto the DFT theory based on the reference fluid density [85, 86] and which address this issue. C. local schemes After reviewing nonlocal approximations for hard-sphere interactions, it may seem a regression to discuss next localapproximations. Nonlocal construction based on weighted densities captures discrete structure of a fluid and is foundto satisfy the contact value theorem sum rule. A local construction, on the other hand, is expressed in terms of localdensity (or a weighted density with a delta weight function), and as such, it does not posses discrete structure of aliquid that is implicit in weight functions, and fails to satisfy the contact value theorem [87]. In other words, densityprofiles produced by the local type of an approximation are unphysical.Then why even bother with local approximations? The first answer is simplicity. But this does not justify a modelas a description of the world. A more reasonable justification may sound like this. For true electrolytes as they arefound in laboratories the exact nature of the excluded volume interactions is not known with precision and there aremany different and complex contributions. The hard-sphere model is itself an idealization. The local approximation,in spite of its shortcomings, can offer a first glance and an estimate of excluded volume effects. One, however, hasto know how to interpret such a local approximation. The structureless density profile cannot be read as physical.The saturation effect of a local density triggered by overcrowding is an artifact of the model. But although densityis unphysical, it does not mean that every other quantity that follows is equally so. For example, the incorrectcontact density does not imply an incorrect contact potential. In fact, the contact potential values were found tobe reasonably well estimated by a local scheme [87]. The local saturation of a density profile, its flattening near acharged surface, captures qualitatively the fact that a double-layer is elongated due to excluded volume effects, andthis in turn reproduces, at least qualitatively, the increase in electrostatic potential that is less efficiently screened.Within local approximation the excess free energy is F hsex = Z d r f ex ( ρ ( r )) , (163)1where the excess free energy density f ex is a function of a local density, and a functional derivative becomes classicalderivative, δF hsex δρ ( r ) = ∂f ex ∂ρ ( r ) = µ ex ( ρ ( r )) , (164)which yields the following density ρ i ( r ) = c i e − β (cid:2) q i ψ ( r )+ µ ex ( ρ ( r )) − µ ex ( ρ b ) (cid:3) , (165)The mean-field Poisson equation becomes, ǫ ∇ ψ ( r ) = − K X i =1 q i c i e − β (cid:2) q i ψ ( r )+ µ ex ( ρ ( r )) − µ ex ( ρ b ) (cid:3) . (166)We assume that all diameters are the same.To complete the model, it remains to choose expression for µ ex . There are several equations of state we can choosefrom. The hard-sphere model (or the quasi-exact Carnahan-Starling equation) is βPρ = 1 + η + η − η (1 − η ) −→ βµ ex = 8 η − η + 3 η (1 − η ) . (167)A cruder van der Waals model for excluded volume interactions is βPρ = 11 − νρ −→ βµ ex = − log(1 − νρ ) + νρ − νρ , (168)where ν denotes the excluded volume. Finally, the lattice-gas model is [24] βPρ = − log(1 − η ) η −→ βµ ex = − log(1 − η ) . (169)In Fig. (18) we compare the lattice-gas equation of state with that for hard-spheres. The two models are completelydifferent. There is no agreement in any limit. The lattice-gas curve is relatively flat and then exhibits a sharp rise as η → η β P / ρ Carnahan-Starlinglattice-gas FIG. 18: The equation of state βP/ρ as a function of a packing fraction η = πσ / e − βµ ex = 1 − η . The result is intuitive and expresses the fraction of an availablevolume not taken up by other particles. This simple result leads to the following density ρ i = c i e − βq i ψ (cid:18) − ν P Ki =1 ρ i − ν P Ki =1 c i (cid:19) , (170)where ν = πσ / ρ i = c i e − βq i ψ ν P Ki =1 c i ( e − βq i ψ − . (171)In the limit ψ → 0, the standard Poisson-Boltzmann equation is recovered, ρ i → c i e − βq i ψ . But if potential becomeslarge a density cannot increase indefinitely as it is bounded from above, ρ ≤ ν − . The modified Poisson-Boltzmannequation that results is ǫ ∇ ψ = − P Ki =1 q i c i e − βq i ψ ν P Ki =1 c i ( e − βq i ψ − . (172)Specializing to the 1 : 1 electrolyte we get [24] ǫ ∇ ψ = − ec s sinh βeψ νc s (cosh βeψ − 1) (173)This is the modified Poisson-Boltzmann equation as derived in [24]. It yields the same boundary condition as thestandard PB equation. Also, the model does not lead to the true contact value theorem, ρ w = P + βσ c ǫ , where inplace of P we use the lattice-gas pressure in Eq. (169). Instead, it obeys another contact value relation, ρ w = 1 ν (cid:20) − (1 − η b ) e − βσ c ν/ ǫ (cid:21) , (174)since, as was said before, the density is not physical. The model introduces a new length scale, νσ c , that correspondsto the width of counterion layer that would form if all counterions were allowed to come to a charged surface and theexceeded volume effect was the only interaction. Now, even for the vanishing screening length, κ − → 0, a double-layer will have thickness νσ c (where κ = √ πc s λ B is the Debye screening parameter). Based on these two competinglength scales it is possible to estimate the importance of the excluded volume effects. If νσ c > κ − , we should expectthe excluded volume effects to play a significant role.In Fig. (19) we compare the results of the modified PB equation with other approximations. The density profile ofthe modified PB equation shows unphysical saturation of a local density. The standard PB equation, in fact, yieldsbetter agreement with the DFT near a wall including a contact density. The modified PB equation, however, yieldsbetter results for electrostatic potential which comes close to the DFT approximation at a wall contact. The modifiedPB equation captures the fact that the surface charge is less efficient screening when the excluded volume interactionsare involved. In Fig. (20) we plot contact potential as a function of the surface charge. The modified PB equationcaptures the influence of the excluded volume effects in relation to the standard PB equation. As σ c becomes large,the agreement with the DFT is less perfect.The question still remains how using the Carnahan-Starling equation of state for the local approximation wouldchange the results. After all, when choosing the lattice-gas equation of state, the only criterion that was followedwas simplicity. It turns out that the equation of state for hard-spheres gives worse agreement with the DFT and itexaggerates overcrowding by yielding too high contact potentials. It somewhat seems a stroke of luck that the lattice-gas equation of state provides both simplicity and relative accuracy, suggesting that some cancellation of errors isbeing involved. VI. CONCLUSION The present review provides a framework for constructing various mean-field models for ions with some sort ofstructure and provides a number of modified PB equations to which this construction leads. All possibilities, ofcourse, cannot be exhausted, but with a large number of detailed constructions it should not be difficult to build a3 x[nm] ρ − [ n m - ] DFTPBMPBsimulation x[nm] e β ψ DFTPBMPB FIG. 19: Counterion density and potential profiles for the conditions as in Fig. (14) and Fig. (14). σ c e β ψ w PBMPBDFT FIG. 20: Electrostatic potential at a wall contact as a function of a surface charge. The system parameters are : λ B = 0 . 72 nm, σ = 0 . c s = 0 . model that fits a given situation. One possible direction to pursue further is to explore in more detail models forions with finite charge distribution, all the way until making contact with polyelectrolytes, whose single configurationis represented with a Brownian walk type of a distribution. What was also left out from the review were modelsthat combine several structures together. For example, a spherical charge distribution could be supplemented withrepulsive Gaussian interactions representing a self-avoiding walk of two polymer chains. This would render a morerealistic representation of polyelectrolytes. But again, such combinations are not difficult to infer from providedconstructions. It is, in fact, one of the goals of this review to motivate such new constructions.The review also leaves some suggestions for the future work. As there is an ever increasing number of new macro-molecules with interactions ranging from ultrasoft to hard-core, particles whose shape is not fixed but flexible, thereis an ever growing demand for their accurate representation. One possible way to proceed is to explore the presentmean-field framework and implement some sort of elastic behavior to allow charge distributions to deform into mostoptimal shape, so that near an interface particles and their mutual interactions are modified. Similar elastic behaviorcould be supplemented to non-electrostatic type of interactions.As far as the treatment of short-range non-electrostatic interactions is concerned, the mean-field is sufficient ifinteractions are soft. The handling of hard-sphere interactions, however, is far more challenging. The most efficienttheory for hard-core interactions, the fundamental measure DFT, shows shortcomings even for the weak-coupling limitconditions. To construct a more accurate theory it is, therefore, necessary to incorporate correlations, thus, to go4beyond the conveniences of the mean-field. But for charged hard-core particles correlations are difficult to implementas they couple hard-core and electrostatic interactions. The treatment of these coupled, short- and long-range,interactions is, in fact, one of the most outstanding problems of soft-matter electrostatics.The idea of a modified PB equation, of representing more physics and more accurately through a model based ona single differential equation, has caught some momentum, and there are models that go beyond the mean-field andattempt to implement effects seen in the intermediate- and strong-coupling limit. An exemplary case is the work byBazant et al. [88, 89] where the authors suggested the modified PB equation with dielectric constant represented asa linear differential operator. The model aims to represent ionic liquids were opposite ions strongly associate andcoexist as a Bjerrum pair rather than as free ions.There is an additional motivation for pursuing various mean-field constructions. Simple models such as chargedhard-spheres are easy to simulate, and for these systems one could simply stick to simulations to cover the entire rangeof electrostatics, from weak- to strong-coupling regime. But there are systems that are not easy to simulate. This isespecially true for polarizable ions and for explicit treatment of water. For these cases the mean-field constructionprovides a real alternative, sometimes the only choice.Finally, it should be reminded and stressed one more time what the boundaries of the mean-field treatment are andwhat type of electrostatics it is capable of representing. Not only it neglects correlational corrections already active inthe weak-coupling regime, but it completely fails in the strong-coupling limit as a predictive tool. The strong-couplinglimit electrostatics calls for different treatments. ACKNOWLEDGMENT The author would like to thank Tony Maggs as well as the ESPCI lab, were bulk of this work was being done, forhelpful and friendly working atmosphere. This work was supported in part by the agence nationale de la recherchevia the project FSCF. [1] D. Ben-Yaakov, D. Andelman and D. Harries, J. Phys. Chem. B , 6001 (2009).[2] F. Hofmeister, Arch. Exp. Pathol. Pharmacol. , 247 (1888).[3] Y. Zhang, P. S. Cremer, Curr. Opin. Chem. Biol. , 658 (2006).[4] A. P. dos Santos, A. Diehl and Y. Levin, Langmuir , 10778 (2010).[5] J. B. Hasted, D. M. Ritson, and C. H. Collie, J. Chem. Phys. , 1 (1948).[6] D. Ben-Yaakov, D. Andelman and R. Podgornik, J. Chem. Phys. , 074705 (2011).[7] A. Abrashkin, D. Andelman, and H. Orland, Phys. Rev. Lett. , 077801 (2007).[8] Y. W. Kim, Y. Yi, and P. A. Pincus, Phys. Rev. Lett. , 208305 (2008).[9] R.R. Netz, H. Orland, Eur. Phys. J. E , 203 (2000).[10] W. Klein and H. L. Frisch, J. Chem. Phys. , 968 (1986).[11] N. Grewe and W. Klein, J. Math. Phys. , 1735 (1977).[12] A. A. Louis, P. G. Bolhuis, and J. P. Hansen, Phys. Rev. E , 7961 (2000).[13] A. Naji, M. Kanduˇc, R. R. Netz, R. Podgornik, ”Exotic electrostatics: unusual features of electrostic interactions betweenmacroions”, in Understanding Soft Condensed Matter via Modeling and Computation (Series in soft condensed matter) Vol. 3, 2010, pp. 265-295.[14] P. Linse and V. Lobaskin, Phys. Rev. Lett. , 4208 (1999).[15] V. I. Perel and B. I. Shklovskii, Physica A , 446 (1999).[16] B. I. Shklovskii, Phys. Rev. E , 5802 (1999).[17] A. G. Moreira and R. R. Netz, Europhys. Lett. , 705 (2000).[18] L. ˇSamaj and E. Trizac, Phys. Rev. Lett. , 078301 (2011).[19] I. Rouzina and V. A. Bloomfield, J. Phys. Chem. , 9977 (1996).[20] Y. Levin, Rep. Prog. Phys. , 1577 (2002).[21] A. Yu. Grosberg, T. T. Nguyen, and B. I. Shklovskii, Rev. Mod. Phys. , 329 (2002)[22] A. Levy, D. Andelman, H. Orland, J. Chem. Phys. , 164909 (2013).[23] J. J Bikerman, Philos. Mag. 384 (1942).[24] I. Borukhov, D. Andelman, and H. Orland, Phys. Rev. Lett. , 435 (1997).[25] C. Azuara, H. Orland, M. Bon, P. Koehl, and M. Delarue Biophysical Journal , 5587 (2008).[26] D. H. Mengistu, K. Bohinc, and S. May, Europhys. Lett. Phys. Rev. Lett. , 087801 (2009).[28] A. Igliˇc, E. Gongadze, and K. Bohinc Bioelectrochemistry , 223 (2010).[29] D. Frydel and M. Oettel, Phys Chem Chem Phys. , 4109 (2011). [30] D. Frydel and M. Oettel, ”Extended Poisson-Boltzmann descriptions of the electrostatic double layer: implications forcharged particles at interfaces”, in ”New challenges in Electrostatics of Soft and Disordered Matter”, Eds. J. Dobnikar, A.Naji, D.Dean and R. Podgornik, Singapore (2014).[31] A. Levy, D. Andelman, H. Orland, Phys. Rev. Lett. , 227801 (2012).[32] F. London, Trans. Farad. Soc. , 8 (1937).[33] R. R. Netz, J. Phys.: Condens. Matter , S2353 (2004).[34] R. R. Netz, Curr. Opin. Colloid Interface Sci. , 192 (2004).[35] L. B. Bhuiyan and C. W. Outhwaite, J. Phys. Chem. J. Chem. Phys. , 234704 (2011).[37] M. M. Hatlo, R. van Roij and L. Lue, EPL J. Chem. Phys. , 174903 (2012).[39] M. Ma and Z. Xu, http://arxiv.org/pdf/1410.4661[40] K. Bohinc, A. Igliˇc and S. May, Europhys. Lett. , 494 (2004).[41] S. May, A. Iglic, J. Rescic, S. Maset and K. Bohinc, J. Phys. Chem. B , 1685 (2008).[42] K. Bohinc, J. Rescic, J. Maset and S. May, J. Chem. Phys. , 07411 (2011).[43] K. Bohinc, J. M. A. Grime, L. Lue, Soft Matter , 5679 (2012).[44] C. N. Likos, Phys. Rep. , 267 (2001).[45] A. A. Louis, P. G. Bolhuis, J.-P. Hansen, and E. J. Meijer, Phys. Rev. Lett. , 2522 (2000).[46] C. N. Likos, M. Schmidt, and H. L¨owen, M. Ballauff and D. P¨otschke, Macromolecules , 2914 (2001).[47] D. Coslovich, J.-P. Hansen and G. Kahl, Soft Matter , 1690 (2011).[48] D. Coslovich, J.-P. Hansen and G. Kahl, J. Chem. Phys. , 244514 (2011).[49] A. Nikoubashman, J.-P. Hansen, and G. Kahl, J. Chem. Phys. , 094905 (2012).[50] P. B. Warren and A. J. Masters, J. Chem. Phys. , 074901 (2013).[51] M. E. Fisher and Y. Levin, Phys. Rev. Lett. , 3826 (1993).[52] C. N. Likos, A. Lang, M. Watzlawek, and H. L¨owen, Phys. Rev. E , 031206 (2001).[53] B. M. Mladek, D. Gottwald, G. Kahl, M. Neumann, and C. N. Likos, Phys. Rev. Lett. , 045701 (2006).[54] C. Marquest and T. A. Witten, J. Phys. France , 1267 (1989).[55] F. H. Stillinger, J. Chem. Phys. , 3968 (1976).[56] F. Cinti, T. Macr´ı, W. Lechner, G. Pupillo and T. Pohl, Nature Communications , 3235 (2014).[57] D. Frydel, Y. Levin, J. Chem. Phys. , 174901 (2013).[58] I. Borukhov, Physica A , 315 (1998).[59] J. Hemmerle, V. Roucoules, G. Fleith, M. Nardin, V. Ball, Ph. Lavalle, P. Marie, J.-C. Voegel, P. Schaaf, Langmuir ,10328 (2005).[60] D. Mijares, M. Gaitan, B. Polk, D. DeVoe, J. Res. NIST , 61 (2010).[61] N. Grewe and W. Klein, J. Math. Phys. , 1729 (1977).[62] W. Klein and N. Grewe, J. Chem. Phys. , 5456 (1980).[63] W. Kunkin and H. L. Frisch, J. Chem. Phys. , 181 (1969).[64] T. Naitoh and K. Nagai, J. Stat. Phys. , 391 (1974).[65] W. Klein, H. Gould, R. A. Ramos, I. Clejan and A. I. Melcuk, Physica A , 738 (1994).[66] C. N. Likos, M. Watzlawek, and H. L¨owen, Phys. Rev. E , 3135 (1998).[67] M. Schmidt, J. Phys.: Condens. Matter , 10163 (1999).[68] F. H. Stillinger and D. K. Stillinger, Physica A , 358 (1997).[69] J. Urbanija, K. Bohinc, A. Bellen, S. Maset, A. Igliˇc, J. Chem. Phys. , 105101 (2008).[70] M. Kanduˇc, A. Naji, R. Podgornik, J. Phys. Condens. Matter. , 424103 (2009).[71] R. I. Slavchov and T. I. Ivanov, J. Chem. Phys. , 074503 (2014).[72] F. J. Dyson and A. Lenard, J. Math. Phys. , 423 (1967); F. J. Dyson and A. Lenard, J. Math. Phys. , 698 (1968); F. J.Dyson, J. Math. Phys. , 1538 (1967).[73] J. P. Hansen and I.R. McDonald, Theory of Simple Liquids , 2nd ed. (Academic Press, London, 1986).[74] R. Evans, Adv. Phys. A , 143 (1979).[75] P. Tarazona, Mol. Phys. , 81 (1984).[76] P. Tarazona and R. Evans, Mol. Phys. 847 (1984).[77] R. Evans, in Fundamentals of Inhomogeneous Fluids, edited by D. Henderson, Chap. 3 (Dekker, New York, 1992), p. 85.[78] Y. Rosenfeld, Phys. Rev. Lett. , 980 (1989).[79] H. Reiss, H. L. Frisch, J. L. Lebowitz, J. Chem. Phys , 369 (1959).[80] M. Heying and D. S. Corti, J. Phys. Chem. B , 19756 (2004).[81] F. H. Stillinger, P. G. Debenedetti and S. Chatterjee, J. Chem. Phys. , 204504 (2006).[82] P. Tarazona, J. A. Cuesta, and Y. Martinez-Raton, Lect. Notes Phys. 247 (2008).[83] R. Evans, Lecture Notes at 3rd Warsaw School of Statistical Physics (Warsaw University Press, Kazimierz Dolny, 2009)pp. 43-85.[84] R. Roth, J. Phys.: Condens. Matter , 063102 (2010).[85] Y. Rosenfeld, J. Chem. Phys. , 8126 (1993).[86] D. Gillespie, W. Nonner, and R. S. Eisenberg, Phys. Rev. E , 031503 (2003).[87] D. Frydel, Y. Levin, J. Chem. Phys. , 164703 (2012).[88] M. Z. Bazant, B. D. Storey, and A. A. Kornyshev, Phys. Rev. Lett. , 046102 (2011). [89] B. D. Storey and M. Z. Bazant, Phys. Rev. E , 056303 (2012)., 056303 (2012).