Electrostatic correlations on the ionic selectivity of cylindrical membrane nanopores

Electrostatic correlations on the ionic selectivity of cylindrical membrane nanopores

Sahin Buyukdagli111email: Sahin.Buyukdagli@iri.univ-lille1.fr and T. Ala-Nissila222email: Tapio.Ala-Nissila@aalto.fi Institut de Recherche Interdisciplinaire USR3078 CNRS and Université Lille I, Parc de la Haute Borne, 52 Avenue de Halley, 59658 Villeneuve d’Ascq, France
Department of Applied Physics and COMP Center of Excellence, Aalto University School of Science, P.O. Box 11000, FI-00076 Aalto, Espoo, Finland
Department of Physics, Brown University, Providence, Box 1843, RI 02912-1843, U.S.A.
July 7, 2019

We characterize the role of electrostatic fluctuations on the charge selectivity of cylindrical nanopores confining electrolyte mixtures. To this end, we develop an extended one-loop theory that can account for correlation effects induced by the surface charge, nanoconfinement of the electrolyte, and interfacial polarization charges associated with the low permittivity membrane. We validate the quantitative accuracy of the theory by comparisons with previously obtained Monte-Carlo simulation data from the literature, and scrutinize in detail the underlying forces driving the ionic selectivity of the nanopore. In the biologically relevant case of electrolytes with divalent cations such as in negatively charged nanopores, electrostatic correlations associated with the dense counterion layer in the channel result in an increase of the pore coion density with the surface charge. This peculiarity analogous to the charge inversion phenomenon remains intact for dielectrically inhomogeneous pores, which indicates that the effect should be observable in nanofiltration membranes or DNA-blocked nanopores characterized by a low membrane permittivity. Our results show that a quantitatively accurate consideration of correlation effects is necessary to determine the ionic selectivity of nanopores in the presence of electrolytes with multivalent counterions.


I Introduction

Electrostatic forces induced by the nanoscale confinement of charged liquids are at the origin of various industrial and biotechnological applications. A prominent example is the water desalination process that aims at separating salt from sea water by making use of the Donnan and dielectric exclusion mechanisms yar1 (); yar2 (); Lang (). The coupling of electrostatic and hydrodynamic forces is also the basis of nanofluidics separation technics that consist in controlling the flow of charged liquids thorough nanochannels. Although a newborn science, nanofluidics has already found many important applications from electromechanical energy conversion nano1 () to DNA sequencing nano2 () and polymerase chain reaction in lab-on-a-chip devices nano3 (). In order to predict and control the functioning of these devices, it is an important task to develop electrostatic formulations of confined charged liquids that can handle ion-ion and ion-substrate correlations in an accurate way; a challenge that remains to be met at present.

One of the major limitations of the current nanofluidic transport theories is clearly the underlying dielectric continuum electrostatics. Field theoretic formulations incorporating the charge structure of solvent dun1 (); orland1 (); buyukdagli13 (); jcpPOL () or fluctuating polyelectrolytes dun2 () have been previously proposed. The additional weak point of nanofluidic approaches is their mean-field (MF) nature. More precisely, these formulations that couple hydrodynamic transport equations with the Poisson-Boltzmann (PB) equation neglect electrostatic correlations, which are known to be highly relevant for charged nanopores confining electrolytes of general composition mcasy (); mcmix (). On the side of the artificial nanofiltration studies, ion rejection models yar1 (); yar2 () that make use of correlation corrected formulations such as the electrostatic self-consistent (SC) equations SCrussian (); netzvar () have been so far one step ahead. However, uncontrolled approximations involved in the derivation of the SC equations and in their solution in cylindrical pore geometries have obscured the quantitative precision of these models.

The standard way to properly account for ionic correlation effects consists of the one-loop expansion of the electrolytic free energy around the PB theory. Unfortunately, the unavailability of analytical solutions for the one-loop level electrostatic equations in curved geometries have restricted the proper consideration of ionic correlation effects to simple planar systems, whereas interfacial curvature effects have been mostly considered at the MF level anan1 (); anan2 (); anan3 (). At one-loop order, electrostatic correlation effects on the interaction energy of charged plates were considered within a semiclassical approximation by Podgornik and Zeks PodWKB () and in terms of special functions by Attard et al. attard (). Correlation induced modifications of ion densities at charged planes were also investigated at the same perturbative order in Refs. netzcoun (); 1loop (); jcp2 () or using strong-coupling expansion technics in Refs. st1 (); st2 (); st3 (); st4 (). However, ionic correlations in cylindrical geometries have been mostly studied within the DH approximation cyl1 (); cyl2 (); cyl3 (), which is more limited than the one-loop theory, and valid exclusively for neutral or very weakly charged cylinders where the surface charge induced electrostatic potential is much lower than the thermal energy .

The main advantage of the SC approach over the one-loop theory is its ability to deal with dielectrically inhomogeneous interfaces where the one-loop expansion is known to fail David (); jcp2 (). Different works introduced approximative solutions of the SC equations beyond one-loop level in slit systems hatlo (); pre (); japan (). Using a midpoint approximation that consists of replacing the local ionic self-energy by its value in the mid-pore, Yaroschuk solved the SC equations in cylindrical nanopores yar1 (). Then, within a restricted variational approach based on a uniform screening parameter in the pore, Buyukdagli et al. improved the solution of Yaroschuk and showed that ion transport through nanoscale pores is characterized by an ionic conductor-insulator (CI) transition PRL (); jcp1 (). Finally, Buyukdagli et al. recently developed a systematic numerical scheme for the exact solution of the SC equations in slit geometries jcp2 (). Comparisons with Monte-Carlo (MC) simulations showed that this SC scheme can handle surface dielectric discontinuities with a better accuracy than the DH approach.

Following these observations, we introduce in the present work an extended one-loop approach that is based on the truncation of the SC equations in such a way that the scheme reduces to the one-loop theory for dielectrically homogeneous systems, but considers surface polarization effects associated with low permittivity membranes in a self-consistent fashion. The latter point is indeed the most important advantage of our approach over numerical simulation schemes, since MC simulation techniques have been so far unable to deal with interfacial dielectric discontinuities in cylindrical pores. The article is organized as follows. We review in Section II.1 the physical framework of the SC formalism, and briefly introduce the extended one-loop theory whose detailed derivation is reported in Appendices B and C. Then, in Section III.1, we make use of the present approach in order to analyze the MC simulation data of Refs. mcasy (); mcmix () for the partition of electrolyte mixtures in dielectrically homogeneous cylindrical pores. Finally, Section III.2 is devoted to dielectric discontinuity effects on the ionic selectivity of the pore, and our main results are summarized in the Conclusion part.

Ii Theory

ii.1 Self-consistent formalism

We review in this part the SC formalism introduced in Refs. SCrussian (); netzvar (). The geometry of the membrane nanopore is depicted in Fig. 1. It consists of a cylindrical pore of infinite length and radius , and a negative wall charge distribution with . The pore is in contact with an ion reservoir at the extremities, which is imposed by assuming the chemical equilibrium condition between the charges in the pore and the reservoir jcp1 (). Furthermore, the electrolyte is composed of ion species, with each species of reservoir concentration and valency . The number density of each ionic species is given by


where the wall potential that imposes the ionic confinement inside the cylinder is introduced as if , and for . Moreover, the function in Eq. (1) stands for the electrostatic potential generated by the fixed membrane charge at the pore wall or ionic charge excesses induced by charge separation effects. Then, the ionic self-energy is defined in the form of an equal point Green’s function


with the electrostatic Green’s function corresponding to the potential induced by an ion located at . In an ion free bulk reservoir, the latter reduces to the simple Coulomb potential . Furthermore, the coefficient Å in Eq. (2) stands for the Bjerrum length at ambient temperature , and is the DH screening parameter. Indeed, the self-energy introduced in Eq. (2) accounts for the spatial variations of the chemical potential of an ion resulting from the deformation of its screening cloud close to boundaries of the system (solvation energy) and its interaction with surface polarization charges (image charge forces).

Figure 1: (Color online) Schematic representation of the pore geometry. The pore radius is , the net negative surface charge , and the pore and membrane permittivities are respectively and .

The external potential and the Green’s function in Eqs. (1) and (2) are solutions of the electrostatic self-consistent equations SCrussian (); netzvar (),


where the dielectric permittivity function accounting for the variations of the background permittivity at the boundaries of the cylindrical pore is given by , with and the water and membrane permittivites, respectively. Moreover, stands for the elementary charge, and the Boltzmann constant. One notes that Eq. (3) is an extended Poisson-Boltzmann (PB) equation for the electrostatic potential . Furthermore, Eq. (II.1) corresponds to a generalized Debye-Huckel (DH) equation for the Green’s function . The difference from the usual PB and DH equations is respectively due to the non-uniform ion density in Eq. (II.1), and the presence of the ionic self-energy term in Eq. (1) that introduces in Eq. (3) an inhomogeneous screening of the external potential. We also note that the bulk limit of Eqs. (3)-(II.1) was investigated in Ref. pre (). Indeed, it was shown that for monovalent ions with concentration M, the solution of the variational equations yields the DH limiting law, whereas for higher concentrations, the electrolytic free energy becomes unstable. Although this instability can be avoided with an ultraviolet cut-off, the regularization is unnecessary since we have previously showed that the instability regime lies well beyond the ion density range where the SC formalism is quantitatively accurate jcp2 ().

To characterize the correlation effects induced by the confinement on the charge selectivity of the pore, it is helpful to identify the key system parameters. To this end, we consider a dielectrically homogeneous pore confining a symmetric electrolyte composed of two ion species of valencies , with . Rescaling the distances in Eqs. (3)-(II.1) according to , and renormalizing the external potential and the Green’s function as and , the SC equations take the form


The rescaled equations (5)-(II.1) show that for dielectrically homogeneous pores and symmetric electrolytes, the ion density in the pore is solely characterized by three parameters. These are the reduced pore radius , the coupling parameter characterizing the strength of the electrostatic potential fluctuations around the MF potential, and the parameter , where stands for the Gouy-Chapman length. Namely, the parameter measures the ratio of the counterion layer thickness at the charged pore wall to the screening cloud radius around a central ion in the bulk reservoir. In Section III.1 on electrolytes in dielectrically homogeneous pores, it will be shown that electrostatic correlation effects are characterized by an interpolation between the parameter regime corresponding to a compact interfacial counterion layer and the regime associated with a diffuse counterion partition close to the pore wall.

In the present work, we will generalize to dielectrically discontinuous pores the one-loop approach developed in Refs. netzcoun (); 1loop (); jcp2 () for dielectrically uniform systems. Our strategy consists in truncating the SC Eqs. (3)-(II.1) in such a way that the extended scheme reduces to the one-loop theory for dielectrically uniform pores with . The details of the derivation of the extended approach that we summarize in Section II.2 are reported in Appendix C.

ii.2 Truncating SC equations

The one-loop theory of a symmetric electrolyte in contact with a dielectrically continuous charged plane was shown in Ref. jcp2 () to follow from the linearization of the SC Eqs. (3)-(II.1) in terms of the electrostatic Green’s function . The extension of the one-loop approach to a general electrolyte mixture is presented in Appendix B. In the present work, we also introduce a generalization of the one-loop approach to dielectrically discontinuous systems where the singular nature of the self-energy on the pore wall does not allow a simple Taylor expansion of Eq. (1David (); jcp2 (). Our strategy that we detail in Appendix C consists in treating in Eqs. (3)-(II.1) the part of the Green’s function induced by the dielectric discontinuity self-consistently, and expanding the rest resulting from solvation forces perturbatively.

To derive the truncated equations, we formally split the self-energy into a solvation and an image charge contribution as , where the potential accounts for the interaction of ions with the interfacial polarization charges, and the potential incorporates the solvation energy induced by the deformation of the ionic cloud in the pore. Then we insert this expansion into Eqs. (3)-(II.1), and linearize the latter in terms of the solvation energy . This truncation results in an expansion of the external potential in the form , where the component of the potential and the Green’s function satisfy the differential equations


and the correction term to the external potential is given by


with the auxiliary number and excess charge densities


One sees that Eq. (7) is a modified PB equation that accounts for the variations of the external potential by the image charge forces. Furthermore, Eqs. (9) and (11) indicate that the term brings corrections from solvation forces. The explicit form of the potentials and as well as the iterative solution scheme of Eqs. (7)-(9) is explained in Appendix C. We finally note that in the limit , Eqs. (7)-(9) reduce to the one-loop equations derived in Appendix B.

Iii Results

iii.1 Dielectrically homogeneous pores

We characterize in this section ionic correlation effects on the selectivity of dielectrically homogeneous pores () confining electrolytes of general mixture. The theoretical ion density results will be compared with MC simulation data from Refs. mcasy (); mcmix () in order to establish the quantitative accuracy of the present theory. The transparency of our theoretical scheme will also allow to probe in detail the underlying physics behind the simulation results. Section III.1.1 is devoted to correlation effects on the partition of asymmetric electrolytes composed of two species, and we will consider in Section III.1.2 the more complicated case of mixed electrolytes.

iii.1.1 Asymmetric electrolytes

In this part, we focus on the MC simulation results of Ref. mcasy () for the partition of asymmetric solutions and scrutinize in detail the underlying electrostatic interactions. The reservoir contains an asymmetric electrolyte composed of two ion species. We first consider in Fig. 2(a) the case of divalent coions and monovalent counterions ( and ), with the pore radius nm and charge . The figure displays the rejection rates defined as , where the ionic partition function corresponds to the pore averaged ion density,


with the local ion densities given by (see Appendix B)


It is seen that over the bulk density regime M, the MC data (red dots) and the one-loop result (solid black curves) in Fig. 2(a) predict a slightly higher rejection rate than the MF theory (dashed black curve). The same effect is also illustrated in Fig. 3(a) that displays the coion rejection rates against the surface charge at fixed bulk concentration M, and for pore radii nm and nm. To understand the underestimation of the rejection rates by the MF theory, one should note that for monovalent counterions, the surface charge densities in Fig. 3(a) correspond to the parameter regime characterized by a diffuse counterion layer next to the pore wall. In this regime where the ion free membrane is responsible for a charge screening deficiency in the vicinity of the charged wall, ions feel a repulsion from the wall towards the mid-pore area where they are more efficiently screened. Thus, correlations decrease the MF density of both anions and cations in the pore. The positive ionic self-energy embodying this repulsive force is reported in Fig. 3(b) (solid black curve).

Figure 2: (Color online) Coion rejection rates against the anion density for (a) divalent anions and monovalent counterions , and (b) monovalent anions and divalent counterions . The pore charge is , pore radius nm, and the membrane permittivity . The dashed and solid black curves are respectively the MF and one-loop results, and the MC data in (a) and (b) denoted by the red circles are respectively taken from Tables III and IV of Ref mcasy ().

An interesting point to be noted in Fig. 3(b) is the increase of the ionic self-energy in the mid-pore with an increase of the surface charge from to . The positive energy barrier in the mid-pore area that appears at finite charges is indeed the second correlation effect behind the deviation of the MF curves from the one-loop and MC results in Fig. 3(a). This peculiarity can be understood in terms of the effective screening parameter

Figure 3: (Color online) (a) Coion rejection rates against the pore charge for divalent anions , monovalent counterions , with the bulk anion density M, and membrane permittivity . The curves and symbols have the same meaning as in Fig. 2. MC data were taken from Table III of Ref mcasy (). (b) Ionic self energy and (c) the effective screening parameter of Eq. (14) for the same system parameters as in (a) and the surface charges given in the legend.

that determines the local screening of the one-loop potential in Eq. (B), with the function defined in Eq. (10). In Fig. 3(c), it is seen that while moving from the charged wall towards the mid-pore, the local screening parameter drops below the bulk one. This originates from the exclusion of divalent coions by the surface charge, which dominates the monovalent counterion excess in the pore. The corresponding charge screening deficiency with respect to the bulk reservoir results in an electrostatic energy barrier for ionic penetration. Then, by further increasing the surface charge from to , one gets into the second regime where the counterions form a dense layer in the vicinity of the charged wall. The resulting interfacial ionic screening excess with respect to the bulk reservoir (see Fig. 3(c)) leads to a negative ionic self-energy close to the pore wall, although the mid pore area is still characterized by the strong coion deficiency resulting in the positive branch of the self-energy (see Fig. 3(b)). As a result, the one-loop rejection rates in Fig. 3(a) drop below the MF result for large pore charges. One can also notice in Fig. 3(a) the good agreement between the one-loop theory and the MC data, which is remarkable if one considers the large coupling parameter for divalent anions at the bulk density M.

Figure 4: (Color online) Coion rejection rates for pore radius (a) nm and (b) nm against the pore charge for monovalent coions , divalent counterions , with the bulk anion density M, and the membrane permittivity . The curves and symbols have the same meaning as in Fig. 2. MC data were taken from Table IV of Ref mcasy (). (c) Ionic self energy and (d) effective screening parameter Eq. (14) for the same system parameters as in (a) and the surface charges given in the legend.
Figure 5: (Color online) Ion density profiles for the electrolyte mixture and in the nanopore with (a) vanishing surface charge and radius nm, and (b) finite surface charge and radius nm. Bulk ion densities are M, M and M, and membrane permittivity is . The curves and symbols have the same meaning as in Fig. 2. MC data in (a) and (b) were respectively taken from Figs. 4 and 12 of Ref mcmix (). (c) External potential (inset) and ionic self energy (main plot) for the same system parameters as in (a) (dashed-dotted curves) and (b) (solid curve).

We next consider an asymmetric electrolyte composed of monovalent coions and divalent counterions. Fig. 2(b) displays for this case the anion rejection rates against the bulk density at the pore charge . The inspection of the figure shows that the one-loop result remains very close to the MC data, significantly improving over the MF result. Then, it is seen that the latter now overestimates the MC and one-loop rejection curves. This effect is also illustrated in Figs. 4(a) and (b) over the whole surface charge regime . To understand this peculiarity, one should note that for divalent counterions, the pore is characterized by a significant counterion excess even for weak surface charges. The resulting screening excess displayed in Fig. 4(d) is shown in Fig. 4(c) to yield a significantly negative ionic self-energy on the order of the thermal energy , which is monotonically lowered with increasing surface charge. As a result of this pure correlation effect, the coion rejection rates in Fig. 4(a) obtained from MC and one-loop approaches are seen to reach a peak at a characteristic surface charge, and to decrease for higher pore charges where the screening excess in the pore takes over the Donnan exclusion. As noted in Ref. mcasy (), this feature may be the precursor of the negative rejection of electrolytes with multivalent counterions observed in ion rejection experiments with porous glass exp0 (). Thus, the enhancement of coion densities with increasing surface charge clearly deserves further investigation, and we will reconsider the effect in Section III.2.3 in further detail.

Finally, in Figs. 4(a) and (b), the one-loop theory is seen to exhibit a quantitative agreement with MC data exclusively at low surface charges. Namely, the turning point where the coion density starts to increase with the surface charge is seen to be underestimated by the one-loop approach, and the deviation increases beyond this value. Although the discrepancy should be mainly due to electrostatic correlation effects, we also expect the hard-core (HC) interactions absent in the one-loop theory to be partly responsible for the disagreement.

iii.1.2 Electrolyte mixtures

We consider in this part the simulation results of Ref. mcmix () for the partition of an electrolyte mixture composed of a symmetric and an asymmetric electrolyte, such as and . Fig. 5(a) compares the one-loop ion densities with the MC data in a narrow pore with radius nm and vanishing surface charge . The ion densities in the reservoir are M, M and M, which are marked in the plots by the dotted horizontal curves. First of all, it is seen that despite the large coion density and strong confinement, the theory agrees well with MC simulations within the numerical uncertainties of the data, except in the vicinity of the pore wall where the theory predicts a stronger ion depletion. This deviation may be due to the fact that our formalism neglects the interfacial ion-free layer associated with the finite ionic radius. This layer could be easily incorporated into the Green’s function computed in Appendix A as it was done in Refs lev2 (); lev3 () in planar geometry and in Ref jcp1 () for cylinders, but we leave this improvement for a future work for the sake of simplicity. We also note in passing that for ions in a neutral pore, the MF theory yields the ideal gas behavior , which is denoted in Figs. 5 (a) and (b) by the horizontal lines.

Two features to be noted in Fig. 5(a) are an overall depletion of and ions, and a weak pore excess of ions. To explain these points, we report in Fig. 5(c) the ionic self-energy (dashed-dotted curve in the main plot) and the external potential (inset). It is seen that the deformation of the screening cloud around the ions results in a positive self-energy, i.e. an electrostatic energy barrier depleting the monovalent anions and divalent cations from the pore. Furthermore, as a consequence of the stronger exclusion of divalent cations with respect to monovalent anions, the bulk electroneutrality is locally perturbated in the pore. The corresponding charge separation effect leads in turn to a weakly negative external potential (see the inset of Fig. 5(c)), resulting in a excess in the mid-pore area.

We next consider in Fig. 5(b) a larger pore of radius nm and surface charge . The bulk ion densities are the same as in Fig. 5(a). The theoretical density curves are again seen to yield a reasonable agreement with MC simulation data that exhibit large uncertainties. In this case, the counterion excess induced by the surface charge at the pore wall (see Fig. 5(b)) is shown in Fig. 5(c) to yield a weakly negative ionic self-energy, except in the vicinity of the pore wall where the screening defficiency takes over and results in a positive branch of the self-energy. In Fig. 5(b), one sees that the repulsive branch leads to an ionic depletion layer close to the wall, and its competition with the surface charge results in a concentration peak for counterions. Moreover, in the same figure, the attractive branch of the self-energy is shown to slightly increase the MF prediction of cation densities. In the Section III.2, we will incorporate into this picture surface polarization charges resulting from the dielectric discontinuity between the pore and the low permittivity membrane.

iii.2 Dielectrically inhomogeneous pores

Figure 6: (Color online) Ion density against the distance from the neutral pore wall for a symmetric electrolyte of two monovalent ion species. The solid red and blue curves are respectively MC simulation data from Ref. jcp2 () and the DH density for ions at a planar dielectric interface, and the square symbols and solid black curve respectively mark for a cylinder of radius nm the ion density profile from the DH approximation and the SC scheme introduced in Appendix C. The bulk salt concentration is M, and the permittivities are and .

We consider in this part ionic partitions in cylindrical nanopores separating the electrolyte from a membrane medium with a low dielectric permittivity . The SC solution of the electrostatic closure equations (7)-(9) for dielectrically inhomogeneous pores is explained in Appendix C. For this more general case, the ion density is given by


where the external potential is the solution of the MPB equation (51), the correction to the external potential is obtained from Eq. (59), and the solvation potential and the image-charge energy resulting from the dielectric discontinuity at are given by Eq. (C).

iii.2.1 Comparison with MC simulations

We note that the present SC scheme that treats solvation forces in a perturbative fashion differs from the approach developed in Ref. jcp2 () for ions at dielectric planes. Although there are no available MC simulation data for ion densities in dielectrically discontinuous cylinders, we can still test the quantitative validity of the present theory for cylinders with a large radius by comparison with MC simulations at planar dielectric interfaces. To this end, we compare in Fig. 6 the MC data from Ref. jcp2 () for a symmetric electrolyte at a neutral dielectric plane (red curve), with the predictions of the present scheme for a cylinder of radius nm (black curve). The bulk ion concentration is M, and the membrane permittivity . In order to illustrate the corrections beyond the WC regime, we also report the DH density profiles for a planar interface jcp2 () (solid blue curve) and a cylinder of radius nm (square symbols). The DH density profile for the cylinder is computed from the relation , with the DH potential given by Eq. (26) of Appendix A in the limit .

First of all, in Fig. 6, one notes that the DH density profiles for the planar interface and the cylindrical pore coincide, which indicates that curvature effects become negligible for the pore radius nm. Then, it is seen that the DH theory overestimates the MC result for the ion density. The failure of the DH theory for large densities in the range M was shown in Ref. jcp2 () to result from the unability of the WC theory to account for the local variations of the ionic screening. Finally, the ion densities computed with the present SC scheme is seen in Fig. 6 to exhibit a very good agreement with the MC data. This point shows that the present approach can accurately handle image charge forces beyond the WC regime. Encouraged by this observation, we make use of our SC scheme to reconsider in the next Section III.2.2 the ionic CI transition that we had discovered in Refs. PRL (); jcp1 () within a restricted variational scheme.

Figure 7: (Color online) Ionic partition coefficient against the pore radius for a symmetric electrolyte of two monovalent ion species with reservoir density M and various values of the membrane permittivity . The transition points where the ion density drops to zero are marked by dotted vertical lines.

iii.2.2 Conductor-insulator transition

We illustrate in Fig. (7) ionic partition coefficients computed within the present approach from Eqs. (12) and (15) against the pore radius for different values of the membrane permittivity. In the case of a low permittivity membrane associated with a strong dielectric jump at the pore wall (solid black curve), the partition coefficient is seen to gradually decrease with the pore size, until it sharply drops from to zero at the pore radius 10 Å. We note that within the restricted variational scheme based on a uniform pore screening parameter, we had already observed in Refs. PRL (); jcp1 () the same ionic CI transition for cylindrical pores. The appearance of this effect within the general SC scheme of the present study shows that the transition is not a simple artefact of the restricted variational choice in Refs. PRL (); jcp1 ().

Then, with an increase of the membrane permittivity from to , the jump in the ion density is seen to be significantly reduced, and the transition becomes a smooth one for the larger permittivity . In Refs. PRL (); jcp1 (), it was shown that the CI transition results from a competition between two opposing effects. On the one hand, the deformation energy of the screening cloud around ions associated with solvation forces is amplified with ionic penetration and thus favors ionic exclusion from the pore. On the other hand, the screening of image interactions that lowers the free energy of the electrolyte favors ionic penetration. The smoothing of the transition stems from the weakening of the competition between these two effects as a consequence of the reduction of image charge forces with increasing . Indeed, the rounding of the transition with a small reduction of the dielectric discontinuity indicates that the effect may be difficult to observe in nanoscale membrane pores where the water permittivity is also expected to be reduced with respect to the reservoir permittivity.

(a) (b)

Figure 8: (Color online) (a) Coion partition coefficient against the pore charge and (b) the effective screening parameter Eq. (14) for an asymmetric electrolyte composed of monovalent coions and divalent counterions . The bulk anion density is M and the pore radius nm. In (b), the surface charge values are (dashed curves) and (solid curves), and the membrane permittivities (blue curves) and (black curves).

Finally, in Fig. 7, the comparison of the lower curves with the result for a dielectrically homogeneous pore shows that regardless of the pore radius, the dielectric exclusion is clearly the dominant rejection mechanism in the nanopore. This is in line with recent nanofiltration experiments where image-charge forces were shown to play the key role in ion rejection from artificial membrane nanopores Lang ().

iii.2.3 Enhancement of coion densities with surface charge

At the end of Section III.1.1 on electrolytes composed of monovalent coions and divalent counterions, it was shown that an increase of the surface charge beyond a characteristic value results in a decrease of the coion rejection rates (see Figs. 4(a) and (b)). This seemingly counterintuitive effect displayed in Fig. 8(a) in terms of the coion partition function (black curve) was shown to result from the strong counterion attraction into the nanopore, transforming the latter into a medium with high screening ability (see the local pore screening parameters displayed by black curves in Fig. 8(b)).

Figure 9: (Color online) Diagram characterizing the enhancement of coion densities with the surface charge for an asymmetric electrolyte composed of monovalent coions and divalent counterions , with the membrane permittivity and the pore radius nm.

Then, in Fig. 8(a) where we illustrate the anion partition functions for dielectrically discontinuous pores with and , it is seen that the surface charge induced enhancement of coion densities is also present for low permittivity membrane nanopores characterized by a strong ion rejection. This peculiarity can be explained in terms of the electroneutrality condition in the pore. More precisely, in Fig. 8(b) where we reported the local screening parameter for , the dielectric exclusion of counterions from the pore wall is seen to be compensated by a local counterion excess in the mid-pore area. The latter effect originates from the global electroneutrality condition that fixes the total number of coions and counterions over a crossection of the charged pore. As a result of this mechanism, the enhancement of the coion attraction with the pore charge survives for dielectrically inhomogeneous pores. Thus, the effect should be observable in real membrane nanopores characterized by a low static dielectric permittivity .

In order to determine the charge density regime where the coion density enhancement is expected, we plot in Fig. 9 the characteristic surface charge versus reservoir concentration that marks the boundary between the area for the surface charge induced anion decrease ( below the curve) and increase ( above the curve). We note that the transition curve also fixes the applicability limit of the WC DH theory that always predicts an enhanced coion rejection with surface charge, i.e. . The diagram in Fig. 9 shows that increasing the ion density from M to M, the characteristic pore charge where the effect is expected to come into play is moved from to . It should be noted that these values correspond to considerably weak charge densities, and most of the charged nanoscale systems are actually located in the parameter regime above the transition curve of Fig. 9. For example, the characteristic surface charge of a DNA molecule in a pore or the wall charge of a polyethylene terephthalate membrane nanopore at are both on the order . This indicates that the effect scrutinized in this part is relevant to a large variety of biological and industrial nanopore systems, and charge fluctuation effects should be properly taken into account beyond the DH theory for an accurate determination of their functioning.

Iv Conclusion

In this work, we have characterized electrostatic correlation effects on the ionic selectivity of cylindrical nanopores. To this aim, we have developed an extended one-loop approach that can account for the charge correlations induced by the strong confinement in the cylindrical pore, the surface charge, and the interfacial polarization charges associated with the low dielectric permittivity of the membrane. We have confirmed the quantitative accuracy of the theory and determined its validity regime by making extensive comparisons with MC simulation data.

In Section III.1, we made use of the present theory in order to analyze MC simulation results of Ref. mcasy (); mcmix () for the partition of asymmetric and mixed electrolytes in dielectrically homogeneous pores. It was shown that for nanoscale pores with surface charges in the range , ionic correlation effects can be categorized into two parameter regimes. For electrolytes with monovalent counterions where the size of the surface counterion layer is larger than the screening cloud radius in the bulk, i.e. , charge screening deficiency close to the pore wall results in an electrostatic barrier decreasing the MF prediction for ion densities in the pore. We emphasize that this effect can be already taken into account in a qualitative way by the weak-coupling DH theory.

In the case of an electrolyte solution with divalent counterions such as where one gets into the second regime , the compact interfacial counterion layer responsible for a charge screening excess in the nanopore was shown to attract both negatively and positively charged particles, increasing their pore density above the MF prediction. This second parameter regime is precisely not covered by the DH theory since the latter is unable to account for the non-uniform screening of the electrostatic potential. Most importantly, we found that beyond a characteristic value of the pore charge where the screening ability of the nanopore overcomes the Donnan rejection mechanism, the coion density rises with the surface charge.

In Section III.2, we incorporated into this picture image charge interactions associated with the dielectric discontinuity between the pore and the membrane. First of all, it was shown that the general SC approach introduced in this work yields the same ionic CI transition as the one already observed by the restricted variational approach of Refs. PRL (); jcp1 (). This suggests that the transition is not an artefact of the restricted trial potential used in the previous variational theory. However, we also showed that the transition is smoothed if the membrane permittivity exceeds , which indicates that the effect may be difficult to observe in membrane nanopores where the strong confinement is also expected to reduce the pore dielectric permittivity. Thus, in order to consider the case of biological pores with subnanometer radius, future works should consider the electrostatic nature of the solvent beyond the dielectric continuum approximation. Moreover, the enhancement of coion densities with surface charge in the presence of divalent counterions was shown to survive for dielectrically inhomogeneous pores. Thus, this effect should be observable in biological and artificial nanopores characterized by a low membrane permittivity. Finally, for electrolytes with submolar bulk concentrations, the effect was shown to come into play at considerably weak surface charge densities on the order . This means that the surface charge induced anion density enhancement should be present in many biological and artificial pore systems where the characteristic surface charges are usually an order of magnitude above this limit. This observation confirms the significance of the present SC approach in the study of nanopores with multivalent counterions.

The approach developed herein presents important advantages over the existing theoretical tools of confined electrolytes. First, for dielectrically uniform pores, it allows to determine the partition of the ions in the cylindrical nanopore at one-loop level, which has been so far limited to single planar interfaces netzcoun (); 1loop (); jcp2 (). Second, the formalism can also accurately consider beyond the WC regime the induced polarization charges resulting from the dielectric jump at the pore wall, a complication that cannot be taken into account by MC simulation technics based on the image charge convention in cylindrical pores. We should also mention that the present theory is less time consuming than MC simulations, since in the most complicated case of a charged nanopore with a strong dielectric discontinuity and divalent counterions (see Fig. 8), the computation of the ion densities for a given parameter set requires less than 20 minutes on a single 2.7 GHz processor.

We would like to discuss as well the limitations of the present theory. As it was shown in Fig. 4, in the presence of divalent counterions, the quantitative agreement with MC results worsens for pore charges above . Since the theoretical results overestimate the coion attraction into the nanopore, we think that the absence of HC interactions in the theory may be partly responsible for this failure. The incorporation of HC interactions is also necessary in order to consider concentrated solutions above the characteristic density M where these interactions were shown to come into play even at simple planar interfaces jcp2 (). For lower ion densities where HC collisions play a minor role, one would expect the finite ion size to modify exclusively the ionic densities in the close vicinity of the pore wall associated with an ion-free Stern layer lev2 (); lev3 (). The latter can be easily incorporated into the electrostatic Green’s function computed in Appendix A, but we expect this improvement to weakly modify the pore average ion densities for channels with nanoscale radius. Moreover, one should note that unlike the biological channels having a finite length, the nanopore model considered in this work is an infinitely long cylinder. The finiteness of the cylinder length is expected to reduce the strength of correlation effects and particularly the dielectric exclusion effect induced by image charge interactions. In a future work, one could extend the present theory by accounting for the finite nanopore length in the same way as it was done in Ref. cyl1 () at the DH level. To conclude, we emphasize that in view of the importance of the cylindrical confinement geometry in nanosystems, the present formalism can find important applications from nanofluidics to artificial nanofiltration technology where a quantitatively accurate consideration of charge correlation and surface polarization effects is still missing.

We thank Ralf Blossey for a critical reading of our work and his valuable comments. SB gratefully acknowledges support from The Academy of Finland through its COMP CoE (no. 251748) and NanoFluid grants and from a postdoctoral grant through the french ANR blanc project ‘Fluctuations in Structured Coulomb Fluids’.

Appendix A Derivation of the reference potential in cylindrical coordinates

We derive in this Appendix the reference Green’s function that will be needed to solve the closure equations (7)-(9) derived in Appendices B and C. The confinement geometry corresponds to a cylinder of charge and radius (see Fig. 1). We choose the reference Green’s function in the form of a DH potential with a general but uniform screening parameter , and satisfying the differential equation


where we introduced the piecewise dielectric permittivity profile . Exploiting the cylindrical symmetry of the system and inserting the Fourier expansion of the potential


into the relation (16), the equation for the Fourier transformed Green’s function follows as


For the solution of the closure equations (7)-(9), we exclusively need the solution of Eq. (A) for the charge sources inside the cylinder, i.e. for . In this case, the homogeneous solution of Eq. (A) is indeed given by


where and stand respectively for the modified Bessel functions of the first and second kind, and we also introduced the auxiliary function .

The integration constants with in Eq. (19) are obtained by imposing the boundary conditions associated with the continuity of the potential and the displacement field. The first set of boundary conditions follows by imposing the continuity of the Green’s function at the position of the source ion and at the cylinder wall ,


The remaining two boundary conditions are in turn obtained by integrating Eq. (A) around and , which yields


Imposing the boundary conditions in Eqs. (20)-(A) to the homogeneous solution (19), the latter finally follows for in the form

where we defined the function


with the parameter and the function . Moreover, we note that the solution of the potential for follows by interchanging in Eq. (A) the variables and . Finally, the self energy associated with the Green’s function (A) follows from Eq. (2) as


For the solution of the SC equations in dielectrically discontinuous pores, we note that the self-energy in Eq. (26) can be separated into a solvation and an image-charge part as , where

with the auxiliary functions


One sees that in the limit , the function (A) and the image-charge contribution to the self-energy in Eq. (A) vanish.

Appendix B One-loop solution of SC equations

We explain in this Appendix the one-loop solution of the SC equations (3) and (II.1) for ions confined in a dielectrically homogeneous interface i.e. , with the surface charge distribution . The one-loop approximation consists of an expansion of the SC equations around the MF theory in terms of the Green’s function and the self-energy . To this end, we split the external potential into a MF and a fluctuating part as , where the MF potential is the solution of the PB equation


where we introduced the MF level ionic number density as


Linearizing the SC equations (3) and (II.1) in terms of the Green’s function and the excess potential , and making use of the MF equation (31), one obtains the following equations


where we introduced in Eq. (B) the excess charge density


With the use of the definition of the Green’s function


and the inverse of the one-loop Green’s function

Eq. (B) can be inverted as


In order to evaluate the potential (38), one has to compute the one-loop Green’s function by solving Eq. (B). Since this equation has no analytical solution in cylindrical coordinates, one has find the solution numerically. To speed up the numerical solution scheme, we opt to solve this equation by choosing as the reference potential the Donnan Green’s function, with the corresponding kernel that we introduce in the form


and the effective screening parameter


Indeed, the approximation in introducing Eq. (39) consisted in replacing the MF potential in Eq. (32) and (B) by the constant Donnan potential , which is solution of the equation


We note that Eq. (41) follows by neglecting in Eq. (31) the spatial variations of the potential, and integrating the rest of the terms over the cross-section of the channel.

Then, using Eqs. (36), (B), and (39), the relation (B) can be formally inverted as


where we introduced the excess number density


Comparing Eqs. (16) and (39), one notices that the Green’s function is obtained from the potential (A) computed in Appendix A by setting .

Moreover, accounting for the cylindrical symmetry and substituting into Eqs. (38) and (42) the Fourier expansion of the Green’s function in the form of Eq. (17), the one loop-potential and the Green’s function in Fourier basis takes the form


We finally note that the one-loop level ion densities given in Eq. (13) follow by expanding Eq. (1) to linear order in and .

To evaluate the one-loop potential correction in Eq. (44), one has to solve Eq. (B) by iteration. The first step consists in numerically solving the PB Eq. (31) with the boundary conditions


where the parameter is the Gouy-Chapman length. Then, at the first iterative step, the MF potential profile is injected into Eqs. (43) and (B), together with the Donnan potential obtained from the solution of Eq. (41), and the Donnan Green’s function from Eq. (A) with used as the input function instead of the function on the r.h.s. of Eq. (B). At the next iterative step, the obtained solution for is reinjected into the r.h.s. of Eq. (B), and this cycle is continued until self-consistency is achieved. In the end, the converged solution for the Fourier transformed potential is used to evaluate the self-energy that follows from Eq. (B) in the form

We also note that in Eq. (B), the self-energy follows from Eq. (A) by setting . Then, the Fourier transformed potential is used with the self-energy (B) in Eq. (44) in order to compute the one-loop correction to the external potential . Finally, the obtained potentials , , and are inserted in Eq. (13) in order to evaluate the ion densities.

Appendix C Extended one-loop theory for dielectrically inhomogeneous pores

In this Appendix, we generelize the one-loop scheme introduced in Appendix B to dielectrically discontinuous systems. To this end, we first split the external potential and the self-energy into a solvation contribution and a singular part resulting from image-charge interactions as


where we choose the component as the solution of the modified PB (MPB) equation,


with the number density


The image-charge potential in Eqs. (50)-(52) will be introduced below.

Inserting the potentials in Eqs. (49)-(50) into the SC Eqs. (3)-(II.1), linearizing the latter in terms of the potentials and , and making use of the MPB equation (51), one gets the following relations,


where we introduced in Eq. (C) the excess charge density


Identifying from Eq. (C) the kernel associated with the Green’s function