Breakdown of electroneutrality in nanopores

TL;DR

This paper investigates how electroneutrality can break down in ultra-narrow nanopores, revealing that classical continuum models fail under extreme confinement and leading to new insights into ionic transport behavior.

Contribution

It develops a theoretical framework explaining electroneutrality breakdown in one-dimensional nanopores with 3D electrostatics, clarifying experimental conductance scaling.

Findings

Electroneutrality can break down in narrow nanopores due to large screening lengths.

Ionic transport behavior deviates significantly from classical models under extreme confinement.

Theoretical predictions match experimental conductance scaling in carbon nanotubes.

Abstract

Ion transport in extremely narrow nanochannels has gained increasing interest in recent years due to its unique physical properties, and the technological advances that allow us to study them. It is tempting to approach this regime with the tools and knowledge developed for wider microfluidic devices and use continuum models like the Poisson-Nernst-Planck equation. However, it turns out that some of the most basic principles we take for granted in a large system, such as electroneutrality, can breakdown under extreme confinement. We show that in a truly one-dimensional system, interacting with three-dimensional electrostatic interactions, the screening length is exponentially large in ionic spacing, and can easily exceed the macroscopic length of a nanotube. Without screening, electroneutrality is broken, and ionic transport can behave in a completely different way. In this work, we build a theoretical framework for electroneutrality breakdown in a one-dimensional nanopore and show how it provides an elegant interpretation for the peculiar scaling observed in experimental measurements of ionic conductance in carbon nanotubes.

Equations

$$\{\begin{array}{ll}{\epsilon }_{\mathrm{in}}{\nabla }^2{\varphi }_{\mathrm{in}}(\mathbf{r})=2c_{0}e\sinh \left[e\beta {\varphi }_{\mathrm{in}}(\mathbf{r})\right]& \forall \mathbf{r}\in \Omega \\ {\epsilon }_{\mathrm{out}}{\nabla }^2{\varphi }_{\mathrm{out}}(\mathbf{r})=0& \forall \mathbf{r}\notin \Omega .\end{array}$$

$${\left[{\varphi }_{\mathrm{out}}(\mathbf{r})-{\varphi }_{\mathrm{in}}(\mathbf{r})\right]}_{\forall \mathbf{r}\in \partial \Omega }$$

$$\mathbf{n}(\mathbf{r})\cdot {\left[{\epsilon }_{\mathrm{in}}\nabla {\varphi }_{\mathrm{in}}-{\epsilon }_{\mathrm{out}}\nabla {\varphi }_{\mathrm{out}}(\mathbf{r})\right]}_{\forall 0<r<R,0<z<L}$$

$${{\varphi }_{\mathrm{out}}(\mathbf{r})|}_{r\to \infty }$$

$$F[\{c_i(\mathbf{r})\}]$$

$$\rho (\mathbf{r})$$

$$F[\{c_i\}]$$

$$\gamma$$

$${\rho }_{\mathrm{ext}}^0$$

$$\frac{\beta F[\{c_i\}]}{V}$$

$$\displaystyle+c_{i}\left(\log(\lambda_{T}^{3}c_{i})-1\right)\Bigg{]},$$

$$\beta {\mu }_i$$

$$\log \left({\lambda }_T^3{c}_i^0\right)$$

$$c_i$$

$${\rho }^0$$

$$\epsilon {\nabla }^2\varphi (\mathbf{r})=-\sum _{i=1}^{M}e{z}_i{c}_i^0\exp \left[-z_{i}e\beta \overline{\varphi (\mathbf{r})}\right].$$

$${\rho }^0=-\frac{{\rho }_{\mathrm{ext}}^0}{1+(e\gamma {\sum }_{i=1}^M{z}_i^2{c}_i^0{)}^{-1}}.$$

$${\int }_V\frac{{\mathrm{d}}^3\mathbf{r}}{V}{\int }_V\frac{{\mathrm{d}}^3{\mathbf{r}}^{\prime }}{{\lambda }_D^2|\mathbf{r}-{\mathbf{r}}^{\prime }|}\gg 1$$

$${\gamma }^{\mathrm{apprx}}$$

$${\rho }_{\mathrm{ext}}^{0,\mathrm{apprx}}$$

$$L\gg (2R)\exp \left(\frac{2{\lambda }_D^2}{R^2}\right),$$

$${\rho }^0$$

$$\tilde{\rho }$$

$${\tilde{\rho }}_{\mathrm{ext}}$$

$$\xi$$

$$\tilde{\rho }=f_{\xi }^{-1}({\tilde{\rho }}_{\mathrm{ext}})-{\tilde{\rho }}_{\mathrm{ext}}=\{\begin{array}{ll}-\frac{{\tilde{\rho }}_{\mathrm{ext}}}{1+{\xi }^{-1}},& {\tilde{\rho }}_{\mathrm{ext}}\ll 1.\\ \log \left(\frac{{\tilde{\rho }}_{\mathrm{ext}}}{\xi }\right)-{\tilde{\rho }}_{\mathrm{ext}}.& {\tilde{\rho }}_{\mathrm{ext}}\gg 1.\end{array}$$

$${\lambda }_D=R\sqrt{2\log \frac{L}{2R}}.$$

$${\lambda }_D=\sqrt{{\lambda }_{\mathrm{GC}}R}{\left(\frac{L}{2R}\right)}^{R/2{\lambda }_{\mathrm{GC}}},$$

$$c_{\mathrm{tot}}=\sqrt{(2c_0{)}^2+({\rho }^0/e{)}^2}.$$

$${\rho }^0=\sum _{i=1}^{N}e{z}_i{c}_i^0\exp \left[-\beta E_i-z_i\gamma ({\rho }^0+{\rho }_{\mathrm{ext}}^0)\right].$$

$${\tilde{\rho }}_{\mathrm{ext}}$$

$$\xi$$

$${\epsilon }_{\mathrm{in}}\frac{\partial {\varphi }_{\mathrm{in}}(r,z)}{\partial r}={\epsilon }_{\mathrm{out}}\frac{\partial {\varphi }_{\mathrm{out}}(r,z)}{\partial r}+q_s$$

$${\partial }_r{\varphi }_{\mathrm{out}}(R,z)\approx -\frac{{\varphi }_{\mathrm{out}}(R,z)}{M_{L/R}R}.$$

$$M_{L/R}=\log \left(\frac{2L}{\pi R}\right)-{\gamma }_{\mathrm{Euler}},$$

Full text

I Introduction

The transport of ions in extreme confinement has applications ranging from physiology to chemical engineeringHille (2001); Doyle et al. (1998); Berneche and Roux (2001); Hess and Tsien (1984); Werner et al. (2001); Vilatela and Marcilla (2015); Yang et al. (2003); Kim et al. (2010a, b); Mani and Bazant (2011); Peters et al. (2016); Schmuck and Bazant (2015). Whether we consider ions traveling through the protein channels in the cell membrane or through pores in an ion-exchange membrane, the underlying physics shares many similaritiesNonner and Eisenberg (1998); Tedesco et al. (2016). A growing interest in ionic transport through nanopores has emerged in recent years owing to nano-fabrication advances that enable us to study pristine nano-channels at the single-channel level, such as carbon nanotubes, boron-nitrite nanotubes or silicon nano-channelsLee et al. (2010); Agrawal et al. (2017); Siria et al. (2013); Guan and Reed (2012); Daiguji et al. (2004). These experiments have revealed that our understanding of even the basic physics is incomplete, and important theoretical knowledge gaps still existFaucher et al. (2019).

Classical theories of ion transport were inspired by the membrane technology that was available at that timeStaverman (1952); Gross and Osterle (1968); Fair and Osterle (1971). The complicated interplay of fluxes and potential gradients (chemical, electrical, and pressure) was naturally modeled with continuum theories. The so-called ``capillary-pore model," based on the Poisson-Boltzmann (PB) equation for the charge distribution normal to the pore walls under the local equilibrium assumption, with Navier-Stokes and Nernst-Planck equations for the fluid flow and ionic flux, is a continuum linear response theory of transport in charged cylindersPeters et al. (2016), which is widely used in different electrochemical applications, from electro-osmotic pumps to energy conversion devicesSasidhar and Ruckenstein (1981); Yao and Santiago (2003); Petsev and Lopez (2006); Kim et al. (2010c); Chang and Yang (2011); Yan et al. (2013). A competing school of thought, rooted in transition state theoryZwolinski et al. (1949), emerged in the biophysical community. Experiments on ion channels show that when open, the transport of the ions is best described by a discrete, single-file reaction model, where ions are attached to specific binding in the channel by a chemical reactionHille (2001); Hess and Tsien (1984); Almers and McCleskey (1984); Doyle et al. (1998); Berneche and Roux (2001); Sather and McCleskey (2003).

Evidently, a new theory is required for nanochannels with a pore diameter of less than $10$nm (single-digit nanopores, or SDNsFaucher et al. (2019)), in order to span between the two limiting regimes, from discrete to continuum behavior. While the two traditional pictures have some merit, neither can exactly capture experimentally measured conductivity curves. Unusual scaling behavior of ionic conductance in carbon nanotubes (CNTs), for example, was recently reported by Secchi et alSecchi et al. (2016) and was subsequently interpreted with a Space-Charge continuum modelBiesheuvel and Bazant (2016). The conductance of a narrow CNT porin, in contrast, was fitted to a Michalis-Menten reaction modelTunuguntla et al. (2017), suitable for a single-file transport mechanism.

In this work, we propose a new theoretical framework for electrolytes in nanopores, consisting of confined ions with three-dimensional (3d) electric fields that extend into the surrounding matrix. In the most interesting and relevant case of SDNs, we construct a truly one-dimensional (1d) mean-field theory of ions confined to a long, thin nanopore in a 3d matrix, in contrast to previous models of ion chains with 1d Coulomb interactions [33-35]. Surprisingly, $1$d electrolytes exhibit several interesting behaviors, most notable is the breakdown of global electroneutrality: a system can have a net charge where the total charge of the ionic solution does not exactly cancel out external charges. When the pore diameter is comparable to the spacing between ions, the system essentially behaves like a $1$d correlated electrolyte. Ions are not necessarily restricted to transport in single-file, but the nature of the electrostatic interactions resembles a $1$d chain. Despite our interest in transport, this paper will only focus on the equilibrium properties of ion channels. The equilibrium properties can in turn be used to understand transport properties of the nanotubeBiesheuvel and Bazant (2016); Gross and Osterle (1968).

Electroneutrality breakdown in nanopores has been observed in Monte-Carlo simulations Sørensen and Sloth (1992); Rivera and Sørensen (1994); Lo and Chan (1995); Lee and Chan (1997); Lee et al. (1999); Boda et al. (2000, 2002, 2006), and was recently even measured experimentallyLuo et al. (2015). However, it was not interpreted as a unique feature of the $1$d geometry. Instead, the breakdown was assumed to occur due to an excess screening of charges outside of the pore. This type of local breakdown of charge neutrality is not suited to most transport problems, where the channel is surrounded by a constant dielectric medium.

Without charge neutrality, electric fields leak out of the confined region into the outer substrate. We derive (section II) a mean-field theory of a confined electrolyte by properly accounting for the outer region as well, and present and illustrate using numerical simulations the emergence of electroneutrality breakdown. In Section III we consider a uniformly charged pore and solve a self-consistent algebraic mean-field equation for the excess charge. Three length-scales govern the accumulated charge: the Debye screening length (ion-ion interactions), the Gouy-Chapman length (ion-wall interactions) and the pore diameter. In Section IV we present a general scaling argument for the enhanced screening length in low dimensions. We show that electroneutrality breakdown is a unique feature of $1$d systems due to their exponentially long screening length that can easily exceed the size of the system. Finally (section V), we solve a full $1$d lattice mean-field equation and observe the emergence of ion-ion correlations at high concentrations. The breakdown of electroneutrality has profound implications on the transport of ions through nanochannels, and in section VI we show that our model can account for the unusual scaling of conductance in CNT.

II Mean-field theory of confined electrolytes

Electroneutrality is often a hidden assumption of continuum models: it hides in the boundary conditions for the Poisson-Boltzmann equation, where the electric fields outside the electrolyte are assumed to vanish Biesheuvel and Bazant (2016); Peters et al. (2016); Gross and Osterle (1968); Fair and Osterle (1971); Sasidhar and Ruckenstein (1981); Secchi et al. (2016); yao2019strong; ren2008slip; cervera2005poisson; van2006electrokinetic; vlassiouk2008ionic; stein2004surface; li2015direct. According to Gauss's law, if there is no electric flux emanating from the electrolyte in the pore, it has zero net charge. Electroneutrality relies on the nanometer-scale screening length, which guarantees macroscopic charge fluctuations are negligible. A rigorous analysis requires us to solve the Laplace equation outside the electrolyte, in addition to the PB equation inside.

It is important to note that electroneutrality is not always assumedLozada-Cassou et al. (1996a, b); Luo et al. (2015); Colla et al. (2016); Schmuck and Bazant (2015). As a recent example, Colla et alColla et al. (2016) considered two charged plates immersed in water, with free ions on both sides of each plate, and solved a density functional theory (DFT) in the entire space. Since the screening is not symmetric, especially when the screening length is large, the accumulated charge between the plates can be small. The authors consider this as an example of a local breakdown of electro-neutrality (LEC), and while it shares some similarities with our approach, LEC is fundamentally different from the pore-wide electroneutrality breakdown which we discuss, and is not a unique property of a $1$d geometry.

II.1 General equations

Let us consider the PB equation for a symmetric binary monovalent electrolyte, fixed at a chemical potential that is set by an external reservoir with ionic concentration $c_{0}$, and confined to a small region in space ($\Omega$, see Fig. 1). We further assume that the electrolyte is embedded in a constant dielectric medium with permittivity $\varepsilon_{\rm out}$, and the boundary is charged with a surface charge $q_{s}$. The electrostatic potential ($\phi$) is determined by a set of PB and Laplace equations:

[TABLE]

where $e$ is the electron charge, $\varepsilon_{\rm in}$ and $\varepsilon_{\rm out}$ are the dielectric constants in the solvent and dielectric matrix, respectively, and $\beta=1/k_{\mathrm{B}}T$ is the inverse temperature and $k_{\rm B}$ is the Boltzmann constant. A similar approach was previously introduced to calculate the effect of image charges on the ionic self energy in confinementJordan et al. (1989), and to study the transport of ions through porous mediaSchmuck and Bazant (2015), but the resulting electroneutrality breakdown was not emphasized.

For a cylindrical geometry, with charge density $q_{s}$ on the pore walls, the boundary conditions for this system are:

[TABLE]

where $\partial\Omega$ is the electrolyte boundary, and ${\bf n}$ is an outward unit vector normal to the boundary.

Solving the set of PB/Laplace equations for a finite cylinder can only be done numerically. In many cases, we will find that accounting for outer electric fields is redundant, since the outer electric fields vanish. In other cases, dramatic differences in the charge profile can be observed. In the the remaining sections of the paper, we will present simplified models that permit analytical results, which will help us explore the implication of Eq. 1. While some of the accuracy of the complete model is lost, we will gain a much better physical understanding, as well as mathematically convenient approximations.

II.2 Numerical simulation

Remarkably, despite numerous theoretical investigations of ionic channels, there is currently no analysis of this simple and fundamental problem, to the best of our knowledge. In many cases, the full set of PB/Laplace equations is oversimplified by disregarding the outer dielectric matrixBiesheuvel and Bazant (2016); Peters et al. (2016); Gross and Osterle (1968); Fair and Osterle (1971); Sasidhar and Ruckenstein (1981).

We solve equation 1 using the COMSOL Multiphysics simulation package. Fig. 2 shows an example of such numerical simulation for a nanopore configuration and the resulting electrostatic potential in a cross section of the cylinder. The simulation box is composed of a membrane dielectric embedding a charged pore that is between two reservoirs. The reservoirs have a thickness of $20\lambda_{D}$ and each side of the box is $5L$. At the top and the bottom of the box, Dirichlet boundary conditions are imposed, to signify a return to bulk reservoir conditions with $\phi=0$ and $c_{i}=c_{0}$. At the lateral surfaces, symmetry conditions are imposed. For simplicity, we analyze the limit of small surface charge density, $q_{s}$, on the pore walls. We also consider the case of $\varepsilon_{in}=\varepsilon_{out}=80\varepsilon_{0}$. Note that there is a significant `leakage' of the electric field into the dielectric matrix, demonstrating the invalidity of the local electroneutrality assumption, as illustrated in Fig. 2b. In Fig. 2c, we quantify the extent of electroneutrality breakdown by calculating the ratio of ionic charges within the cylinder, $Q_{in}$, to the number of charges on the pore walls, $Q_{ext}$. For $\mid Q_{in}/Q_{ext}\mid\rightarrow 1$, the pore walls and the pore are electroneutral, while for $\mid Q_{in}/Q_{ext}\mid\rightarrow 0$, the pore walls and the pore domain experience electroneutrality breakdown.

Note that the continuum simulation as a whole (pore walls, pore, and reservoirs) satisfies global electroneutrality. In fact, one way to check that the reservoir domains are sufficiently large is to ensure that the overall system satisfies the global electroneutrality constraint. When electroneutrality breakdown occurs within the pore domain, the excess screening charge exists at the membrane-reservoir interfaces. Even so, the breakdown in electroneutrality within the pore domain remains an unexpected phenomenon, especially as the pore length goes to macroscopic scales many times longer than the bulk Debye length. Through a series of analytical approximations, we will show that electroneutrality breakdown arises especially strongly in 1d confinement, where the screening length becomes exponentially long.

III Uniform Embedded Pore model

The numerical solution of Eq. 1 in a long and narrow cylinder showed that overlapping double layers are accompanied by a net charge of the pore. With this result in mind, let us now introduce an analytical model of a uniformly charged pore embedded in a constant dielectric medium. This Uniform Embedded Pore" (UEP) model is closely related to the widely used Uniform Potential" (UP) modelSchlögl (1955); Sonin (1976); Cwirko and Carbonell (1989); Verbrugge and Hill (1990); Bowen and Welfoot (2002); Kim et al. (2010b); Peters et al. (2016), also known as Teorell-Meyers-Sievers theory Cowan et al. (1937); Teorell (1953); Meyer and Sievers (1936), which takes advantage of the narrow pore geometry to approximate a constant charge distribution in the radial direction. This approximation is further adapted in the ``Leaky membrane model", which describes the concentration polarization in porous media, based on similar microscopic assumptionsMani and Bazant (2011); Yaroshchuk (2012); Dydek and Bazant (2013). While we make the same assumption about the charge distribution within the pore, we recognize that the sourrounding dielectric matrix cannot be neglected. At first, to keep matters simple, we only consider the case where the dielectric constants inside and outside of the pore are equal.

III.1 Derivation

The free energy functional of $M$ ionic species immersed in a dielectric continuum with permittivity $\varepsilon$, assuming ideal mixing entropy, with thermal de Broglie wavelength $\lambda_{T}$, reads:

[TABLE]

where $c_{i}(\mathbf{r})$ and $z_{i}$ are the concentration and valency of the $i$-th ionic species, respectively, and $\rho_{\rm ext}$ is an external charge distribution. For a uniform ionic density we write the free energy as a function of the mean ionic concentrations:

[TABLE]

The ``mean-interaction" integrals in the above equation describe the electrostatic energy associated with uniformly distributed charge inside a volume, $V$. They depend only on the external charge and geometry, so it is convenient to define the following integrals:

[TABLE]

The free energy density is now simplified:

[TABLE]

where we neglect constant contributions to the free energy. The chemical potential is the derivative of the free energy density with respect to concentration:

[TABLE]

The first term in the LHS of Eq. 7 describes the excess free energy due to electrostatic interactions, while the second term is the ideal gas entropy. If our system is in chemical equilibrium with a bulk reservoir, the chemical potential will only have the second term ($\beta\mu_{i}^{0}=\log\left(\lambda_{T}^{3}c_{i}^{0}\right)$), with the bulk values of ionic concentrations. Equating the chemical potential in Eq. 7 to the bulk value, we find the following equation for the average ionic density, $c_{i}$:

[TABLE]

The average charge density, $\rho^{0}=\sum ez_{i}c_{i}$, therefore, satisfies the following equation:

[TABLE]

We can incorporate non-idealities to the system by adding an excess chemical potential to Eq. 7. If this were the case, the $c_{i}^{0}$ would be replaced by the ionic activities, rather than the concentrations.

The equation for the mean charge in the UEP model, Eq. III.1, is written in terms of the average charge density, $\rho^{0}$, which makes it easier to solve. To see its relation with the set of PB/Laplace equations discussed in section II, we write Eq. III.1 in terms of the electrostatic potential. Note that the term in the exponent is the electrostatic potential, averaged over the volume of the electrolyte, and multiplied by $z_{i}e\beta$. Furthermore, the Poisson equation relates charge density and electrostatic potential, $\rho^{0}=-\varepsilon\nabla^{2}\phi$, and thus Eq. III.1 becomes a partial differential equation:

[TABLE]

Hence, the uniform embedded pore model is an approximation of the standard PB equation, where we replace the electrostatic potential with its volume average. The complete set of equations for the potential in the UEP model includes in addition the Laplace equation in the outer dielectric medium, and the boundary conditions of Eq. II.1. This set of equations is mathematically equivalent to Eq.III.1, as long as the permittivity is everywhere the same.

Under electroneutrality the charge density, $\rho^{0}$, exactly cancels out the external charge distribution, $\rho_{\rm ext}^{0}$. A criterion for electroneutrality to be satisfied, for weakly charged systems, is obtained by linearzing equation III.1:

[TABLE]

Hence, the condition for electroneutrality is $e\gamma\sum_{i=1}^{M}z_{i}^{2}c_{i}^{0}\gg 1$. Using the definition of $\gamma$ (Eq. III.1), we can write this condition in terms of the bulk Debye length- $\lambda_{D}^{-2}=4\pi l_{{\rm B}}(\sum_{i=1}^{M}z_{i}^{2}c_{i}^{0})$:

[TABLE]

In typical $2$ and $3$ dimensional systems, this condition is satisfied if the characteristic size of the system is larger than the Debye length, but in $1$d we find a very different result.

III.2 Cylindrical nanopore geometry: a phase diagram

Our numerical solution suggests that electroneutrality breaks down for cylindrical nanopores if the Debye length is greater than the pore radius. We can now see that this is a property of the cylindrical mean-interaction integral. Let us consider a cylinder with radius $R$, length $L$, and a surface charge density $q_{s}$ ($\rho_{\rm ext}=q_{s}\delta(r-R)$). The mean-interaction integrals, $\gamma$, and $\rho_{\rm ext}$, can be approximated to a good precision by considering a test charge located at the center of the pore:

[TABLE]

Interestingly, this approximation deviates by only a few percents from an accurate numerical evaluation of the mean-interaction integrals. Based on this approximation, the electroneutrality condition ($e\gamma^{\rm apprx}\sum z_{i}^{2}c_{i}^{0}\gg 1$) in nanopores can be expressed in terms of the natural system length scales:

[TABLE]

which will be derived via scaling arguments in section IV. To extend the electroneutrality condition beyond weakly charged systems, the full solution of Eq. III.1 is required. For a monovalent binary electrolyte, Eq. III.1 reads:

[TABLE]

We note that $\gamma^{\rm apprx}$ has dimensions of inverse charge density, so it is instructive to study the dimensionless charge density, $\tilde{\rho}=\gamma^{\rm apprx}\rho^{0}$. According to Eq. III.2, a dimensionless charge density with a value of $1$ describes a system where the distance between ions is proportional to the Bjerrum length. Hence, we expect strong ion-ion interaction when the dimensionless charge density is large, while for weak charge densities ($\tilde{\rho}\ll 1$), we expect weak electrostatic interactions, that would result in the breakdown of electroneutrality.

In terms of the dimensionless charge density we obtain the following algebraic equation:

[TABLE]

where the two parameters, $\xi$ and $\tilde{\rho}_{\rm ext}$, are related to two important length-scales of the system, the Debye-length and the Gouy-Chapman length:

[TABLE]

The Gouy-Chapman length ($\lambda_{\rm GC}=e/2\pi q_{s}l_{{\rm B}}$) characterizes the distance from a charged surface at which its electrostatic and thermal energies are equal. If this distance is much smaller than the pore radius, the system is effectively three-dimensional, and electroneutrality holds. As pore charge decreases, $\lambda_{\rm GC}$ becomes comparable with the pore radius and the one-dimensional geometry is recovered. Our parameter space is therefore described by three length-scales: the Debye screening length, Gouy-Chapman length, and the pore diameter. The breakdown of electroneutrality into a disordered phase is further augmented by the aspect ratio of the system, which effectively re-scales the pore diameter.

The solution to Eq. 16 can be written as $\tilde{\rho}=-\tilde{\rho}_{\rm ext}+f_{\xi}^{-1}(\tilde{\rho}_{\rm ext})$, where $f_{\xi}(x)=x+\xi\sinh(x)$. Approximated solutions can be found for highly and weakly charged pores:

[TABLE]

This solution is illustrated in Fig. 3 as a function of $\lambda_{GC}$ and $\lambda_{D}$, for a pore with dimensions $R=1$nm and $L=100$nm. We identify four different regimes, as shown in the four panels of Fig. 3. At low surface charge (large $\lambda_{GC}$) and short screening lengths (thin double layer), the system can be described by the standard DH approximation, and the electroneutrality assumption is valid. This theory also covers the beginning thick double layer regime, where the double layers begin to overlap. For high surface charges the linearized Poisson-Boltzmann equation is no longer valid. In this strong coupling regime the full non-linear equation is required, but interestingly, it also ensures the ions will completely screen any surface charge. Our solution, however, predicts a fourth regime, of low surface charge and small concentration. Under these conditions, electroneutrality is broken, and external fields must be accounted for.

Two theoretical curves mark the boundaries of the electroneutral phase. The vertical line is derived from the weak charge approximation, as the Debye length at half screening:

[TABLE]

The second curve mark the transition to electruneutrality due to high surface charge, and is obtained by requaring half screening in the strong coupling limit of Eq. III.2:

[TABLE]

which asymptotes to a horizontal line at $\lambda_{GC}\approx R$ as the concentration decreases. Note that the transition to electroneutrality is slow, and spans roughly an order of magnitude change in the parameters. Electroneutrality breakdown can thus play a major role in the physics of nanopores.

The total ionic concentration inside the nanopore, depicted in Fig. 4, is related to the pore charge by Donnan equilibrium:

[TABLE]

In electroneutral systems the concentration reaches a plateau in the dilute limit, where the only ions inside the channel balance the surface charge. Accounting for electroneutrality breakdown, however, significantly alters the behavior in the dilute limit.

III.3 Sub-nanometer nanopores: dehydration and images forces

Ion specific effects have important consequences on the behavior of nanopores, especially in the sub-nanometer scale. In extremely narrow pores ions have to remove their hydration shell, which creates a large energy barrier for entering the pore Epsztein et al. (2019); Sahu and Zwolak (2017). Denoting this energetic cost for the $i$th specie by $E_{i}$, we generalize the self-consistent equation for the mean charge distribution:

[TABLE]

In the monovalent case the pore-charge equation takes a similar form to Eq. (16), with re-scaled coupling parameter $\xi$ and external charge $\tilde{\rho}_{\rm ext}$:

[TABLE]

Hence, the phase diagram (Fig. 3) remains similar, but skewed: the x-axis is rescaled by the average Boltzmann weight, while the y-axis is shifted by dimensionless energy difference. Any asymmetry in the dehydration energy will result in an excess charge within the pore since the dehydration energy plays a similar role to that of the surface charge.

Even though we incorporate in our model an energy barrier, it is important to emphasize again that we are considering a model with a constant permittivity everywhere. The differences in self-energies are only one aspect of a dielectric mismatch, that can change ion-ion and ion-wall interactions as wellKamenev et al. (2006); Zhang et al. (2005, 2006); Levin (2006); Bordin et al. (2012).

III.4 Towards a general theory: the electric leakage boundary conditions and validation the UEP model

We examined in this section a quantitative model for electroneutrality breakdown, and obtained a simple algebraic equation that determines the net charge inside a nanopore. In agreement with the numerical analysis of section II, we find that breakdown occurs only in the thick doubly layer limit. Yet, our model does not allow us to resolve the spatial distribution of ions inside the channel. As we show in section II, the full charge density can be found by solving the combination of Laplace and PB equations (Eq. 1) with the BC of Eq. II.1.

For narrow and long pores ($L\gg R$) it is tempting to simply take $L\rightarrow\infty$ and solve a $1$d ordinary differential equation. However, as we have shown here, even in very long channels the lengths plays an important role. Hence, solving a cumbersome set of two coupled nonlinear $2d$ PDEs to get the full charge density is required. For weakly charged surfaces (DH regime) an exact analytical solution is available in the form of an infinite series (see Appendix A), but it is difficult to gain any insight from this complicated expression.

It turns out, however, that for long and narrow channels we can go one step beyond the UEP model, and obtain approximated boundary conditions for PB equation. With this approximated BC, we can again solve a simple $1$d ODE, but a one that inherently captures electroneutrality breakdown.

Recall that the in cylindrical coordinates boundary conditions are:

[TABLE]

In Appendix B we show that the electric field and electrostatic potential in the outer region are related by a simple expression:

[TABLE]

where $M_{L/R}$ is defined as:

[TABLE]

and $\gamma_{\rm Euler}\approx 0.577$ is Euler's constant. Using the electrostatic potential continuity condition, the new Robin type ``electric leakage" boundary condition is:

[TABLE]

In Fig. 5 we compare a numerical simulation to the exact analytical solution (Eq. A) for the total charge in a weakly charged nanopore, as well as to the UEP model and the new BC. All approximated models are compared with the numerical COMSOL simulations as well. To satisfy the weak potential assumption ($e\beta\phi\ll 1$) of the Debye-Hückel approximation, we focused only on the limit where the charge density approaches [math]. As expected, in the dilute limit the UEP model, as well as the approximated boundary condition, perfectly match the exact analytical result. As concentration increases, the uniform approximation assumption is no longer valid.

The excellent match of the tractable UEP model to more elaborate schemes allows us to use it when we are interested in fitting experimental measurements. In section VII we will explore some available data, and show the importance of electroneutrality breakdown to transport properties. First, however, we would like to take a small detour and discuss why electroneutrality was obtained. As we mentioned in the introduction, we claim that this is a unique feature of the $1$d geometry. In the next two sections we discuss how the basic physics of $3$d electrostatic interactions between particles confined to a $1$d narrow channel, and show that the UEP model is best understood as an example of $1$d theory. The impatient reader, however, can skip directly to section VI.

IV Screening length for ions confined to lower dimensions

A central ion in an aqueous solution is surrounded by a screening cloud meaning that ions in its vicinity tend to be of opposite charge. In three dimensions, we can find an analytical expression for the shape of the screening cloud, achieved by solving the linearized Poisson-Boltzmann equation. The density profile of screening ions decays exponentially with a characteristic length, the ``Debye length" ($\lambda_{{\rm D}}$). When ions are restricted to different dimensionality, still interacting with a $1/r$ pair-wise potential, there is no similarly tractable equation to determine the shape and size of the screening cloud. We will show in the following sections, that an equivalent screening length exists and plays an important role in the physics of lower dimension electrolytes, but it formally requires a cumbersome derivation. Before we delve into a more rigorous formalism, we first present a simple scaling argument for the screening length that holds in any dimension.

The screening length, in essence, is the distance at which entropic and electrostatic forces balance each other. The electrostatic force pulls the screening cloud closer to the central ion, while entropy favors uniform charge distributions and pushes the screening cloud away. If we consider a screening cloud of radius $\lambda_{s}$ with $N_{s}$ ions around a negatively charged univalent ion, the probability of each ion to have charge $+e$ equals $p^{+}_{i}=1/2+1/2N_{s}$. Maximal entropy is achieved if there is an equal probability for an ion to be positive or negative. By forcing the ions to screen the central charge, the probability of positively charged ions slightly increases, and the entropy is reduced. Note that the ions outside the screening cloud have equal probabilities to have positive/negative charge. The entropy associated with this screening cloud equals (assuming $N_{s}\rightarrow\infty$, which corresponds to the dilute limit):

[TABLE]

where $N$ is the total number of ions in the system, both inside and outside the screening cloud. The first term is constant (does not depend on $N_{s}$) and can be ignored. Note that this is only the entropy associated with the possibility of each ion to be positive or negative, and positional entropy is neglected. Relating the entropy to the $d$-dimensional sphere volumeMcLaren-Young-Sommerville (1958), $V=A_{d}\lambda_{s}^{d}/d$, where $A_{d}$ is the surface area of a unit sphere, and average ionic pair concentration ($2c=N_{s}/V$), we find that the entropy of ions inside the sphere is reduced by:

[TABLE]

The electrostatic energy of a uniformly distributed screening cloud equals:

[TABLE]

where $a$ is the ion size, assumed to be much smaller than $\lambda_{s}$. Here, we make the critical assumption of a $3$d Coulomb potential decaying as $1/r$, which effectively spills out of the confining geometry, and neglect for now any dielectric response of the surrounding matrix, which modifies the result but does not alter the basic scaling arguments (as shown below).

Minimizing the free energy ($F=U-T\Delta S$) with respect to $\lambda_{s}$ , and ignoring the numerical prefactors, yields the following scaling behavior for the screening length:

[TABLE]

where $l_{{\rm B}}=e^{2}/4\pi k_{\mathrm{B}}T\varepsilon$ is the Bjerrum length. In three dimensions, we recover the standard Debye screening: $\lambda_{s}\propto\lambda_{D}=(8\pi cl_{B})^{-1/2}$. In a nanoslit geometry of two dimensions, the screening length is proportional to $(cl_{B})^{-1}$. Forcing ions to reside on a two-dimensional plane reduces electrostatic interactions, and slightly increases the size of the screening clouds. A much more dramatic effect is observed for one-dimensional confinement in a long, thin nanopore, where the screening length is exponentially large, as shown in Fig. 6.

If we consider our one-dimensional system to be a cylinder with radius $R$, length $L$, and a $3$d concentration (per volume) of $c_{0}$, the corresponding $1$d concentration (per length) in the axial direction is $c^{\rm 1d}=\pi R^{2}c_{0}$. In terms of these bulk properties, we find that the dimensionless factor, $c^{\rm 1d}l_{B}=\pi R^{2}c_{0}l_{B}$, is related to the ratio of Debye length to pore radius, which enters the exponential scaling of the screening length in one dimension: $\lambda_{s}\propto a\exp\left[(\lambda_{{\rm D}}/R)^{2}\right]$. For narrow pores in dilute electrolytes, the Debye length can easily be greater than the pore radius, which is the traditional limit of ``thick double layers" ($\lambda_{D}\gg R$), but in contrast to classical continuum models, we predict that this results in an extremely large screening length, easily exceeding all geometrical scales in the problem. In particular, when the screening length exceeds the total length of the channel($\lambda_{s}\gg L$), the central ion is no longer screened, and electroneutrality breaks down.

V Mean-field theory on a one-dimensional lattice

In the next section we compare the UEP model with experimental data, and show it can be used to interpret conductance curves. Before doing so, we would first like to show quantitatively how electroneutrality breakdown is a fundamental property of $1$d confinement, by considering a system of ions restricted to reside along a line. This will also allow us to explore ion-ion correlations along the pore axis, and observe a transition from order to disorder.

We study a lattice-gas model, and not a continuum model, for two reasons. First, it will enable us to discuss packing constraints at the high concentration limit. A more fundamental reason was hinted in the previous section: there is no equivalent continuum PB model in $1$d. The scaling argument showed that the screening length depends on the minimal distance between ions ($a$), and this will remain valid in the analysis here as well. As a result, we cannot find a corresponding differential equation that describes the system in the continuum limit. Note that in contrast to many previous $1$d models of electrolytes (Edwards and Lenard (1962); Lenard (1961); Démery et al. (2012)), the electrostatic interaction is three dimensional ($1/r$): we study point-like ions along a line, and not parallel charged sheets.

For mathematical convenience, we consider the free energy functional of a 1d periodic lattice model (a ring), with lattice spacing equal to the ionic size, $a$ (see Fig. 7). The $i^{\rm th}$ site can be occupied by a positive ion, negative ion or a vacancy, with probabilities $p_{i}^{+}$, $p_{i}^{-}$ and $1-p_{i}^{-}-p_{i}^{+}$, respectively. The electrostatic energy, $U$, is given by:

[TABLE]

where $\phi^{\rm ext}_{i}$ is an external electrostatic potential. We denote the dimensionless average charge vector $q_{i}=(p^{+}_{i}-p^{-}_{i})$ and the dimensionless interaction matrix $\tilde{\Phi}_{ij}=|i-j|^{-1}$, so the electrostatic energy reads:

[TABLE]

To calculate the free energy, we add the entropy of mixing:

[TABLE]

If connected to a particle reservoir, the chemical potential is set externally and is calculated by the functional derivative of the free energy density $f=(U-TS)/a$ with respect to concentration ($c^{\pm}_{i}=p^{\pm}_{i}/a$):

[TABLE]

Rearranging the terms, we obtain the $1$d analog of the Bikerman model (Bikerman (1942); Borukhov et al. (1997); Kornyshev (2007); Kilic et al. (2007); Bazant et al. (2009)):

[TABLE]

where the fugacity $\Lambda={\rm e}^{\beta\mu}$ is proportional to the bulk reservoir ion activity.

The general form of the $1$d ring equation is not very illuminating, and specific examples are required to show how ion-ion correlations play an essential role in this model. We study two systems: a charged homogeneous nanopore and the charge distribution around a central ion. In the dilute limit, we recover the behavior described in the previous section. As ionic concentration increases, the model naturally predicts a transition to an ordered structure, including the short-range over-screening phenomenon in intermediate concentrations.

V.1 A Uniformly Charged 1-d Nanopore

Let us assume we have a homogeneous charge density in our system: $q_{i}=q$, $\phi^{\rm ext}_{i}=\phi^{\rm ext}$. For periodic boundary conditions, a uniform charge distribution is an eigenvector of the interaction matrix $\tilde{\Phi}$, where the eigenvalue is twice the harmonic number $H_{L/a}$. For a long chain ($L\gg a$) the harmonic number can be expanded:

[TABLE]

If the external potential is due to a uniform charge distribution on the pore walls ($Q$ per site), Eq. 36 becomes the following algebraic equation:

[TABLE]

In the dilute limit, this equation has a similar form to Eq. 16. As the external charge increases, a maximal charge density of one charge per lattice site prevents an unphysical accumulation of charges inside the pore.

Limiting the $1$d-ring equation to only uniform distributions explicitly neglects any ion-ion correlations. When ion-ion interactions are strong enough, the system will form a crystal structure. In contrast to the $3$-dimensional case, where the PB equation cannot predict ordered structures, our $1$d model can be easily extended to include the expected phase transition. To account for patterns of alternating signs, we generalize our argument and replace the constant charge density with the form: $q_{i}=q_{0}+q_{1}(-1)^{i}$. Both uniform charge density and an alternating pattern are eigenvectors of the interaction matrix $\tilde{\Phi}$, where the eigenvalue of $(-1)^{i}$ equals $-2\log 2$. With this functional form we obtain a set of coupled non-linear algebraic equations, for the average charge at even/odd locations:

[TABLE]

Exploring the space of solutions of Eq. V.1 is shown as a phase diagram in Fig. 7, for $a=5$Åand $L=100$Å. As long as the ion-ion correlations are weak, we recover the same behavior found in the continuum model, including electroneutrality breakdown in the dilute limit. However, we find another breakdown of electroneutrality in the high concentration limit. This is the oscillatory regime ($q_{1}>0$). Since $q_{1}+q_{0}<1$, an increase of oscillating term, $q_{1}$, has to come at the expense of the average term, $q_{0}$, and electroneutrality is broken. With an increase of external charge, oscillations are suppressed and electroneutrality is again favored. As shown in Fig. 7, the order-disorder phase boundary only weakly depends on the surface charge and can be evaluated analytically based on the $Q=0$ limit. In this limit the average charge is $q_{0}=0$, and the number of solutions is determined by a single algebraic equation:

[TABLE]

The RHS of Eq. (V.1) is monotonically increasing, starting from [math], and has a monotonically decreasing slope. Hence, a second solution is available only if the slope at $q_{1}=0$ is greater than $1$, which leads to a critical fugacity of:

[TABLE]

The results in Fig. 7 are displayed in terms of the standard $3$d system parameters. With a pore radius $a$, the surface charge density equals $q_{s}=Q/(2\pi a^{2})$. The bulk ionic concentration is related to the fugacity $c^{\rm 3D}=\Lambda/\pi a^{3}$. The Gouy-Chapman length and the Debye length are defined as usual.

V.2 Charge distribution around a central ion

Bulk oscillations are maintained only for concentrations beyond a critical concentration, with persisting long-range order. Slightly below the critical density, we expect temporary short ranged decaying oscillations, that will eventually be replaced by monotonic decaying fluctuations in the dilute limit. We show how this behavior emerges with a standard Debye-Hückel approach, by solving for the charge distribution around a central ion.

We take advantage of the periodic boundary conditions, and constrain the $i=0$ site to have a charge $Qe$ by adding a term $\alpha e^{2}/a\varepsilon(q_{0}-Q)$ to the free energy functional (Eq. 33), where $\alpha e^{2}/a\varepsilon$ is a Lagrange multiplier. The resulting mean-field equation reads:

[TABLE]

which we solve in the linear Debye-Hückel regime:

[TABLE]

where the coupling parameter $\chi=\frac{2\Lambda l_{B}}{a}$ is defined as the ratio of ionic spacing ($a/2\Lambda$) and the Bjerrum length. By the translation symmetry of the matrix $\tilde{\Phi}_{ij}=\tilde{\Phi}_{|i-j|}$, we find a closed-form expression in the discrete Fourier space:

[TABLE]

where the normalization constant $\alpha$ is set such that $q_{0}=Q$. Note that the solution is only valid in the disordered phase, where $\chi^{-1}>2\log 2$, which coincides with the critical fugacity (Eq. 41) in uniformly charged nanopores if size effects are neglected ($a\ll l_{{\rm B}}$). Fig. 8 shows charge density profiles for different concentration, illustrating how Eq. 44 is able to capture both the dilute Coulomb gas limit and the onset of long-range ordering, and predict the transition from screening, to over screening and oscillations.

So far we assumed an infinitely long chain, that allowed us to get a closed-form solution in Fourier space (Eq. 44). As a consequence, the electronutrality is guaranteed, and the total charge along the pore accumulates to [math]:

[TABLE]

However, the decay rate is very slow, and the screening cloud extends many lattice sites. If we look at $\tilde{q}(\omega)$ , the Fourier transform of $q_{i}$, we find a steep increase near $\omega=0$, on its way to its maximal value at $\omega=\pi$. We evaluate its width by finding the frequency at which $\tilde{q}(\omega)$ reaches half of its maximal height:

[TABLE]

Invoking the uncertainty principle, we conclude that the width of the screening cloud scales as ${\rm e}^{\frac{a}{4\Lambda l_{{\rm B}}}}$, in agreement with our scaling analysis (Eq. 31) and the electroneutrality condition (Eq. 12).

Finally, we use our explicit solution of the screening cloud to evaluate the activity coefficient. For this purpose, the electroneutrality breakdown plays only a minor role. Ignoring the finite length of the system, the Debye-Hückel activity coefficient can be written as an integral expression:

[TABLE]

where $\mu^{\rm ex}$ is the excess chemical potential. In dilute systems we can expand the activity coefficient to lowest order in the coupling parameter $\chi$, and get:

[TABLE]

The activity coefficient in confinement is much smaller than the bulk value. As ions cross over to the nanopore, they effectively shed off their ionic screening cloud. For a nanopore of radius $a$, with bulk ionic concentration ($c=\Lambda/\pi a^{3}$), the $1$d activity is only a fraction of the bulk one:

[TABLE]

where we used the standard Debye-Hückel activity coefficient, $\beta\mu^{3D}=l_{B}/2\lambda_{D}$.

As shown in Fig. 8, this approximation is only valid in the dilute limit. As concentration increases, screening and over-screening dominates the electrostatic interactions, and reduce the activity coefficient further.

VI Comparison with experiments

VI.1 Single Digit Nanopores

Our model predicts a non-trivial charge accumulation within the nanopore, due to the breakdown of electroneutrality, which has a direct implication on the measured conductivity of the pore. Assuming equal mobilities for all ions ($\upmu_{D}$), the conductance of the pore is given by:

[TABLE]

where the concentration, $c$, is related to the accumulated charge by Eq. 21. We compare our results with conductance measurements in carbon nanotubes (CNT), taken from Ref. Secchi et al. (2016). It is important to note that the conductance behavior can be explained by different models. The CNT data was originally assumed to have a concentration-dependent surface charge and was later fitted by predicting the adsorption of hydroxyl ions to CNT pore walls (Biesheuvel and Bazant (2016)). Our goal is not to underestimate the importance of a charge regulation mechanism, that can lead to concentration-dependent surface charge by affecting the adsorption rate (Dukhin et al. (2005)), but to suggest a plausible alternative with a constant surface charge. We show that the $1$d geometry by itself can lead to the variety of conductance curves observed in experiments.

Fig. 9 shows the conductance curves for CNTs with varying pore size and surface charge, as a function of KCl concentration. The experimental data were fitted according to Eq. III.3, with three fitting parameters- the ionic mobility, surface charge, and energy barrier.

For the larger pore sizes, shown in Fig. 9(b-d), the energy barrier for entering the pore ($E_{\pm}$) was neglected, and the fitted ionic mobilities were on the order of the bulk mobility of KCl, and varied from $5-10\times 10^{-8}{\rm m}^{2}/{\rm Vs}$, compared to the bulk value of $7.62\times 10^{-8}{\rm m}^{2}/{\rm Vs}$ (see Fig. 11, bulk value of ionic mobility was taken from Lide (2004)). This effectively only leaves one fitting parameter to determine the shape of the curves- the surface charge. As surface charge increases, electroneutrality is maintained especially for the larger pore radii. The curve approaches a constant in the dilute limit (see top curves in figures 9c and 9d). However, for smaller surface charges, the apparent decrease in conductance in the dilute limit is due to electroneutrality breakdown.

In the narrow nanopore (Fig. 9a) the behavior is more complicated. First, a small energy barrier of $\approx 2k_{\mathrm{B}}T$ is needed to obtain the correct trend. This small energy barrier is expected due to a lower dielectric constant inside the nanopores, related to the confinement of water. More importantly, the fitted mobility, in this case, is significantly higher: $12$ times higher than bulk value for the lower two surface charges, and more than a $100$ times greater for the high surface charge (Fig. 11). The large mobility in the high surface charge limit can be due to ion-ion correlations, where the positive charges push each other to move faster. It can also be related to enhanced water flow, due to an increased slip length Holt et al. (2006); Whitby et al. (2008).

VI.2 Sub-nanometer channels

We conclude by applying our model to sub-nanometer channels, which is the relevant limit for ion channels in nature. The most prominent example of a biological nanochannel is the Gramicidin-A channel. With a pore diameter of about $4$Å, it is truly a one-dimensional system. It is often described using Michaelis-Menten type conductance: ions which travel through a channel that connects reservoirs A and B are transitioning between three possible states (A", B" and ``in channel"). This framework is successful since it naturally ignores any charge neutrality constraints. It predicts a linear dependence in dilute solutions and a saturation at high concentrations, limited by the maximal occupancy of the pore. Despite its good agreement with experimental data, it can only describe systems with a handful of ions. Continuum models, on the other hand, which are much better suited to handle numerous ions, cannot predict the dilute limit linear behavior as long as electroneutrality is assumed.

We consider two experimental datasets: a conductance measurement Gramicidin-A channel (taken from Finkelstein and Andersen (1981)), and more recent conductance measurement in a CNT porin experiment (taken from Tunuguntla et al. (2017)). We fit the data according to Eq. 50 as before, with one important modification. As water molecules are excluded from these channels, the dielectric constant is now much smaller and was chosen to be $\varepsilon=5$.

As shown in Fig. 10, our model is able to capture both neutral (linear conductance) and charged (Michalis-Menten conductance) sub-nanometer nanopore behavior. The plateau at high concentrations is not predicted, as our model fails in the concentrated regime. For higher concentrations, a more detailed picture of the (coupled) fluxes has to be accounted for and is beyond the scope of this paper. As shown by the lattice model (see Fig. 7), ion-ion correlations can lead to a decrease in the total charge with increased concentration and eventually to an overall reduction in conductance.

The predicted energy barriers for both experiments were similar ($5k_{\mathrm{B}}T$). This energy is much smaller than the Born solvation energy in the vacuum, which might imply that dehydration is not complete, and is compensated by the interactions of ions with the pore walls. It is also smaller compared with energy barriers estimated by Michalis-Mentan type theories, which are of the order of $10k_{\mathrm{B}}T$Hille (2001). We note, however, that in order to keep the model simple, we assumed a constant energy difference between the pore and its surrounding. Entrance effects were smeared throughout the system, so a smaller energy barrier is expected. A more careful derivation is required to accurately separate the pore mouth contribution to transport.

Another interesting result is that the ionic mobility in GA is much lower compared to bulk water, while the CNT porins have higher mobility. The biological channel has fitted mobility of $5.7\times 10^{-9}{\rm m}^{2}/{\rm Vs}$ which is an order of magnitude less than the bulk ionic mobility (see Fig. 11). The CNTP fitted mobility is closer to bulk KCl and equals $11\times 10^{-9}{\rm m}^{2}/{\rm Vs}$ and $50\times 10^{-9}{\rm m}^{2}/{\rm Vs}$ for the neutral and charged pores, respectively. The reduced mobility compared to isolated ions in a solvent (water) can be understood from the strong attraction of ions and the opposite fixed wall charge and hence larger friction for relative motion. Our expression for conductance is derived based on a normal Nernst-Einstein relation with uncoupled fluxes of each ion. This is only valid for “pseudo binary” transport where each specie only interacts with an abundant solvent. Generally, the Stefan-Maxwell picture has a diffusion tensor with coupled fluxes between each pair of species Donev et al. (2019); Balu and Khair (2018), and can be extended to a ``dusty gas model," if wall molecules are treated as a fixed specie Evans III et al. (1961); Fu et al. (2015).

VII Conclusions

We have analyzed ions in a $1$d electrolyte interacting with $3$d electrostatic interactions, which is the appropriate limit for single-digit nanopores. In $3$d systems the fluctuations around electroneutrality are limited only to the microscopic length-scale. The strong Coloumbic cost of large-scale deviations is much greater than the thermal energy. In contrast, when ions are forced to reside along a line, even macroscopically-long charged chains can spontaneously form.

We first showed directly, using numerical simulations, that a net charge inside long and narrow nanopores is easily achieved when the radius of the pore is small compared to Debye length. We then developed an analytical mean-field model and and predicted a phase diagram for the accumulated charge inside the pore. We found that the pore behavior depends on two length-scales: the Debye screening length and the Gouy-Chapman length. If both length-scales are larger than the pore diameter, the ions do not know about the $3$d nature of the system, and electroneutrality is broken. This behavior is best understood by examining how the competition between electrostatic forces and entropy determines the screening length. In three dimensions we recover the classical short-range Debye screening length, but in $1$d we find an exponentially large screening length.

In our efforts to provide analytical results with a clear physical meaning, we neglected several important aspects of the problem. Most notably, our models fail to account for the polarization charge induced by a discontinuity of the dielectric constant. A large mismatch in the dielectric constant can alter ion-ion interactions inside the poreLevin (2006); Bordin et al. (2012), and can lead, for example, to a $1$d Coulomb interactions in short nanoporesZhang et al. (2005). Molecular dynamics studies have shown that the interactions with images forces are especially relevant for selectivity in ion channels Boda et al. (2006, 2007a, 2007b). Selectivity is also sensitive to the size of the ions, and a proper theory of confined electrolytes must include finite size effects.

Within the uniform embedded pore model, we derived approximated but accurate closed-form expressions for the expected charge and ionic concentration inside a charged pore, and the resulting ionic conductance. With two fitting parameters, the surface charge and ionic mobility, we were able to fit a wide range of conductance curves in narrow nanopores. We interpret the unusual scaling behavior observed as a consequence of the breakdown of electroneutrality.

While transport measurements are an extremely useful tool for studying nanopores, inferring the ionic concentrations is a difficult task. Ion-ion, ion-water and water-pore interactions all play a role in the complicated transport phenomena. A more complete description that explicitly accounts for water flow is required to correctly predict the conductance. Moreover, the mobility of the ions may vary under confinement and compositionEsfandiar et al. (2017); Duan and Majumdar (2010), and the linear relation with concentration is only appropriate in infinite dilution.

Another important aspect of ionic transport, especially relevant to short nanopores, are entrance effects and access resistance. The transition from a microchannel to a nano-channel adds additional resistance to the system, and the access resistance decreases with increasing concentration, which is an alternative explanation for the scaling observed in the conductance, as argued in recent papersGreen et al. (2014, 2016). The entrance effects are of even greater importance if charges are added to the pore mouth, for example, to increase selectivityCorry (2011); Chan et al. (2013); Gravelle et al. (2013).

With current available data, we cannot rule out alternative explanations to describe the scaling and shape of the conductance curves. However, as concentration decreases, our predictions deviate substantially from other models. For example, we predict a linear conductivity in the very dilute limit, and not a plateau. Experiments with a wider range of dilute concentrations are therefore needed to correctly identify the key mechanism.

Added complexities are surely required to adequately describe the transport of ions through nanopores. Yet, as extremely long and narrow nanopores become technologically accessible, determining their net charge is a crucial first step. In our quest to understand the physics of single-digit nanopores, we highlight a simple but consequential observation on the nature of geometrical confinement: it breaks charge neutrality.

Acknowledgements.

This research was supported as part of the Center for Enhanced Nanofluidic Transport (CENT), an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award $\#$ DE-SC0019112 (continuum modeling), and by an Amar G. Bose Research Grant ($1$d lattice-gas model). JPD acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. 1122374.

Appendix A Exact solution of a charged cylinder in the Debye-Hückel regime

In section II.A we found an analytical solution for a weakly charged sphere in the DH regime. We show here the solution for a weakly charged cylinder, in the form of an infinite sum. We consider a cylinder of length $L$, radius $R$ filled with an electrolyte with permittivity $\varepsilon_{\rm in}$ and immersed in a dielectric medium dielectric constant $\varepsilon_{\rm out}$. The set of Debye-Hückel and Laplace equations, derived by linearizing Eq. 1, are:

[TABLE]

The boundary conditions (Eq. II.1) in a cylindrical geometry read:

[TABLE]

Following Barcilon (1992), we also imposed that the potential vanishes at the ends of the cylinder. Taking advantage of the azimuthal symmetry, we can write the inner and outer solutions as an infinite sum:

[TABLE]

where $I_{0}(x)$ and $K_{0}(x)$ are the $0^{\rm th}$ order modified Bessel function of first and second kind, respetively, and $\omega_{n}=\pi n/L$. Note that we wrote the solution in a manner that satisfies the last three boundary conditions. To find $A_{n}$ and $B_{n}$, we plug in the expansion to the first two boundary conditions, where the second boundary condition is expanded in a similar way:

[TABLE]

Clearly, for even $n$ we have $A_{n}=B_{n}=0$. For odd $n$, we obtain the following equations for the expansion coefficients:

[TABLE]

Finally, we get:

[TABLE]

Appendix B Approximated boundary conditions for long and narrow nanopores

In the limit of small aspect ratios ($L\gg R$) we find an approximated Robin type boundary condition for the Poisson-Boltzmann equation, that obviates the need for outer fields solution. We do so by relating the electrostatic potential on the cylinder walls ($\phi_{\rm in}(R,z)=\phi_{\rm out}(R,z)$) to the outer electric field, $\varepsilon_{\rm out}\partial_{r}\phi_{\rm out}(R,z)$. This will allow us to express the jump in the electric field as a function of inner electric field and inner electrostatic potential, instead of solving the outer region as well.

Since the outer region satisfies the Laplace equation, and as long as we have azimuthal symmetry, the electrostatic potential at the surface can be written as the following sum:

[TABLE]

Note that we are not limited to the linearized DH region as in in Appex. A. In the limit $L\gg R$, the modified Bessel function is approximated to lowest order:

[TABLE]

where $\gamma_{\rm Euler}\approx 0.577$ is Euler's constant. We note that to lowest order in $\log(L/R)$ we find:

[TABLE]

which leads to the following approximation for the electrostatic potential:

[TABLE]

where $M_{L/R}$ is defined as:

[TABLE]

Our next step is to relate the derivative of the potential to the potential itself, which will allow us to mask the solution as a boundary condition. Taking the first derivative of the electrostatic potential we obtain:

[TABLE]

Luckily, we find that to lowest order ($K_{1}(x)\approx 1/x$) it is indeed proportional to the electrostatic potential:

[TABLE]

Finally, assuming a surface charge on the pore walls ($q_{s}$), we get a new boundary condition for nanopores:

[TABLE]