Charge neutrality breakdown in confined aqueous electrolytes: theory and simulation
We study, using Density Functional theory and Monte Carlo simulations, aqueous electrolyte solutions between charged infinite planar surfaces, in a contact with a bulk salt reservoir. In agreement with recent experimental observations [Z. Luo et al., Nat. Comm. 6, 6358 (2015)], we find that the confined electrolyte lacks local charge neutrality. We show that a Density Functional Theory (DFT) based on a bulk-HNC expansion properly accounts for strong electrostatic correlations and allows us to accurately calculate the ionic density profiles between the charged surfaces, even for electrolytes containing trivalent counterions. The DFT allows us to explore the degree of local charge neutrality violation, as a function of plate separation and bulk electrolyte concentration, and to accurately calculate the interaction force between the charged surfaces.
In physics, as well as in chemistry, one often finds situations in which electrolyte is confined between charged surfaces. These surfaces may belong to electrodes, macromolecules, charged colloidal particles, polymers, etc., and can give rise to complex ionic distributions generally known as Electrical Double Layers (EDLs). Presence of electrolyte between charged surfaces strongly affects their interaction and can lead to fascinating phenomena such as like-charge attraction Levin (2002); Linse and Lobaskin (1999); Ŝamaj and Trizac (2011) and charge reversal Pienegonda, Barbosa, and Levin (2005); Fernandez-Nieves et al. (2005). The interaction between EDLs is fundamental for understanding colloidal stability and efficient energy storageLevin (2002); Likos (2001). Electrolytes confined by porous walls show promising application as supercapacitors, since carbon-based electrodes are known to increase the storage performance of these devices due to their high specific surface area Zhang and Zhao (2009). The porosity brings about strong ionic confinement within the nano-sized pores. Despite intense exploration, many questions are yet to be elucidated regarding the behavior of strongly correlated Coulomb systems in a confined environment Zhang and Zhao (2009).
From a theoretical perspective, the earliest attempts to describe the properties of electrolytes in a close vicinity of charged surfaces go back to the pioneering work of Gouy, Chapman, and Stern (GCS) Gouy (1910); Chapman (1913); Stern (1924). In the original approach, the double layer structure was described as a thin layer of counterions, condensed onto a charged surface – the so-called Stern-layer – followed by a diffuse region in which ionic distributions rapidly decayed to their bulk values Zhang and Zhao (2009). This simple picture was based on the mean-field Poisson-Boltzmann (PB) theory. It provided a clear physical explanation for a number of interesting phenomena inherent to electrical double layers – ranging from the charge renormalization of charged surfaces Levin (2002); Trizac, Bocquet, and Aubouy (2002) to the capacitance of simple electrodes Fedorov and Kornyshev (2014); Zhang and Zhao (2009). For this reason the GCS theory has remained very popular in both electrochemistry and biophysics literature Fedorov and Kornyshev (2014). Care, however, is required when extrapolating this simple picture to the case of strongly correlated electrolytes – such as aqueous electrolytes made of multivalent ions, organic electrolytes, or ionic liquids. In these cases, the structure of the EDL is dominated by ionic correlations Paillusson and Trizac (2011), giving rise to complex phenomena such as charge reversal in multivalent ionic electrolytes and layering in ionic liquids Wu et al. (2011). Even in situations of low electrostatic couplings, the packing effects in strongly confined electrolytes can alone be responsible for oscillatory profiles which cannot be captured at a mean-field level of GCS theory Frydel and Levin (2012).
Over the last decades, a number of different approaches have been put forward to quantitatively describe the ionic correlations in EDLs beyond the mean-field theory Attard (1996). A great progress has been achieved with computer simulations, in part due to rapidly increasing computational power, but also because of the development of new techniques which allow to efficiently handle the long-range Coulomb interactions in slab geometry Yeh and Berkowitz (1999); Arnold, de Joannis, and Holm (2002); dos Santos, Girotto, and Levin (2016). The theoretical advances have relied on extending the integral equations (IE) methods developed for bulk systems Henderson and Blum (1978); Henderson, Blum, and Smith (1979); Vericat, Blum, and Henderson (1982); Ballone, Pastore, and Tosi (1986) to inhomogeneous ones Kjellander and Marčelja (1986); Kjellander (1988a, b); Kjellander et al. (1992); Greberg and Kjellander (1994), as well as to interacting EDLs Lozada-Cassou (1984); Lozada-Cassou and Henderson (1986); Lozada-Cassou and Díaz-Herrera (1990); Alejandre et al. (1990). Other approaches explored Modified Poisson-Boltzmann (MPB) equation Outhwaite, Lamperski, and Bhuiyan (2011); Bhuiyan et al. (2007); Lamperski, Outhwaite, and Bhuiyan (2009), and a classical Density Functional Theory (DFT) Jaing et al. (2014). The main advantage of using DFT is that various contributions to free energy can be handled separately, allowing distinct approximations in their calculations. For instance, it is well known that the excluded volume interactions can be described to a very high degree of accuracy using the Fundamental Measure Theory (FMT) Rosenfeld (1989, 1990, 2002); Roth (2010), and can be conveniently decoupled from the electrostatic interactions. The DFT can be constructed using relatively simple arguments based on truncation of functional expansions, introduction of coupling parameters, local or weighted density approximations, or even combinations of these Jaing et al. (2014); Yang and Liu (2015). All these approximations are easy to control and possess quite transparent physical interpretations, allowing for their systematic improvement. This should be contrasted with IE techniques, for which the approximations are based on closure relations which appear in diagrammatic cluster expansions for the correlation and bridge functions Attard (1996); Lozada-Cassou (1992), the degree of accuracy of which is, in general, not known a priori.
It is usual to take for granted the local charge neutrality between charged surfaces dos Santos and Levin (2015); Bakhshandeh et al. (2015); Zhang, Davis, and White (1993) when dealing with EDL in both theoretical and computational approaches. This is not to be confused with the overall charge neutrality, which clearly has to be always satisfied for a charged system as whole. In order to avoid any confusion, we emphasize that from now on the term charge neutrality will refer to the local electroneutrality of electrolyte confined between charged walls. The possibility that electrolytes in contact with a bulk reservoir might lack charge neutrality in situations of strong confinement has been already discussed both theoretically Lozada-Cassou, Olivares, and Sulbarán (1996); Lozada-Cassou et al. (1996) and observed in computer simulations Lo et al. (1998); Lee and Chan (1997); Lee et al. (1999). However, these findings have only very recently been confirmed experimentally in the work of Luo et al. Luo et al. (2015), who demonstrated that electrolytes confined in narrow pores do not in general exhibit electroneutrality. To explore this issue Luo et al. used Nuclear Magnetic Resonance (NMR) techniques to demonstrate that electrolytes confined by nanoporous carbon with graphite-like internal surfaces can violate charge neutrality. They have observed a dependence of charge neutrality breakdown on ion specificity which followed the Hofmeister series. This is consistent with the recent theoretical works which show that ionic specificity arises from a combination of hydrophobic, dispersion, and polarization interactions, which are very different for , , , and near hydrophobic surfaces Levin, dos Santos, and Diehl (2009); dos Santos and Levin (2011); Levin and dos Santos (2014); Jungwirth and Tobias (2006, 2002); Ghosal et al. (2005); Jungwirth and Tobias (2001); Kunz, Nostro, and Ninham (2004). Although this ionic specificity is extremely interesting, we will not address it in this paper. Instead we will explore the role that electrostatic correlations Levin (2002) play in charge neutrality breakdown. A similar line of investigation has been taken by Lozada-Cassou and co-workers Lozada-Cassou, Olivares, and Sulbarán (1996); Lozada-Cassou et al. (1996), who used a three-point extended version of the traditional IE formalism to address the problem of electroneutrality violation in electrolytes confined by parallel charged surfaces. In this approach, the two-plate system was modeled as “dumbbell-like” molecules at infinite dilution, and the resulting ionic profiles were calculated as thee-body wall-wall-ion correlations Lozada-Cassou (1984); Lozada-Cassou and Henderson (1986); Lozada-Cassou and Díaz-Herrera (1990); Lozada-Cassou and DíazâHerrera (1990); Alejandre et al. (1990); Feller and McQuarrie (1994); Jiménez-Ángeles, Odriozola, and Lozada-Cassou (2009). The ionic correlations were treated altogether at the level of the Mean Spherical Approximation (MSA). It is, however, well known that in situations in which electrolyte is bounded by narrow pores – where violation of electroneutrality is expected to take place – size correlations between the confined ions may strongly influence the structure of the EDLs Zhou and Stell (1989); Henderson, Sokolowski, and Wasan (1998); Smith and Vörtler (2003); Kamalvand, Keshavarz, and Mansoori (2008). Furthermore, electrostatic correlations in electrolytes containing multivalent ions usually require use of more accurate approximations. It is, therefore, unclear to what extend can strong correlations be captured at the MSA level. In order to further explore the question of how important are the correlation effects in determining electroneutrality violation we shall, therefore, apply a somewhat different approach. The size correlations will be described through the FMT approach, whereas the electrostatic correlations will be taken into account at the level of the hyppernetted chain (HNC) approximation. Furthermore, a recently developed MC simulation technique, well suited to simulate ionic systems in slab geometries at very low computational cost will be applied dos Santos, Girotto, and Levin (2016) to this system. We will show that for large separations between surfaces the macroscopic charge neutrality is restored and the traditional Donnan approach, in which a potential difference across the system-reservoir interface is established to force the overall electroneutrality (Donnan equilibrium) Donnan (1924); Ohshima and Ohki (1985), can be applied. On the other hand, for narrow pores the interior of the pore is not electroneutral, in which case the internal charge of the pore will be determined by the chemical equilibrium with a bulk reservoir. These findings are in qualitative agreement with earlier theoretical predictions obtained in Refs. [Lozada-Cassou, Olivares, and Sulbarán, 1996] and [Lozada-Cassou et al., 1996]. We will also investigate to what extent the absence of electroneutrality between the confining plates influences the net forces between the charged surfaces. It is important to emphasize that such effects from local violation of electroneutrality on the interaction between double layers can have important implications in a number of biological and physical systems composed by charged objects immersed in an electrolyte solution Levin (2002). Examples range from colloidal self-assembly Odriozola and Lozada-Cassou (2013, 2016); Furst (2015); Cademartiri and Bishop (2015) to membranes Jiménez-Ángeles and Lozada-Cassou (2004) and biological cells Gelbart et al. (2000), as well as confined phase separations Lozada-Cassou and Yu (1996), ionic correlations across charged membranes Lozada-Cassou and Yu (1997) and colloidal systems confined by charged walls Jiménez-Ángeles, Odriozola, and Lozada-Cassou (2009). For a recent overview on the way EDL interaction might influence the assembly of soft-matter systems, we refer the reader to Ref. [Editorial, 2015] and references therein. In spite of these numerous applications, the role that the charge neutrality violation plays in the EDL interactions has been mostly neglected after the original contributions of Lozada-Cassou et al. Lozada-Cassou (1984) which addressed this issue. The aim of the present work is to investigate the effects that addition of salt and electrostatic correlations have on the charge neutrality violation, and explore how it influences the interaction between planar EDLs.
The work is organized as follows. In Sec. II, a description of the system under consideration is briefly given. Some details about the MC techniques and the DFT are presented in Secs. III and IV, respectively. Results are shown in Sec. V and the Conclusions are outlined in Sec. VI. Technical details regarding some of the specific calculations employed in this work can be found in the Appendices.
Ii Model System
Throughout the paper we will use the framework of the Primitive Model (PM) of electrolytes, in which the ions are represented as hard spheres with point charge embedded at the center. The solvent is modeled as a uniform dielectric continuum with its molecular structure disregarded. We will consider a system composed of two positively charged walls located along the direction at positions and . Both walls have a surface charge density uniformly distributed over their infinite surfaces. The space between the walls is occupied by an aqueous electrolyte (permittivity with the dielectric constant of vacuum) composed of cations of valence and anions of valence . The confined electrolyte is in contact with a bulk reservoir of concentration . In simulations, we will consider that the whole system (electrolyte plus reservoir) is confined by neutral walls symmetrically located at positions , such that (see Fig. 1). Due to the symmetry in the plane, all the system inhomogeneity takes place along the direction. Although the net charge in the space between the walls is not explicitly taken to be zero, it is important to point out that the system as a whole does obey the global charge neutrality, which is expressed as:
where represents the ionic profile of component at position , and is its corresponding charge (with being the charge of a proton).
Iii Model and Monte Carlo Simulations
The difficulty of simulating systems with long-range interactions is that the interaction potential can not be cutoff and one can not use simple periodic boundary conditions. Instead one must consider an infinite set of replicas of the simulation cell, so that the ions in the principal cell feel the electric field produced by the ions from all the replicas. To efficiently sum over the replicas one must use some form of Ewald summation D. Frenkel and B. Smit (2002). There are a number of very efficient implementations of Ewald summation for isotropic 3d Coulomb systems. The situation, however, is more complicated for systems with reduced symmetry, such as our two infinite charged plates, in which case the system is periodic only in two out of three dimensions. Recently we have introduced a method which allows for a very efficient adaptation of 3d Ewald summation to systems with reduced slab geometry, with confining charged surfaces. The fundamental idea of the method is to separate the electric field produced by the infinite charged plates from the rest of the system. The difficulty, however, is that when this is done, the residual systems is no longer charge neutral, so that the electrostatic energy of the infinitely replicated system will diverge. We showed, however, that this divergence can be renormalized away, allowing us to calculate finite renormalized electrostatic energy. The details of all the derivations can be found in Ref. [dos Santos, Girotto, and Levin, 2016]. Here we will just present the applications of the algorithm to the present problem.
Our simulation cell has the volume , with Å and Å. The electrolyte is in the region (with ), and the pair of charged plates are positioned at and . In the regions and there is pure water, see Fig. 1. Water is treated as a continuum of dielectric constant . The Bjerrum length, defined as , is Å. This value is appropriate for water at room temperature. All ions have radius Å. The plates are equally charged with a positive surface charge density in the range . The effect of counterion valence is explored by changing between 1:1, 2:1, and 3:1 electrolyte.
To efficiently sum over all the replicas the electrostatic potential is split into long and short-range contributions, in addition to the potential produced by the charged surfaces, see Ref [dos Santos, Girotto, and Levin, 2016], and a residual potential of self-interaction. The short-range electrostatic potential at position is
The dumping parameter is set to . The long-range electrostatic potential is
The number of -vectors defined as , where are integers, is set to around in order to achieve fast convergence. The electrostatic potential produced by the charged surfaces is
The self-interaction potential has the form
We can now easily compute the total ionic electrostatic energy
To perform MC simulations we use Metropolis algorithm with MC steps to achieve equilibrium. The profile and force averages are performed with uncorrelated samples. During the equilibration, we adjusted the length of the particle displacement to achieve an acceptance of trial moves near . We are particularly interested in the net charge of electrolyte in between the confining surfaces. The external electrolyte acts as a reservoir for the internal region, so that the ions are allowed to freely move across the charged surfaces. The interaction between two surfaces is modulated by the external and internal electrolyte, and has both electrostatic and entropic contributions. To calculate the mean electrostatic force we use the method of virtual displacement in which one of the plates is moved while the other plate and all the ions remain fixed, which implies that the electrostatic force per unit area in the direction is
where is the area of the plate, , and
The first term on the right-hand side of Eq. (7) is the mutual force between the charged wall surfaces, while the second term represents the ionic-averaged electrostatic forces on the first wall. Note that a positive force means repulsion between the walls. Eq. (7) can be simplified to yield,
where and represent the total charge per unit of area located in the regions and , respectively.
To calculate the entropic force, which arises from the transfer of momentum in the collisions between ions and plates, we use the method introduced by Wu et al. Wu et al. (1999) It consists of performing a small virtual displacement of the plates along the direction – while all the ions remain fixed – and counting the number of resulting virtual overlaps between the plates and the ions. The entropic force per unit area can then be written as
where is the number of virtual overlaps between the plates with the ions after a small displacement Å that brings plates closer together (superscript stands for closer) and is the number of overlaps of the plates with the ions after a displacement that moves the two plates farther apart (superscript stands for farther) Colla, dos Santos, and Levin (2012).
The profiles were made counting the average number of particles in a volume range running over the direction. The value of was set to Å.
Iv Density Functional Theory
The DFT assumes that the functional of a set of densities ,
where and are respectively the chemical potential and the external potential acting on the ion of type , achieves a minimum at equilibrium Evans (1992). This minimum corresponds to the system’s grand potential (or the ground-state energy of a quantum system Mermin (1965)). The functional is the intrinsic free energy, since it contains all the information about the particle interactions, regardless of any particular external potential acting upon the system. It can be split into an ideal gas contribution plus an excess free energy , in which all contributions resulting from the particle interactions are included Evans (1979); Singh (1997); Löwen (2002); Evans (1992); Wu and Li (2007). The ideal gas contribution is:
where is the usual de Broglie wavelength. Use of this exact relation together with a straightforward application of the Euler-Lagrange minimization condition to Eq. (11) provides the following equilibrium distributions:
where , and defines the first-order direct correlation function. If the density profiles corresponding to a given external potential are known, the equilibrium distributions (13) can be easily written in terms of such reference profiles as:
where are the reference density profiles, with and . Unfortunately, the intrinsic energy, , is not known for arbitrary density distributions. However, it can be very accurately calculated, for example, for a bulk system which can then be used as the reference state. In general, if a suitable approximation for the reference fluid with densities is known, Eq. (14) can then be applied in a perturbative scheme to obtain improved corrections for the desired distributions , as we will see later on.
So far no approximations have been made in writing Eqs. (13) and (14). The excess free energy can be further split in accordance with different particle interactions. In the present case of ionic systems in the framework of the PM approach, the functional is a combination of hard-core and electrostatic contributions, . Accordingly, the single-particle direct correlation function can be written as a sum of hard-core and electrostatic contributions, . This possibility of treating separately the different contributions for the particle correlations is one of the major advantages of using a DFT approach. The hard-core contribution can, for example, be treated in the framework of the Fundamental Measure Theory (FMT) Rosenfeld (1989, 1990, 2002), which has proven to be extremely accurate in describing both thermodynamic and structural properties of hard- sphere system for a variety of packing fractions and size asymmetries Roth (2010); Davidchack et al. (2016). We will, therefore, apply the FMT for calculating the ionic exclusion volume interactions, using its more recent White-Bear formulation Roth et al. (2002); Yu and Wu (2002), the details of which are outlined in Appendix A. For a deeper analysis of the hard-sphere FMT we refer the reader to Roth’s recent review on this topic, Ref. [Roth, 2010].
The difference between the electrostatic single-particle direct correlation function and the reference state, , can be formally evaluated by introducting a parameter to continuously interpolate between the equilibrium and reference states. The simplest choice is to use a linear path, such that , with and . When varies from zero to unity, the profiles smoothly change between the reference state and the desired equilibrium distributions. The corresponding variation in the single-particle electrostatic correlations with respect to their reference counterpart can be written as:
where represent the direct correlation functions evaluated when the profiles are . The derivative on the right-hand side of the last equality can be obtained by noting that the single-particle correlations depend on only implicitly, through the density profiles . It can, therefore, be related to the changes in the density distributions as J. -P. Hansen and I. R. McDonald (2006):
where in the second equality the definition of the direct pair correlation function has been used. With the particular choice of the linear dependence of on , the derivative on the right-hand side of Eq. (16) reduces to , and the changes in the direct correlation functions described in Eq. (15) become:
In the above relation, we have defined the mean pair direct correlation function as:
Even though the relations (17) and (18) are formally exact, they require knowledge of of the direct pair correlation functions for inhomogeneous systems with densities between reference and the equilibrium state, which in general are not known. Further progress can be achieved by removing the long-distance Coulomb interaction from the pair direct correlation functions in Eq. (18), which are expected to behave as at large particle separations. The resulting short ranged direct correlation can be used together with Eq. (18) to rewrite Eq. (17) as:
where is the electrostatic potential arising when the ionic density profiles in the reference state are perturbed by an amount . In the second term on the right-hand side of Eq.(19), we have defined the residual direct correlation function as,
where the mean short range correlation is defined in Eq. (18), with the integral in performed over . Clearly, the first term on the right-hand side of (19) corresponds to the mean-field electrostatic potential, while the residual electrostatic single-particle correlations contain all the correlation effects beyond the mean-field. Combining the above results with the Euler-Lagrange equation, Eq. (14), the ionic profiles corresponding to an arbitrary external potential acting on the ionic system can be finally written as:
In the case of planar double layers, where the external potential is provided by the electrostatic and hard core ion-wall interactions, the density profiles change only along the direction perpendicular to the parallel plates, and Eq. (21) reduces to:
where now represents the (dimensionless) total electrostatic potential change with respect to the reference state (including the ionic interactions with the wall). Here, it is assumed that the hard core ion-wall interaction is the same for both equilibrium and the reference state. Eq. (22) represents the basic relation for the planar electric double layer structure in the framework of the DFT approach. In practice, it has to be coupled with the Poisson equation for the ionic distributions. In the case depicted in Fig. 1 of two parallel charged plates located at positions , Poisson equation is:
where is the surface density of charge carriers on the wall. Once the single-particle hard-core and residual electrostatic correlations are known, Eqs. (22) and (23) can be self-consistently solved (e. g. via a Picard iteration method) to provide the ionic distributions around the charged plates. As already mentioned, the exclusion volume contributions can be evaluated with a high degree of accuracy in the framework of the FMT approach. What essentially determines the ability of different DFT approaches to capturing the fine structure of the EDL are, therefore, the approximations employed in the evaluation of both reference state and the residual electrostatic correlations. Several different approaches have been proposed over the last years to accurately calculate these quantities in a number of different contexts Yang and Liu (2015). Since very little is actually known about the pair correlations of inhomogeneous ionic system, it is in practice usual to take a homogeneous system, with the reservoir bulk concentrations , as a representative reference state. In most cases, a further bulk-like approximation is invoked by taking the effective pair correlation functions in Eq. (20) to coincide with the corresponding bulk ones:
where the correlations on the right-hand side are calculated in the charge-neutral bulk fluid. This approximation, usually referred to in the literature as the bulk expansion approximation, can be easily identified as a second-order truncation in the Taylor functional expansion of the excess coulomb free energy around the reference bulk state Jaing et al. (2014). It was first employed by Rosenfeld Rosenfeld (1993) more than twenty years ago, and has been extensively used since then to describe EDL structures. The bulk expansion approximation is also equivalent to the hyppernetted chain (HNC) approximation for the residual wall-ion correlation J. -P. Hansen and I. R. McDonald (2006). It is employed in a combination with the Mean Spherical Approximation (MSA) for the bulk direct pair correlations (bulk-MSA) for ion-ion interactions, for which an analytical solution is known exactly Blum and Henderson (1981). This leading-order bulk-MSA approximation has been applied with a good success in a number of situations involving both planar Patra (2010); Henderson et al. (2011); Wu et al. (2011); Medasani et al. (2014); Zhou, Lamperski, and Zydorczak (2014) as well as spherical Yu, Wu, and Gao (2004); Li and Wu (2004) double layers, provided the electrostatic coupling inside the EDL is not to high Jaing et al. (2014); Yang and Liu (2015).
Going beyond the bulk-MSA approach, Gillespie et al. Gillespie, Nonner, and Eisenberg (2002, 2003); Gillespie, Valiskó, and Boda (2005) proposed a Reference Density Fluid (RDF) approximation in which the residual electrostatic correlations are expanded around an inhomogeneous reference fluid represented by density profiles that are themselves functionals of the real equilibrium profiles. Having recognized that these profiles do not need to represent real equilibrium states of the system, they have chosen them in such a way as to satisfy electroneutrality at each position, allowing for the use of a local approximation for the evaluation of the pair correlation functions, again in the context of the MSA approach. Another accurate approach for the residual contributions is the so-called Weighted Correlation Approach (WCA), which has been recently developed by Wang et al. Wang, Liu, and Neretnieks (2011a, b); Jiang et al. (2014). The basic idea relies on the interpretation of the effective direct correlations in Eq. (18) as the correlations which are averaged over the path of inhomogeneous densities connecting the reference and the final state. Since the residual correlations are short ranged, the resulting mean correlations in Eq. (20) can be approximated by direct pair correlations weighted over a region whose typical size is the range of the residual correlations Wang, Liu, and Neretnieks (2011a). The question of which approximation – among bulk-MSA, RDF and WCA – is best to describe the EDL structure seems to depend on the particular problem at hand Jaing et al. (2014), although the bulk-MSA is in all cases computationally cheaper than the other ones. Recently, Yang and Lui performed a careful analysis on the accuracy of these approaches within a large range of electrostatic couplings. They conclude that the optimal choice is in general a combination of these methods Yang and Liu (2015), in which the reference profiles used in Eq. (22) and the ones used for calculating in Eq. (20) are decoupled from one another, leading to different approximations for the zeroth and first order residual correlations Yang and Liu (2015). Besides all these expansion-based approaches, in a very recent contribution Roth and Gillespie Roth and Gillespie (2016) proposed a first-principle approach in which the residual single-particle correlations are obtained via a generalization of the MSA bulk excess free energy to an inhomogeneous system, showing promising results with yet relatively low computational cost.
All the aforementioned approximations for the residual contributions involve the application of the MSA in the calculation of the direct pair ion-ion correlations. Since the MSA is based on a linear response approximation for such correlations, it is reasonable to expect that such MSA-based approaches should fail at sufficiently high electrostatic couplings – where the MSA proves to be inaccurate even for bulk solutions Levin and Fisher (1996). In these limits, strong non-linear effects such as ionic association start to take place Levin (2002); Paillusson and Trizac (2011), therefore requiring more sophisticated approaches for the ionic correlations. One easy way to circumvent this problem is to adopt an approximation in which the electrostatic bulk correlations are evaluated in the framework of the HNC relation, whereby the total and direct pair correlations are related by
where is the radial total correlation function. The MSA relation is clearly recovered upon linearization of the exponential factor on the right-hand side. This relation is complemented by the Ornstein-Zernike (OZ) equation for the homogeneous bulk electrolyte, which in the Fourier space reads as:
These relations can be numerically evaluated to obtain the ionic direct correlation functions for the bulk system. Although remarkably more accurate than its MSA counterpart, the HNC approach has the disadvantage of not possessing an analytic solution, which makes its direct application in the context of RDF or WCA approaches quite difficult from a computational perspective. However, the implementation of the HNC approximation for the electrostatic correlations together with bulk expansion, Eq. (24), can be readily accomplish with not significant additional increase in computational effort. Since we aim to investigate the EDL at high electrostatic couplings – namely high ionic valencies and small surface charges – we will, in what fallows, employ (unless otherwise specified) the bulk expansion in combination with the HNC solution obtained from Eqs. (25) and (26) for computing the residual contributions in Eq. (20). While for moderate ionic correlations the MSA-based approach should provide results comparable to the HNC-based approximation, we expect HNC to be much more accurate in the limit of very strong electrostatic correlations – where the validity of even bulk MSA approach must to be put in doubt.
Having established the theoretical basis as well as the simulation methods to be employed, we now briefly outline some aspects of their numerical implementation before proceeding to analyze the properties of the model electrolyte system described in Sec. II.
In practice, the regions outside the charged plates — the bulk — are taken to be large enough to guarantee that the ionic profiles relax to their bulk values sufficiently far away from the plates. Using these bulk concentrations as the reference state, the ionic profiles in Eq. (22) read as:
The bulk densities are set by the reservoir salt concentration, and should satisfy the overall electroneutrality condition . The potential represents the exclusion volume interaction between the ions and the hard walls located at positions :
Making use of the one-dimensional Poisson equation, Eq. (23), the difference between inhomogeneous and bulk potentials can be conveniently written as a simple functional of the ionic profiles as:
where is the electrostatic potential produced by the charged walls, as defined in Eq.(8), and represents a position deep inside the bulk solution, . When the single-particle correlations are set to zero in Eq. (27), combination with Eq. (29) recovers the traditional Poisson-Boltzmann equation for the ionic distributions. Since both residual and hard sphere contributions are also explicitly written as functionals of the ionic profiles, Eqs. (27) and (29) can be numerically solved with a standard Picard-like iteration method. Starting with guess functions , the corresponding electrostatic profiles are calculated from Eq. (29), along with the FMT hard sphere contributions , as described in Appendix A. As for the residual single-particle correlations , the bulk fluid expansion resulting from the application of Eqs. (24) and (20) is employed. Here, the direct pair correlations for the bulk system are calculated via the solution of the OZ equation, Eq. (26), together with either MSA or the HNC approximations, Eq. (25). While for the MSA case analytic expressions for the pair correlations are available, the HNC approach requires numerical solutions, which are here obtained using the methods described in Ref. [Ng, 1974]. For the computation of , one has to explicitly remove the short range hard sphere contributions from the HNC solution. This is accomplished by taking for the HNC electrostatic correlations. Once all the relevant contributions are calculated, new estimate density profiles are obtained using Eq. (27). The whole process is then repeated until convergence is achieved. In practice, improved estimates for the ionic profiles are obtained by properly combining the corresponding input and output resulting from consecutive iterations. Finally, we notice that the imposed boundary conditions – namely that the bulk concentrations should be recovered far away from the charged walls, , are naturally satisfied within this iterating scheme, provided they are fulfilled by the initial profile guess.
We are now going to use the DFT theory presented above to study ionic distributions between two like-charged surfaces in a contact with a bulk salt reservoir, focusing on the equilibrium ionic profiles, degree of charge neutrality breakdown between the charged walls, and the net force on each surface.
v.1 Ionic profiles
Ionic distributions around two charged surfaces separated by distance Å are shown in Fig. 2, for several salt concentrations and electrolyte asymmetries. Open symbols represent the MC results, the solid lines are profiles obtained from the bulk-HNC approximation, while dashed lines stand for the bulk-MSA approach. We note that the recently developed efficient MC method dos Santos, Girotto, and Levin (2016), allows to very accurately perform statistics while using small grid sizes over the simulation box. The resulting ionic profiles are able to precisely capture the fine details of the double layer structure in the vicinity of the charged surfaces, providing an excellent basis for testing the validity of different theoretical approaches.
Ionic density profiles are strongly inhomogeneous in the region between the plates, and assume their bulk values shortly beyond the charged walls. Exception is the 1:1 electrolyte at small salt concentration (Fig. 2a), in which case the ionic inhomogeneities extends farther away from the walls, which is clearly a consequence of a larger screening length. We can also see that both bulk-HNC and bulk-MSA show excellent agreement with the simulations in the case of monovalent electrolyte (Fig. 2a and Fig. 2b). Since the curves overlap each other, it is actually difficult to distinguish between the two approximations. At such weak electrostatic couplings, the hard core correlations strongly dominate over the electrostatic ones, and the ionic density profiles are well described by the mean-field electrostatic approximation in addition to hard-core corrections Frydel and Levin (2012). As the electrostatic coupling increases – changing the anionic valencies (2:1 in Fig. 2c and Fig. 2d and 3:1 in Fig. 2e and Fig. 2f) – the deviations between different theoretical approaches begin to appear, the bulk-HNC being always the closest to the simulation results. In particular, at higher salt concentrations and strong electrostatic coupling (Fig. 2d and Fig. 2f) the bulk-MSA results fail to reproduce the MC counterion and coion distributions, while the bulk-HNC approximation shows perfect agreement with the simulations. The situation is slightly different at smaller salt concentrations (see Fig. 2c and Fig. 2e), where still the bulk-MSA deviates from the MC data at strong couplings. Although closer to the MC results, the bulk-HNC approximation in this situation is less accurate in comparison with the high salt case. This loss of accuracy is a consequence of the local character of the bulk expansion. While the HNC approach can correctly account for the stronger magnitude of ionic correlations in the case of divalent and trivalent ions, at small ionic strengths the range of the electrostatic correlations becomes larger and, therefore, more important, rendering the bulk expansion less accurate. It is likely that in this limit the non-local approaches such as the RDF or the weighted densities approximations – in which the inhomogeneities are averaged over a distance that scales with the range of the electrostatic correlations – will become more reliable for describing electrostatic correlations. Unfortunately, the implementation of such methods in a combination with the HNC is very difficult. A detailed comparison of the different DFT approaches against the MC results at several salt concentration and ionic asymmetries goes beyond the scope of this work, and will be the subject of future investigations. For now, we simply emphasize that the bulk-HNC is far more accurate than the bulk-MSA approximation and, therefore, in what follows we will exclusively use the bulk-HNC approach in the calculation of the system properties.
It is not clear from the density profiles in Fig. 2 whether or not the electrolyte confined between the charged plates obeys charge neutrality. Notice that nowhere this condition has been explicitly imposed in our calculations. We now investigate under what circumstances will the electroneutrality breakdown takes place. To this end, we define the total charge of the confined electrolyte per unit area as:
The, electroneutrality between the plates, is satisfied if , see Fig. 1. In Fig. 3, the ratio between the magnitude of such “internal” net surface charge and the bare plate surface charge is plotted as a function of the wall separation for different electrolyte asymmetries and salt concentrations. In all cases, breakdown of charge neutrality is found at small wall separations. As the distance between the walls grows larger, the ratio converges rapidly to unity, meaning that the overall electroneutrality in the region between the charged walls is naturally recovered. The crossover distance at which charge neutrality starts to take place is dependent on the amount of salt in the bulk reservoir. Clearly, at small salt concentrations (full curves) larger wall separations are required to guarantee charge neutrality in the interior region. Again, this can be understood in terms of the larger inhomogeneity region across the charged walls at smaller ionic strengths, as has been observed in the left panels of Fig. 2. As the salt concentration increases, the ionic profiles rapidly converge to the bulk value in the outer region, leading to a decrease in neutrality breakdown between the plates. This effect is also consistent with the Donnan approach in which the electrostatic potential difference between the system and the reservoir becomes smaller as the salt concentration increases Tamashiro, Levin, and Barbosa (1998); Ohshima and Ohki (1985).
There are several physical mechanisms responsible for electroneutrality violation in narrow pores. First, electrostatic effects lead to strong counterion condensation at the charged walls, in an attempt to neutralize their surface charge. On the other hand, entropic effects try to induce a homogeneous particle distribution all over the system, thereby limiting counterion condensation. Electrostatic correlations between anions and coions are responsible for complex ionic association, which also effectively reduce the amount of counterion condensation. Moreover, the bulk counterion-counterion correlations favor additional condensation. Finally, exclusion volume correlations strongly constrain the number of ions which can be present between the walls, significantly affecting the charge neutrality of the interior region. Clearly, this effect is much more important for small wall separations, when only few layers of ions can be accommodated between the surfaces.
Fig. 3 also reveals the effect of ionic correlations on electroneutrality. As the counterion valency changes, different qualitative behaviors are clearly observed. In the case of monovalent ions (Fig. 3a), the overall charge between the surfaces grows monotonically as the wall separation increases. Strong violation of charge neutrality is found in this regime, the net charge between the walls being always positive (). Remarkably, at the smallest salt concentration M and at the narrowest separation studied, the net internal charge density is of the bare wall surface charge. Upon salt addition, the magnitude of electroneutrality breakdown becomes weaker, and the separation between the surfaces at which it occurs shorter. A quite different behavior takes place for higher charge asymmetries, as can be observed in Fig. 3b and Fig. 3c. For both 2:1 and 3:1 electrolyte, the magnitude of the net charge between the plates reaches a maximum before charge neutrality is achieved. Furthermore, the phenomenon of charge reversal – whereby the internal double layer has a net charge with sign opposite to the surface charge – is observed in these cases and is a manifestation of strong electrostatic correlations between the multivalent counterions. In the case of 2:1 electrolyte (Fig. 3b) for very narrow slits the internal charge density is smaller than the surface charge. However, when the separation between the surfaces is increased the internal charge overcompensates the surface charge, resulting in a charge reversal. For larger separation between the plates the charge neutrality is restored. The situation is remarkably different for trivalent counterions (Fig. 3b), for which the internal region is always overcharged for very narrow pores – the internal net charge has a sign opposite to the bare surface charge () at salt concentrations M and M. On the other hand, when M (dashed curve) the internal region is undercharged for intermediate separations – the net EDL charge is slightly positive when . Another interesting feature in the case of trivalent counterions is that the breakdown of charge neutrality is significantly more pronounced for very low salt concentration, M, and narrow slits.
To investigate the effect of increasing the surface charge at constant salt concentrations on charge neutrality violation, the ratio as a function of the wall separation for the different charge asymmetries is plotted in Fig. 4. Once again, a different qualitative behavior is observed in the case of monovalent and multivalent ions (Fig. 4a). While a system with 1:1 electrolyte always remains undercharged, and restores the internal charge neutrality only at asymptotically large separations, both 2:1 and 3:1 electrolytes exhibit strong charge reversal for very narrow pores (Fig. 4b and Fig. 4c). For monovalent electrolytes the degree of electroneutrality breakdown becomes weaker as the surface charge increases. On the other hand for very narrow slits, we see that while a system with 2:1 electrolyte is undercharged, for very low surface charge densities, it actually becomes overcharged for large charge density. For 3:1 electrolyte, the overcharging occurs even at very low surface charge densities, and becomes very pronounced as the surface charge is increased.
v.3 Wall forces
We now investigate the interaction between the two surfaces. According to the Contact Value Theorem (CVT), the net force acting on each wall has two contributions, an electrostatic and an entropic. The electrostatic force felt by one of the surfaces is produced by the interaction with the other surface and with the electrolyte, both external and internal ones. The entropic force arises from direct collisions of ions with the surface and the resulting momentum transfer.
Consider the surface located at . The electric field it produces is , (where denotes the sign of ). Newton’s third law then requires that the net force per unit of area, , is
where and represent the total charge per unit of area located in the regions and , respectively. Since for equally charged surfaces the ionic distributions are even functions, , it is easy to see that these quantities are constrained by the condition . The total force on the plate located at can therefore be rewritten as:
Notice that we have defined in such a way that a positive force corresponds to repulsion between the charged surfaces. While the first term on the right-hand side of this equality represents the net electrostatic force per unit of area, the second term is the pressure resulting from the transfer of momentum from the ions to the surface. We notice here that in the present situation the same result is obtained from a direct thermodynamic calculation of the osmotic pressure across the walls, and that it does not depend on the particular approximations employed in building up the excess functional (see Appendix B).
The net force, given by Eq. (32), shows a very non-trivial influence of the electroneutrality on the interaction between the surfaces. When , the net charge between the walls is positive, leading to a repulsive electrostatic force. This, however, implies a net negative charge on the bulk side of the walls, which leads to a larger accumulation of counterions at the surfaces facing the bulk electrolyte. This imbalance of ions on both sides of the surface provides a net force that pushes the plates towards each other, leading to an attractive entropic contribution. In the opposite case where , the inverse situation takes place: the electrostatic interactions become attractive, while the balance of ions between confined and bulk electrolytes leads to a repulsive force. It is the fine balance between these contributions – described by the two terms at the right-hand side of Eq. 32 – that dictates the net, ion mediated, force between the charged surfaces. Fig. 5 shows the net force per unit of area acting on a charged wall at different electrolyte charge asymmetries and salt concentrations. The parameters are the same as in Fig. 3. Although in this case the degree of electroneutrality for monovalent and multivalent ions has quite different behaviors, the net force between the walls has similar form for monovalent and divalent ions (Fig. 5a and Fig. 5b). In both cases, the force is always repulsive and decays monotonically when M and M. At the highest salt concentration M, the force becomes slightly attractive at wall separations in the case of 2:1 electrolyte. In this region, the net charge within the walls is slightly positive, as can be verified in Fig. 3b which shows that . The attractive force induced by the net momentum transfer by the bulk side ions, in this case, overcomes the weakly repulsive electrostatic force. The fact that such effect takes place at higher salt concentrations is due to stronger screening of the electric field on the bulk side of solution, which forces the ions to be accumulated within a narrow region in the vicinity of the wall surfaces (see Fig. 2). The absence of attractive force in the case of surfaces surrounded by monovalent electrolyte has been proven in the framework of the Poisson- Boltzmann theory Trizac and Raimbault (1999); Trizac (2000). Since for 1:1 electrolyte the electrostatic correlations play only a minor role, it is not surprising that no attraction is found between the charged plates under such conditions. It is important to stress, however, that a steric-driven attractive force as well as charge reversal have been already reported in the case of size asymmetric, strongly confined, monovalent electrolytes Guerrero-García, González-Mozuelos, and de la Cruz (2011); Guerrero-García, E.González-Tovar, and de la Cruz (2011). In the case of trivalent counterions (Fig. 5c), the force changes sign at all concentrations under consideration. Again, the low salt situation M differs qualitatively from the other ones. In this case, the force achieves a well defined minimum at , and is attractive even at large wall separations. This result is consistent with the large ratio observed in Fig. 3c, which implies a net negative charge and an electrostatic attraction that strongly overcomes the repulsive pressure provided by the confined electrolyte.
We have investigated the electrostatic properties of an electrolyte confined between charged surfaces in contact with a bulk salt reservoir. A DFT approach combined with a bulk-HNC expansion was employed to calculate the density profiles and the forces acting on the surfaces. Special attention was paid to the local charge neutrality violation in a confined electrolyte. Contrary to the traditional Donnan approach – in which electroneutrality is enforced by the introduction of a potential difference across the system boundaries – in our calculations charge neutrality has not been assumed a priori. The model system can be used to describe the situation in which an electrolyte is confined in carbon nanoporous, for which experimental evidences of local electroneutrality violation has been recently reported Luo et al. (2015). The breakdown of electroneutrality occurs naturally when confined electrolyte is able to exchange particles with a bulk reservoir. Furthermore, the net charge within the confined region can be controlled by electrolyte properties other than ionic specificity, such as the salt concentration and charge asymmetry, as well as the surface charge of the confining walls. In particular, we find that the degree of charge neutrality violation is much more pronounced in the limit of small ionic strengths. These results are in line with the ones obtained previously using a three-point extended HNC-MSA approach Lozada-Cassou, Olivares, and Sulbarán (1996); Lozada-Cassou et al. (1996). The results of the present theory were compared to simulations based on a recently introduced efficient implementation of Ewald summation method in a slab geometry, Ref [dos Santos, Girotto, and Levin, 2016]. The agreement between the theory and simulations is excellent.
The system under investigation can also be used to study interactions between colloidal particles Ohshima (1974a, b, c, 2006); Paillusson, Barbi, and Victor (2009). In the traditional Derjaguin, Landau, Verwey, and Overbeek (DLVO) theory it is assumed that charge neutrality is satisfied in the inter-particle region. On the other hand, the discussion in the present paper shows that at short separations between the colloidal surfaces charge neutrality is violated. In order to asses the effect of charge neutrality violation, we have used the contact value theorem to calculate the force between charged planar surfaces inside an electrolyte solution. We find that at large salt concentrations, and in some range of wall separations, the net force on each surface is attractive in the case of multivalent electrolytes. Using Derjaguin approximation, it should now be possible to construct the effective interaction potential between charged colloidal particles of radius . The work along these lines is currently in progress.
This work was partially supported by the CNPq, FAPERGS, CAPES, INCT-FCx, and by the US-AFOSR under the grant FA9550-12-1-0438.
Appendix A Fundamental Measure Theory
Here we discuss briefly the FMT. The basic assumption of the FMT is that the excess free energy of a hard sphere fluid has the general form
where the excess free energy density is a local functional of weighted densities , which are in turn defined as:
The simple form of Eq. (33) is strongly suggested by the leading term of in its low-density diagrammatic expansion. In this limit, the set of weight functions can be inferred from the deconvolution of the underlying Mayer functions . It is composed of the four scalar functions,
together with two vector weight functions ,
In the context of the Scaled-Particle Theory (SPT) these functions characterize the fundamental measures of hard spheres Rosenfeld (1988, 2002). Note that the corresponding weighted functions in (34) have dimensions of . Since the free energy density has clearly dimensions of , it follows that a it can be quite generally written as the following combination of weighted densities Rosenfeld (1990):
with the coefficients being all functions of the dimensionless weighted density . The next step in the FMT approach is to impose some conditions to be fulfilled in the limit of homogeneous density distribution, leading to a differential equation from which the coefficients can be determined. The condition used determine the version of FMT. In Rosenfeld’s original formulation, a SPT condition is applied in the limit where , resulting in the Percus-Yevick (PY) compressibility equation of state for the bulk fluid. The White-Bear version, on the other hand, requires that the (more accurate) Mansoori-Carnahan-Starling-Leland (MCSL) equation of state for hard-sphere mixtures be recovered in the bulk limit Roth et al. (2002); Yu and Wu (2002); Mansoori et al. (1971). The resulting White-Bear excess free energy density reads:
where . For a given set of density profiles , the White-Bear hard-sphere excess free energy follows then from the direct applications of Eqs. (33), (34) and (38). We finally note that the corresponding hard-sphere excess single-particle correlations used in Eq. (21) are given by:
Appendix B Thermodynamic route for the wall forces
We now show that the forces between the walls given by Eq. (32) also follows from the direct thermodynamic calculation of the wall osmotic pressure, regardless of the underlying approximations invoked for the free energy calculation. The osmotic pressure on the walls separated by a distance can be written as:
According to Eq. (11), the grand canonical potential for the two plate system is:
where represents the system size across the direction (see Fig. 1). The last term on the right-hand side is the electrostatic energy resulting from the direct wall-wall interaction. It has to be included since the remaining two terms only contain ion-ion and wall-ion interactions. The resulting osmotic pressure is:
Now, since the intrinsic free energy only depends on the ionic interactions, it can not have any explicit dependence on the particular external potential the particles are subjected to, and has therefore no explicit dependence on . All changes in come exclusively from the corresponding changes in the ionic profiles as is varied. Applying the usual rule for the derivative of a functional, we find:
In the last equality, the Euler-Lagrange equilibrium condition has been employed. With the above result, Eq. (42) for the osmotic pressure can be simplified to:
Since the chemical potentials have the purpose of fixing the ionic bulk concentrations – and these are in the present formulation kept constant as changes – the first derivative on the right-hand side can be set to zero. We further notice that the wall-ion interaction potential can be split into hard-core and electrostatic contributions, , which are given by Eqs. (8) and (28), respectively. It follows that the osmotic pressure can be written as:
We now turn to the calculation of each contribution on the right-hand side of (45) separately. First, it is convenient to rewrite the first integral as:
The so-called cavity functions are either zero at particle overlap with the hard walls or one otherwise. Similarly, the function goes to infinity at ion-wall overlap and to unity anywhere else. It follows from Eq. (28) that
Performing the straightforward differentiation with respect to provides:
Using the above result in Eq. (46) and performing the integration, we find:
where in the last equality the parity of the distribution functions was used. As for the electrostatic contributions in Eq. (45), we first notice that Eq. (8) can be used to write the ion-wall electrostatic potential as:
Differentiation with respect to provides:
The electrostatic contribution to the osmotic pressure then becomes:
Now, using the overall electroneutrality condition for the system as a whole, Eq. (1), the term in parenthesis can be clearly identified as , and the above relation can be simplified to:
which precisely recovers the result (32) for the force per unit of area exerted on the wall at . It is important to emphasize that no assumption has been made on other than its fulfillment the Euler-Lagrange equilibrium condition. The result above is, therefore, quite general and independent on the particular set of approximations employed in the construction of . This is to be contrasted with the situation in which electrolytes are present at only one single side of the surface, where then the contact condition does depend on the particular free energy functional Frydel and Levin (2012). When the charged surface lies in between two electrolytes, the mutual interaction between these electrolytes will also play a role on the net force acting upon the surface Lozada-Cassou and Yu (1996). The numerical derivative in Eq. (40) can in practice be performed as an accuracy check of the calculated ionic density profiles in Eq. (54), since the corresponding values are usually sensitive to numerical precision Wang, Liu, and Neretnieks (2011b).
- Levin (2002) Y. Levin, Rep. Prog. Phys. 65, 1577 (2002).
- Linse and Lobaskin (1999) P. Linse and L. Lobaskin, Phys. Rev. Lett. 83, 7208 (1999).
- Ŝamaj and Trizac (2011) S. Ŝamaj and E. Trizac, Phys. Rev. Lett. 106, 078301 (2011).
- Pienegonda, Barbosa, and Levin (2005) S. Pienegonda, M. C. Barbosa, and Y. Levin, Europhys. Lett. 71, 831 (2005).
- Fernandez-Nieves et al. (2005) A. Fernandez-Nieves, A. Fernandez-Barbero, F. J. de las Nieves, and B. Vincent, J. Chem. Phys. 123, 054905 (2005).
- Likos (2001) C. N. Likos, Phys. Rep. 348, 267–439 (2001).
- Zhang and Zhao (2009) L. L. Zhang and X. S. Zhao, Chem. Soc. Rev. 38, 2520â2531 (2009).
- Gouy (1910) G. Gouy, J. Phys. 9, 457 (1910).
- Chapman (1913) D. L. Chapman, Philos. Mag. 6, 475 (1913).
- Stern (1924) O. Stern, Z. Elektrochem. 30, 508 (1924).
- Trizac, Bocquet, and Aubouy (2002) E. Trizac, L. Bocquet, and M. Aubouy, Phys. Rev. Lett. 89, 248301 (2002).
- Fedorov and Kornyshev (2014) M. V. Fedorov and A. A. Kornyshev, Chem. Rev. 114, 2978â3036 (2014).
- Paillusson and Trizac (2011) F. Paillusson and E. Trizac, Phys. Rev. E. 84, 011407 (2011).
- Wu et al. (2011) J. Wu, T. Jiang, D. en Jiang, Z. Jin, and D. Henderson, Soft Matter 7, 11222 (2011).
- Frydel and Levin (2012) D. Frydel and Y. Levin, J. Chem. Phys. 137, 164703 (2012).
- Attard (1996) P. Attard, Adv. Chem. Phys. 92, 1–159 (1996).
- Yeh and Berkowitz (1999) I. Yeh and M. L. Berkowitz, J. Chem. Phys. 111, 3155 (1999).
- Arnold, de Joannis, and Holm (2002) A. Arnold, J. de Joannis, and C. Holm, J. Chem. Phys. 117, 2496 (2002).
- dos Santos, Girotto, and Levin (2016) A. P. dos Santos, M. Girotto, and Y. Levin, J. Chem. Phys. 144, 144103 (2016).
- Henderson and Blum (1978) D. Henderson and L. Blum, J. Chem. Phys. 69, 5441 (1978).
- Henderson, Blum, and Smith (1979) D. Henderson, L. Blum, and W. R. Smith, Chem. Phys. Lett. 63, 381 (1979).
- Vericat, Blum, and Henderson (1982) F. Vericat, L. Blum, and D. Henderson, J. Chem. Phys. 77, 5808 (1982).
- Ballone, Pastore, and Tosi (1986) P. Ballone, G. Pastore, and M. P. Tosi, J. Chem. Phys. 85, 2943 (1986).
- Kjellander and Marčelja (1986) R. Kjellander and S. Marčelja, Chem. Phys. Lett. 127, 402–407 (1986).
- Kjellander (1988a) R. Kjellander, J. Chem. Phys. 88, 7129 (1988a).
- Kjellander (1988b) R. Kjellander, J. Chem. Phys. 88, 7138 (1988b).
- Kjellander et al. (1992) R. Kjellander, T. Åkesson, B. Jönsson, and S. Marčelja, J. Chem. Phys. 97, 1424 (1992).