The effect of ionic size on polyion-small ion distributions in a cylindrical double layer

C.N.Patra 1,3 , L.B.Bhuiyan 2 1 Theoretical Chemistry Section, RC & CD Division, Chemistry Group, Bhabha Atomic Research Centre, Mumbai 400 085, India 2 Laboratory of Theoretical Physics, Department of Physics, Box 23343, University of Puerto Rico, San Juan, Puerto Rico 00931–3343 3 Department of Materials Science and Engineering University of Utah 122 S. Central Campus Drive Salt Lake City, Utah 84112–0560


Introduction
The study of the interactions of a polyion with the ionic atmosphere surrounding it is quite important in diverse areas such as colloidal systems [1], industrial polyelectrolytes [2], and biologically significant systems like DNA [3][4][5].The structure and thermodynamics of all these systems manifest large salt dependencies.A detailed knowledge of the spatial distribution of small ions in the vicinity of a polyelectrolyte is therefore fundamental to a microscopic understanding of these polyelectrolyte solutions.
The polyion together with its charge cloud constitute what is known as an electric double layer in the literature with the geometry of the polyion shaping the double layer geometry.In the planar symmetric situation one has the well known planar double layer (PDL), which has been extensively studied over the past three decades using formal statistical mechanical methods (see for example, reviews [6][7][8][9]).In contrast, the double layer in other geometries, for instance the cylindrical double layer (CDL) and the spherical double layer (SDL), has received comparatively less attention.
Theoretical advances in the last few decades have made possible the calculation of the ionic distribution and mean electrostatic potential profiles, and consequently the characterization of various other thermodynamic properties of the polyion -small ion systems.In almost all of these studies, the polyion is modelled as an infinitely long, rigid charged cylinder, while the small, mobile ions are treated as charged hard spheres in a continuum dielectric.The situation when the polyion is at infinite dilution, viz., the primitive model cylindrical double layer system, and the situation when the polyion has a finite concentration, viz., the cylindrical cell model, have both been treated.The theoretical works are mainly based on four distinct approaches.The first one is the counterion condensation (CC) theory of Manning (see for example, reference [10]), which has been quite extensively applied in several biological problems [11] including analyses of competitive binding equilibria [12].The second approach corresponds to methods based on the classical Poisson-Boltzmann (PB) theory which has been applied in different forms, viz.linear [13], nonlinear [14,15], and the cell model [16,17].The third approach is based on more rigorous formal statistical mechanical methods which includes functional expansions [18], liquid structure integral equations [19][20][21][22][23][24][25][26], the modified Poisson-Boltzmann theory (MPB) (including the cell model) [27][28][29][30][31][32][33], and the density functional theory (DFT) [34,35], all of which have been found to be successful in varying degrees in predicting the structural and dynamical information of such systems.The fourth one involves Monte Carlo (MC) simulations (again including the cell model) [30,32,33,[36][37][38][39][40][41][42], which gives interesting insights into the exact energetics, density distributions, and colligative properties of such systems.The qualitative and quantitative comparisons of these approaches have been attempted from time to time (see for example, references [24,40,43]), but it is beyond the scope of the present paper.We intend to concentrate more on some of the specific issues such as the effect of ionic size on double layer structure, involving the application of the DFT and MPB theories to the cylindrical double layer.
The DFT is a relatively new analytical approach to the electric double layer phenomenon and has been found to be useful in some earlier structural studies of the PDL [44][45][46][47].In recent years these techniques have been utilized in more exhaustive studies of the PDL system by Boda et al. [48][49][50][51] and Patra and Ghosh [52], while Yu et al. [53] have applied the methods to the SDL.The studies have confirmed the usefulness of the density functional approach for these systems.Such detailed studies for the CDL, however, lacking in the literature.An earlier DFT study of the CDL was limited, viz., the zeta potentials were not calculated [34,35].Furthermore, it is tempting to see how the approach compares with earlier and more established theories of the CDL such as the MPB, which has already seen successful applications to the cylindrical geometry [27][28][29].The DFT technique applied in this work is a partially perturbative one [34,54] where the hard sphere contribution to the first order correlation function is evaluated from the weighted density approach (WDA) of Denton and Ashcroft [55] and the electrical contribution is approximated by a quadratic expansion with respect to the corresponding bulk electrolyte [56].Henceforth in this paper the term "density functional theory" or its acronym DFT, will refer to this version.
The MPB theory belongs to the family of potential based theories in the statistical mechanics of charged fluids.It has evolved from the attempts to improve upon the classical PB theory by incorporating (i) inter-ionic correlations through a fluctuation potential term, and (ii) ionic exclusion volume effects -features that are missing in the classical formulation.
For a cylindrical double layer situation the MPB theory was first derived by Outhwaite in the mid-1980s [27].Subsequently the theory was applied to a primitive model CDL by Bhuiyan and Outhwaite [28,29].In these applications the value of the ionic diameter was held fixed at d = 4.25 • 10 −10 m.It is important to point out that with few exceptions, most of the double layer studies in planar, cylindrical, and spherical symmetries have been at this value of the ionic diameter.The very recent MC and DFT works by Boda et al. [50,51] on PDL have for the first time examined the dependency of double layer properties on ionic size with diameters d = 2 • 10 −10 and 3•10 −10 m being used besides the habitual d = 4.25•10 −10 m.An interesting outcome is that the classical Gouy-Chapman-Stern theory of the PDL [57][58][59] becomes increasingly less accurate as the ionic size decreases.In a later detailed comparative study by Bhuiyan and Outhwaite [60] of the DFT results versus the corresponding MPB predictions vis-a-vis the MC, the two theories were found to be remarkably consistent over a broad range of physical parameters including different ionic sizes.
In this paper we will explore the DFT and MPB theories for a primitive model CDL with a particular emphasis on the effect of variation of ionic size on structural properties of the double layer.It is further noted that no comparison of the DFT and the MPB theory for the CDL exists in the polyelectrolyte literature.

Theoretical formulation 2.1. Molecular model
The polyion is modelled as an infinite, isolated, rigid cylinder bearing a uniform axial charge density given by where e is the proton charge, b is the length per unit charge (inverse of linear charge density), β = (k B T ) −1 , k B is the Boltzmann's constant, T is the temperature, and is the dielectric constant of the pure solvent modelled as a uniform dielectric continuum, which is taken to be = 78.358(characterizing water).Throughout this work we set T = 298.15K, ξ = 4.2, and b = 1.7 • 10 −10 m, which are the accepted values for a double-stranded DNA [43].The small ions (with α denoting the species) are modelled using the restricted primitive model (RPM), viz.they are charged hard spheres of equal diameter, d = 2, 3, or 4 • 10 −10 m., and charge q α = Z α e, with Z α being the ionic valency.The polyion was ascribed a radius R = 8 • 10 −10 m so that the distance of closest approach of an ion varied from 9 • 10 −10 m to 10 • 10 −10 m.In another case, we have taken R = 20•10 −10 m, which should resemble the distributions corresponding to that in a planar double layer.

Density functional theory
In density functional theory [61], one starts with an expression for the grand potential, Ω, as a functional of the singlet density profiles, ρ α (r), of each of the species, α.At equilibrium the grand potential is minimal with respect to variations in the density profiles, viz., for each α and this condition is used to determine the density profiles and the free energy.Without going into details, we present here the relevant equations for the nonuniform density distribution of small ions in the region r > (R + d/2).In the above expression, ψ(r) represents the mean electrostatic potential and the quantities c (1)hs α (r; [{ρ α }]) and c (1)el α (r; [{ρ α }]) denote, respectively, the hard sphere and electrical contributions to the first order correlation function.The electrostatic potential ψ(r) is obtained from a solution of the corresponding PB equation in cylindrical geometry, and is given as In the absence of explicit expressions, c (1) hs α (r; [{ρ α }]) and c (1) el α (r; [{ρ α }]) have been approximated as where "tilde" (∼) refers to the uniform fluid, with ρ 0 α being the mean ionic number density of the αth species.The quantity ρ(α) (r) is the total weighted density, which has been calculated using the WDA scheme of Denton and Ashcroft [55].The uniform second order correlation functions, c(2 ) are taken from the analytical expressions [62] within the mean spherical approximations for a uniform mixture of charged hard spheres.Thus the density functional theory is completely specified through equations ( 3)-( 6) for calculation of the ionic distribution and the mean electrostatic potential profiles.It ought to be pointed out that the DFT reduces to the nonlinear PB theory [15] once the correlation between the small ions is neglected, i. e., by setting c

Modified poisson Boltzmann theory
As indicated in the Introduction, the MPB equation appropriate for the cylindrical geometry as applied to the present polyion-ion system was formulated by Outhwaite [27].Without delving into the details, which have been described in references [28,29], for completeness we outline here the relevant equations for the ionic density distributions.
In cylindrical symmetry, the MPB polyion-ion singlet density distribution reads where, u(r) = √ rψ(r).In equation ( 7), F and F 0 are given by with and κ and κ 0 are the local and bulk Debye-Hückel screening parameters The exclusion volume term ξ α (r) is approximated using the Bogoliubov-Born-Green-Yvon hierarchy for a planar symmetry [63] where Θ(x) is the Heaviside function and φ(y, x) is the fluctuation potential evaluated on the exclusion surface of the discharged ion.This approximate ξ α (r) was found to be reasonable in the earlier applications [28,29], and as we will see later, is adequate for the physical parameters apropos of DNA utilized in the present work.It should also be noted at this point that MPB theory reduces to the classical nonlinear PB theory upon taking ξ α (r) = 1, F = F 0 and d → 0 [28].We remark however, that the PB theory used here is in the Stern [59] spirit, that is, the distance of closest approach of a small ion to the polyion continues to be R + d/2 so that ξ α (r) = Θ(r − (R + d/2)).

Results and discussion
The DFT and MPB equations have been solved using numerical iterative methods, the details of which can be found elsewhere in the literature (see for example, references [44,47] for DFT and references [64][65][66] for MPB).We have also obtained numerical solution of the PB theory for the cylindrical polyion-ion system in question for comparison purposes.Some of the physical parameters utilized in the calculations have been mentioned earlier in sub-section 2.1.We have treated 1:1, 2:1/1:2, and 2:2 valency electrolyte systems in the concentration range 0.01 mol/dm 3 -2 mol/dm 3 , while the surface charge density on the polyion varied from 0 to 0.3 C/m 2 .At the ionic diameter d = 4 • 10 −10 m some calculations have also been performed at a higher polyion radius R = 20 • 10 −10 m.
We first discuss the zeta potential profiles for a cylindrical double layer system as a function of polyion surface charge density σ, bulk salt concentration c, ionic valency Z α , and ionic size (diameter) d.The zeta potential is defined as the mean   electrostatic potential ψ(r) at the closest approach between a small ion and the polyion and is given as The DFT, MPB, and PB zeta potentials are presented as a function of σ in figures 1-3 for 1:1, 2:1/1:2, and 2:2 valency electrolytes, respectively.The sub-figures a-c in each figure correspond to d = 2, 3, and 4 • 10 −10 m, respectively, while in each sub-figure the results for different c are shown.Some general trends of the results are clear from the figures.Firstly, at low concentrations, the dependence of ζ on ionic diameter is rather weak.For instance, the curves at c = 0.01 mol/dm 3 corresponding to different d are quite similar, both qualitatively and quantitatively.Secondly, the magnitude of ζ decreases with salt concentration for all ionic valencies and all ionic diameters.Thirdly, the classical mean field ζ generally tends to overestimate the corresponding statistical mechanical results.Fourthly, the deviation between the PB ζ and the DFT, MPB ζ increases as the ionic size decreases.These trends are consistent with similar observations for the PDL [50,51,60].The relative magnitude of the deviation mentioned above is however, the least for 1:1 systems (cf.figure 1) and increasing for asymmetric and higher valency systems (cf.figures 2 and 3).This is of course a well known feature of the mean field result that the theory is more useful for monovalent systems where the effect of the neglected inter-ionic correlations is relatively smaller than in the presence of multivalent ions.In figure 1c (1:1 salt, d = 4 • 10 −10 m) the zeta potential corresponding to the bigger polyion radius R = 20 • 10 −10 m (curves b) is seen to lie above that for R = 8 • 10 −10 m (curves a) for all the theories.This is expected since for higher R the system approaches the PDL.We note further that the MPB results for d = 4 • 10 −10 m are reminiscent of results seen in the earlier MPB application [28,29] with very similar physical parameters, viz., d = 4.25 • 10 −10 m.One of the more noteworthy features of figure 1 is the consistency of the DFT and MPB predictions for all ionic sizes and all polyion radii, the agreement between the two theories increasing as d increases from 2 • 10 −10 m through 3 • 10 −10 m to 4 • 10 −10 m.At the highest d treated here, the DFT and MPB curves are virtually indistinguishable at c = 1 mol/dm 3 .A similar level of agreement between the two approaches has been reported for the PDL [60].
Turning now to 2:1/1:2 electrolytes in figure 2, we notice immediately that the behaviour of the zeta potential for 2:1 valency systems (monovalent counterions) is very similar to that in 1:1 valency systems.This is a well established result in the double layer literature, which testifies to the fact that it is the electrode-counterion interaction that shapes the system characteristics.In particular, the DFT and MPB results again reveal consistency although the differences between them for the two smaller ionic diameters at the highest concentration are comparatively larger than that for the earlier 1:1 system.The situation, however, changes for divalent counterions (1:2 valency systems) when the zeta potential is seen to pass through a minimum.This is true at almost all the concentrations and at all the ionic diameters studied here.Although the theoretical predictions now reveal more deviations between them, the overall agreement between the two theories again improves as the ionic diameter increases.In the absence of exact simulation data, it is difficult, however, to make any definitive statement about the relative accuracy of either the DFT or the MPB theory.The PB theory does not capture the above effects and is not even in qualitative correspondence with the DFT and MPB.As we move on to 2:2 electrolytes as shown in figure 3, the DFT and MPB zeta potential shows a maximum and then starts decreasing.Expectedly, the 2:2 results now show similarities with the 1:2 results and similar considerations as for the earlier seen latter results apply.Such maxima or minima have been previously predicted theoretically through simulations for the PDL [6,50,51,60] and the SDL [53,[67][68][69], and have also been observed earlier for the CDL [19,22,28,29].The maximum or minimum in the ζ for an electrolyte with a divalent counterion implies that the differential capacitance can become infinite before becoming negative [50].In the next two figures (figures 4 and 5) we compare the theoretical results with some published MC data [42] for ionic density distributions (in the form of local concentration profiles conc α (r)[= cρ α (r)/ρ 0 α ]) for 1:1 and 2:1 salts at c = 0.016 mol/dm 3 and σ = 0.187 C/m 2 .Figure 4 shows the results for a 1:1 salt where the reduced mean electrostatic potential ψ * (r) [= βeψ(r)] is also given in the top panel.The DFT and MPB local concentration profiles reproduce the MC data quite accurately.Note also that these profiles show little structure.Not surprisingly, for such a low concentration 1:1 system the PB results are also in good agreement with the simulations.The predicted ψ * (r)s also show similar consistency.A slightly more structure begins to appear in the polyion-ion distributions for a 2:1 valency system shown in figure 5. Again the DFT and MPB results generally agree well with the simulation data, with the MPB co-ion profile tending to follow the MC data more closely.The PB results now show a greater deviation than in the 1:1 case and predict a more diffuse, thicker double layer.To investigate the effect of ionic size on the double layer structure, we present, in figures 6 and 7, the local concentration and potential profiles for 1:1 and 2:1 electrolytes, respectively, at different ionic diameters.These calculations are at c = 1 mol/dm 3 , and in each figure the panels a, b, and c correspond to d = 4, 3, and 2 • 10 −10 m , respectively.A glance at the figures across the panels from left to right reveals the principal effect of ionic size on structure -the double layer becomes thicker and diffuse as the ionic diameter decreases.For example, in figure 6a, the concentration profiles reach their bulk values at around r ∼ 5d from the polyion axis.In figures 6b, 6c this occurs at r ∼ 7d and r > 8d, respectively.For 2:1 salts, these numbers are r ∼ 5d for d = 4 • 10 −10 m (figure 7a), r ∼ 6.5d for d = 3 • 10 −10 m (figure 7b), and r ∼ 8d for d = 2 • 10 −10 m (figure 7c).This feature is corroborated by the corresponding potential profiles shown.Decreasing the size of the ions leads to smaller accumulation near the polyion, which in turn yields a more diffuse double layer.The DFT and MPB ψ * (r) for the 2:1 salt systems at d = 4, and 3 • 10 −10 m shown in figures 7a and 7b dip below zero near the polyion surface.This is indicative of polyion overcharging, that is, more counterions accumulate in the vicinity of the polyion than is necessary to screen the surface charge.This has been observed in earlier double layer studies involving, besides the DFT and MPB, other theoretical and numerical methods [22,26,24].We refer the reader to a recent excellent review on macroion overcharging by Quesada-Perez et al. [70] and references therein.At a given electrolyte concentration overscreening is also a function of ionic size as it is not observed at d = 2 • 10 −10 m and c = 1 mol/dm 3 .This is not to say that overscreening or charge inversion does not occur for a 2:1 salt at this value of d or for that matter for a 1:1 salt.Indeed there is evidence in the literature to suggest that overscreening in these systems is possible at sufficiently high concentrations [28,34,46].The DFT and MPB results again show remarkable consistency for both 1:1 and 2:1 salts at all ionic diameters.As expected, the PB profiles are closer to their statistical mechanical counterparts in 1:1 systems but deviate more in 2:1 systems.Further, in the latter case the PB ψ * (r) are monotonic in variance to the DFT and MPB ψ * (r) and do not show any overcharging effect.Evidently the phenomenon has origins in the inter-ionic correlations and cannot be accounted for within the traditional PB framework unless counterion binding is invoked [6,70].Increased electrostatic attraction between the polyion and a divalent counterion can play a major role in characterizing the polyion-ion distributions and mean elec-trostatic potential profiles.This is clearly noticeable in figures 8 and 9, which show the structure of a double layer containing 1:2 and 2:2 electrolytes, respectively, at c = 1 mol/dm 3 and d = 4 • 10 −10 m.The conspicuous structural features here are the pronounced oscillations in the polyion-ion distributions and the appearance of a substantial well in the mean potential around the distance where the local concentrations of the ions cross each other.These would imply a rather large charge inversion.This is examined in some details in figures 10 and 11 where the 1:2 and 2:2 ψ * (r) are shown for the three ionic sizes.Once again we notice that a decreasing ionic size leads to a thickening of the double layer.Parallely, the magnitude of the potential minimum also decreases.It is relevant to note that the comparative behaviour of the DFT and the MPB theory seen earlier [60] extends to these situations.Although the predicted results are generally consistent over the range of the ionic size studied in this work, the agreement diminishes a little at the lower diameters.The PB theory on the other hand deviates substantially from the formal theories in all cases with the deviation increasing with decreasing ionic size.

Concluding remarks
The density functional and the modified Poisson-Boltzmann theories have been employed in an in-depth study of the equilibrium properties of a primitive model cylindrical double layer with a particular focus on the effect of ionic size on double layer structure.The principal findings are (i) the double layer tends to become more diffuse, thicker as the ion-size decreases, and (ii) the classical Poisson-Boltzmann description of structure becomes progressively poorer with decreasing ion-size.Indeed the PB results are not even qualitatively similar to the statistical mechanical results in some instances, for example, at higher concentrations and/or when asymmetric or higher valency electrolytes are present.These results re-affirm similar conclusions reached earlier for the corresponding planar double layers.[50,51] Other relevant results include (i) the appearance of a maximum or a minimum in zeta potentials in a double layer containing electrolytes with multivalent counterions, and (ii) the phenomenon of polyion overcharging in certain cases.The former clearly indicates a change in sign in the differential capacitance in presence of multivalent counterions.These features are tied to the inter-ionic correlations since such features cannot be captured by the PB theory.
The relative behaviour of the DFT and MPB approaches is interesting.Overall the theoretical predictions are reasonably consistent over the range of physical parameters studied.Similar behaviour was observed for these theories in their application to the PDL [50,51].Although in contrast to the planar situation the limited availability of simulation data at the present parameters makes a more detailed assessment of the theories difficult, the general consistency of the polyion -small ion distribution functions and the mean electrostatic potentials out of the two very different theories adds to the reliability of the results.Another detailed simulation study of the CDL would be quite useful nonetheless, and such a study is planned for the future.

Figure 1 .
Figure 1.Zeta potentials for a 1:1 electrolyte as a function of surface charge density of the polyion at different electrolyte concentrations (in Molar (mol/dm 3 )) with the ionic diameter (a) d = 2 • 10 −10 m, (b) d = 3 • 10 −10 m, and (c) d = 4 • 10 −10 m.The dot-dashed, solid, and dashed curves are predictions of the DFT, MPB, and PB theories, respectively.The legend a and b in the bottom sub-figure corresponds to two values of polyion radius, R = 8•10 −10 m and R = 20•10 −10 m, respectively.

Figure 4 .
Figure 4. Concentration profiles (lower panel) of counterions and coions around a polyion of radius R = 8•10 −10 m and surface charge density σ = 0.187 C/m 2 for a 0.016 mol/dm 3 1:1 electrolyte.The plots in the upper panel are the corresponding reduced mean electrostatic potentials.The symbols are the MC simulations, and the dot-dashed, solid, and dashed curves are predictions of the DFT, MPB, and PB theories, respectively.MC data from reference [42].

Figure 5 .
Figure 5. Concentration profiles (lower panel) of counterions and coions around a polyion of radius R = 8•10 −10 m and surface charge density σ = 0.187 C/m 2 for a 0.016 mol/dm 3 2:1 electrolyte.The plots in the upper panel are the corresponding reduced mean electrostatic potentials.Notation as for figure 4. MC data from reference [42].

Figure 6 .
Figure 6.Concentration profiles (lower panels) and the reduced mean electrostatic potential profiles (upper panels) of a 1 mol/dm 3 1:1 electrolyte around a polyion of radius R = 8 • 10 −10 m and surface charge density σ = 0.187 C/m 2 .The dot-dashed, solid, and dashed curves are predictions of the DFT, MPB, and PB theories, respectively.The legend a, b, and c corresponds to the three values of the ionic diameters, d = 4, 3, and 2 • 10 −10 m, respectively.

Figure 8 .
Figure 8. Concentration profiles (lower panel) and the reduced mean electrostatic potential profiles (upper panel) of 1 mol/dm 3 1:2 electrolyte around a polyion of radius R = 8 • 10 −10 m and surface charge density σ = 0.187 C/m 2 with the value of the ionic diameter d = 4 • 10 −10 m.Notation as for figure 6.

Figure 9 .
Figure 9. Concentration profiles (lower panel) and the mean electrostatic potential profiles (upper panel) of 1 mol/dm 3 2:2 electrolyte around a polyion of radius R = 8 • 10 −10 m and surface charge density σ = 0.187 C/m 2 with the value of the ionic diameter d = 4 • 10 −10 m.Notation as for figure 6.