# Breakdown of electroneutrality in nanopores

Amir Levy[1], J. Pedro de Souza[2], and Martin Z. Bazant[2,3]  [1] Department
of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts
02139 USA  [2] Department of Chemical Engineering, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139 USA  [3] Department of Mathematics,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139 USA

###### Abstract

Ion transport in extremely narrow nanochannels has gained increasing interest
in recent years due to unique physical properties at the nanoscale and the
technological advances that allow us to study them. It is tempting to approach
this confined regime with the theoretical tools and knowledge developed for
membranes and microfluidic devices, and naively apply continuum models, such
as the Poisson-Nernst-Planck and Navier-Stokes equations. However, it turns
out that some of the most basic principles we take for granted in larger
systems, such as the complete screening of surface charge by counter-ions, can
break down under extreme confinement. We show that in a truly one-dimensional
system of ions interacting with three-dimensional electrostatic interactions,
the screening length is exponentially large, and can easily exceed the
macroscopic length of a nanotube. Without screening, electroneutrality breaks
down within the nanotube, with fundamental consequences for ion transport and
electrokinetic phenomena. In this work, we build a general theoretical
framework for electroneutrality breakdown in nanopores, focusing on the most
interesting case of a one-dimensional nanotube, and show how it provides an
elegant interpretation for the peculiar scaling observed in experimental
measurements of ionic conductance in carbon nanotubes.

##  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:

$\displaystyle\begin{cases}\varepsilon_{\rm in}\nabla^{2}\phi_{\rm
in}(\mathbf{r})=2c_{0}e\sinh\left[e\beta\phi_{\rm
in}(\mathbf{r})\right]&\forall\mathbf{r}\in\Omega\\\ \varepsilon_{\rm
out}\nabla^{2}\phi_{\rm
out}(\mathbf{r})=0&\forall\mathbf{r}\not\in\Omega.\end{cases}$ (1)

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:

$\displaystyle\left[\phi_{\rm out}(\mathbf{r})-\phi_{\rm
in}(\mathbf{r})\right]_{\forall\mathbf{r}\in\partial\Omega}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle{\bf n}(\mathbf{r})\cdot\left[\varepsilon_{\rm
in}\nabla\phi_{\rm in}-\varepsilon_{\rm out}\nabla\phi_{\rm
out}(\mathbf{r})\right]_{\forall 0<r<R,0<z<L}$ $\displaystyle=$ $\displaystyle
q_{s}$ $\displaystyle\left.\phi_{\rm
out}(\mathbf{r})\right|_{r\rightarrow\infty}$ $\displaystyle=$ $\displaystyle
0,$ (2)

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

Figure 1: Sketch of a confined electrolyte, in chemical equilibrium with a
distant external reservoir (in the grand canonical ensemble), as indicated by
an ideal ``salt bridge," which does not otherwise perturb the system. The
governing equation in the inner region ($\Omega$) is the usual Poisson-
Boltzmann equation. The outer region has a fixed dielectric constant and is
described by a Laplace equation. The surface charge on the boundary layer
($\partial\Omega$) determines the jump in the normal component of the electric
field.

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.

Figure 2: Electroneutrality breakdown for a nanopore in equilibrium with an
external reservoir. The surface charge is not fully screened by the confined
electrolyte, as electric fields extend into the solid matrix (assumed here to
have the same dielectric constant as the pore). (a) simulation box in COMSOL,
as well as field intensity. (b) an illustration of the electric field
intensity for different concentrations. As the concentration increases, the
electric fields are screened from the outer region and are concentrated in a
narrow region near the surface of the sphere. (c) the accumulated charge
inside the sphere, relative to the total surface charge, as a function of the
ratio of the radius to the Debye length

###  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:

$\displaystyle F[\\{c_{i}(\mathbf{r})\\}]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{V}{\rm d}^{3}\mathbf{r}\int_{V}{\rm
d}^{3}\mathbf{r}^{\prime}\frac{\rho(\mathbf{r})\rho(\mathbf{r}^{\prime})}{4\pi\varepsilon|\mathbf{r}-\mathbf{r}^{\prime}|}$
$\displaystyle+$ $\displaystyle k_{\mathrm{B}}T\int_{V}{\rm
d}^{3}\mathbf{r}\sum_{i=1}^{M}c_{i}(\mathbf{r})\left[\log\left(\lambda_{T}^{3}c_{i}(\mathbf{r})\right)-1\right]$
$\displaystyle\rho(\mathbf{r})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{M}ez_{i}c_{i}(\mathbf{r})+\rho_{\rm ext},$ (3)

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:

$\displaystyle F[\\{c_{i}\\}]$ $\displaystyle=$
$\displaystyle\frac{e^{2}}{2}\sum_{i=1}^{M}\sum_{j=1}^{M}z_{i}z_{j}c_{i}c_{j}\int_{V}{\rm
d}^{3}\mathbf{r}\int_{V}{\rm
d}^{3}\mathbf{r}^{\prime}\frac{1}{4\pi\varepsilon|\mathbf{r}-\mathbf{r}^{\prime}|}$
(4) $\displaystyle+$ $\displaystyle\int_{V}{\rm d}^{3}\mathbf{r}\int_{V}{\rm
d}^{3}\mathbf{r}^{\prime}\frac{\rho_{\rm ext}(\mathbf{r})\rho_{\rm
ext}(\mathbf{r}^{\prime})}{4\pi\varepsilon|\mathbf{r}-\mathbf{r}^{\prime}|}$
$\displaystyle+$ $\displaystyle\sum_{i=1}^{M}ez_{i}c_{i}\int_{V}{\rm
d}^{3}\mathbf{r}\int_{V}{\rm d}^{3}\mathbf{r}^{\prime}\frac{\rho_{\rm
ext}(\mathbf{r})}{8\pi\varepsilon|\mathbf{r}-\mathbf{r}^{\prime}|}$
$\displaystyle+$ $\displaystyle
Vk_{\mathrm{B}}T\sum_{i=1}^{M}c_{i}\left[\log\left(\lambda_{T}^{3}c_{i}\right)-1\right]$

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:

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle\frac{l_{{\rm
B}}}{eV}\int_{V}{\rm d}^{3}\mathbf{r}\int_{V}{\rm
d}^{3}\mathbf{r}^{\prime}\frac{1}{|\mathbf{r}-\mathbf{r}^{\prime}|},$
$\displaystyle\rho_{\rm ext}^{0}$ $\displaystyle=$ $\displaystyle\frac{l_{{\rm
B}}}{e\gamma V}\int_{V}{\rm d}^{3}\mathbf{r}\int_{V}{\rm
d}^{3}\mathbf{r}^{\prime}\frac{\rho_{\rm
ext}(\mathbf{r})}{|\mathbf{r}-\mathbf{r}^{\prime}|}.$ (5)

The free energy density is now simplified:

$\displaystyle\frac{\beta F[\\{c_{i}\\}]}{V}$
$\displaystyle=\sum_{i=1}^{M}\Bigg{[}\gamma
z_{i}c_{i}\left(\sum_{j=1}^{M}\frac{ez_{j}c_{j}}{2}+\rho_{\rm ext}^{0}\right)$
(6) $\displaystyle+c_{i}\left(\log(\lambda_{T}^{3}c_{i})-1\right)\Bigg{]},$

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

$\displaystyle\beta\mu_{i}$ $\displaystyle=$ $\displaystyle\gamma
z_{i}\left(\sum_{j=1}^{M}ez_{j}c_{j}+\rho_{\rm
ext}^{0}\right)+\log\left(\lambda_{T}^{3}c_{i}\right).$ (7)

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}$:

$\displaystyle\log\left(\lambda_{T}^{3}c_{i}^{0}\right)$ $\displaystyle=$
$\displaystyle\gamma z_{i}\left(\sum_{j=1}^{M}ez_{j}c_{j}+\rho_{\rm
ext}^{0}\right)+\log\left(\lambda_{T}^{3}c_{i}\right)$ $\displaystyle c_{i}$
$\displaystyle=$ $\displaystyle c_{i}^{0}\exp\left[\gamma
z_{i}\left(\sum_{j=1}^{M}ez_{j}c_{j}+\rho_{\rm ext}^{0}\right)\right].$ (8)

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

$\displaystyle\rho^{0}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{M}ez_{i}c_{i}^{0}\exp\left[-z_{i}\gamma(\rho_{\rm
ext}^{0}+\rho^{0})\right],$ (9)

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:

$\displaystyle\varepsilon\nabla^{2}\phi(\mathbf{r})=-\sum_{i=1}^{M}ez_{i}c_{i}^{0}\exp\left[-z_{i}e\beta\overline{\phi(\mathbf{r})}\right].$
(10)

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:

$\displaystyle\rho^{0}=-\frac{\rho_{\rm
ext}^{0}}{1+(e\gamma\sum_{i=1}^{M}z_{i}^{2}c_{i}^{0})^{-1}}.$ (11)

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})$:

$\displaystyle\int_{V}\frac{{\rm d}^{3}\mathbf{r}}{V}\int_{V}\frac{{\rm
d}^{3}\mathbf{r}^{\prime}}{\lambda_{D}^{2}|\mathbf{r}-\mathbf{r}^{\prime}|}\gg
1$ (12)

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.

Figure 3: Left: Phase diagram for electroneutrality, as a function of the
Gouy-Chapman length ($\lambda_{\rm GC}$) which is related to the surface
charge, and Debye screening length ($\lambda_{D}$), which is related to the
bulk ionic concentration. The color intensity indicates the ratio of total
accumulated charge inside a nanochannel to the surface charge and is obtained
by solving Eq. 16. In electroneutral systems, this ratio is $1$. The dashed
mark the transition to electroneutrality Eq. (19-20). Right: illustration of
the 4 regimes of the phase diagram. The blue circles mark the Debye screening
length.

###  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:

$\displaystyle\gamma^{\rm apprx}$ $\displaystyle=$ $\displaystyle\frac{l_{{\rm
B}}}{eV}\int_{V}{\rm d}^{3}\mathbf{r}\int_{V}{\rm
d}^{3}\mathbf{r}^{\prime}\frac{1}{|\mathbf{r}|}$ $\displaystyle=$
$\displaystyle\frac{l_{{\rm B}}}{e}\int_{0}^{R}\int_{-L/2}^{L/2}\frac{2\pi
r\,drdz}{\sqrt{r^{2}+z^{2}}}=\frac{2\pi R^{2}l_{{\rm
B}}}{e}\log\left(\frac{L}{2R}\right),$ $\displaystyle\rho_{\rm ext}^{0,{\rm
apprx}}$ $\displaystyle=$ $\displaystyle\frac{l_{{\rm B}}}{e\gamma^{\rm
apprx}V}\int_{V}\frac{{\rm d}^{3}\mathbf{r}}{|\mathbf{r}|}\int_{V}{\rm
d}^{3}\mathbf{r}^{\prime}q_{s}\delta(r^{\prime}-R).$ (13) $\displaystyle=$
$\displaystyle\frac{2\pi RLq_{s}}{V}=\frac{2q_{s}}{R}.$

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:

$\displaystyle L\gg(2R)\exp\left(\frac{2\lambda_{D}^{2}}{R^{2}}\right),$ (14)

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:

$\displaystyle\rho^{0}$ $\displaystyle=$ $\displaystyle
ec_{i}^{0}\exp\left[-\gamma(\rho_{\rm
ext}^{0}+\rho^{0})\right]-ec_{i}^{0}\exp\left[\gamma(\rho_{\rm
ext}^{0}+\rho^{0})\right]$ (15) $\displaystyle=$
$\displaystyle-2ec_{i}^{0}\sinh\left[\gamma(\rho_{\rm
ext}^{0}+\rho^{0})\right].$ $\displaystyle\approx$
$\displaystyle-2ec_{i}^{0}\sinh\left[\gamma^{\rm apprx}(\rho_{\rm ext}^{0,{\rm
apprx}}+\rho^{0})\right].$

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:

$\displaystyle\tilde{\rho}$ $\displaystyle=$
$\displaystyle-\xi\sinh(\tilde{\rho}_{\rm ext}+\tilde{\rho}),$ (16)

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:

$\displaystyle\tilde{\rho}_{\rm ext}$ $\displaystyle=$
$\displaystyle\gamma^{\rm apprx}\rho_{\rm ext}^{0,{\rm
apprx}}=\frac{2R}{\lambda_{\rm GC}}\log\left(\frac{L}{2R}\right)$
$\displaystyle\xi$ $\displaystyle=$ $\displaystyle 2ec_{0}\gamma^{\rm
apprx}=\frac{1}{2}\left(\frac{R}{\lambda_{D}}\right)^{2}\log\left(\frac{L}{2R}\right).$
(17)

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:

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

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:

$\displaystyle\lambda_{D}=R\sqrt{2\log\frac{L}{2R}}.$ (19)

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:

$\lambda_{D}=\sqrt{\lambda_{\rm
GC}R}\left(\frac{L}{2R}\right)^{R/2\lambda_{\rm GC}},$ (20)

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:

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

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.

Figure 4: The total ionic concentration in a nanotube, as a function of
concentration for different surface charges. Dashed lines show concentration
calculated for electroneutral pores.

###  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:

$\displaystyle\rho^{0}=\sum_{i=1}^{N}ez_{i}c_{i}^{0}\exp\left[-\beta
E_{i}-z_{i}\gamma(\rho^{0}+\rho_{\rm ext}^{0})\right].$ (22)

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}$:

$\displaystyle\tilde{\rho}_{\rm ext}$ $\displaystyle=$
$\displaystyle\frac{2R}{\lambda_{\rm
GC}}\log\left(\frac{L}{2R}\right)+\frac{E_{+}-E_{-}}{2k_{\mathrm{B}}T}$
$\displaystyle\xi$ $\displaystyle=$ $\displaystyle
2\left(\frac{R}{\lambda_{D}}\right)^{2}\log\left(\frac{L}{2R}\right)\exp\left(-\frac{E_{+}+E_{-}}{2k_{\mathrm{B}}T}\right).$
(23)

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:

$\varepsilon_{\rm in}\frac{\partial\phi_{\rm in}(r,z)}{\partial
r}=\varepsilon_{\rm out}\frac{\partial\phi_{\rm out}(r,z)}{\partial r}+q_{s}$
(24)

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

$\partial_{r}\phi_{\rm out}(R,z)\approx-\frac{\phi_{\rm out}(R,z)}{M_{L/R}R}.$
(25)

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

$M_{L/R}=\log\left(\frac{2L}{\pi R}\right)-\gamma_{\rm Euler},$ (26)

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:

$\partial_{r}\phi_{\rm in}(R)=\frac{q_{s}}{\varepsilon_{\rm
in}}-\frac{\varepsilon_{\rm out}}{\varepsilon_{\rm in}RM_{L/R}}\phi_{\rm
in}(R),$ (27)

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 $0$. 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.

Figure 5: Charge accumulated in a weakly charged pore as a function of inverse
Debye length, evaluated based on exact analytical solution (solid line),
numerical solution of the PB equation with the Robin boundary condition of Eq.
27 (dashed line), and the UEP model (dash-dotted line, Eq. 16. Results are
compared with COMSOL simulation (circles). The pore dimensions are: $R=5$nm
and $L=1\mu$m.

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):

$\displaystyle S$ $\displaystyle=$ $\displaystyle-k_{\rm
B}\sum_{i=1}^{N}\left[p_{i}^{+}\log p_{i}^{+}+p_{i}^{-}\log p_{i}^{-}\right]$
(28) $\displaystyle=$ $\displaystyle-N_{s}k_{\rm
B}\left[\left(\frac{N_{s}+1}{2N_{s}}\right)\log\left(\frac{N_{s}+1}{2N_{s}}\right)\right.$
$\displaystyle+$
$\displaystyle\left.\left(\frac{N_{s}-1}{2N_{s}}\right)\log\left(\frac{N_{s}-1}{2N_{s}}\right)\right]+(N-N_{s})k_{\rm
B}\log 2$ $\displaystyle\approx$ $\displaystyle k_{\rm B}\left(N\log
2-\frac{1}{2N_{s}}\right),$

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:

$\displaystyle\Delta S=-\frac{dk_{\rm B}}{4A_{d}c\lambda_{s}^{d}}.$ (29)

The electrostatic energy of a uniformly distributed screening cloud equals:

$\displaystyle U$ $\displaystyle=$
$\displaystyle-\frac{e^{2}A_{d}}{4\pi\varepsilon
V}\int_{a}^{\lambda_{s}}r^{d-2}{\rm d}r$ (30) $\displaystyle=$
$\displaystyle-\frac{e^{2}}{4\pi\varepsilon\lambda_{s}}\times\begin{cases}\frac{d}{d-1},&d>1.\\\
\log\frac{\lambda_{s}}{a},&d=1,\end{cases}$

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:

$\displaystyle\lambda_{s}\propto\begin{cases}(cl_{{\rm
B}})^{\frac{1}{1-d}},&d>1\\\ a\exp(\frac{1}{8cl_{{\rm B}}}),&d=1,\end{cases}$
(31)

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.

Figure 6: Illustration of the standard Debye screening in 3d(a), and the
extended screening length of 1d(b), according to Eq. 31.

##  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:

$\displaystyle U$ $\displaystyle=$ $\displaystyle\frac{1}{8\pi}\sum_{i\neq
j}\frac{e^{2}(p^{+}_{i}-p^{-}_{i})(p^{+}_{j}-p^{-}_{j})}{a\varepsilon|i-j|}$
(32) $\displaystyle+$ $\displaystyle\sum_{i}e(p^{+}_{i}-p^{-}_{i})\phi^{\rm
ext}_{i},$

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:

$U=\frac{e^{2}}{8\pi
a\varepsilon}\left(\mathbf{q}^{T}\tilde{\Phi}\mathbf{q}\right)+e\mathbf{q}^{T}\phi^{\rm
ext}.$ (33)

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

$\displaystyle S$ $\displaystyle=$ $\displaystyle-k_{\rm
B}\sum_{i}\left[p^{+}_{i}\log p^{+}_{i}+p^{-}_{i}\log p^{-}_{i}+\right.$ (34)
$\displaystyle\left.(1-p^{+}_{i}-p_{i}^{-})\log(1-p^{+}_{i}-p_{i}^{-})\right].$

Figure 7: Solution of the 1D-ring equation with uniform surface charge (Eq.
V.1). The solution, $q_{i}=q_{0}+q_{1}(-1)^{i}$, has an average part ($q_{0}$,
left figure) that screens the external charge ($Q$) and an oscillatory part
($q_{1}$, right figure). Solutions are presented as a function of the bulk
Debye length ($\lambda_{D}$) and the Gouy-Chapman length ($\lambda_{\rm GC}$).
The dashed vertical line is the critical Debye length, calculated by Eq. 41.

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$):

$\displaystyle\mu_{\pm}=k_{\mathrm{B}}T\log\left(\frac{\mathbf{p}_{\pm}}{1-\mathbf{p}_{+}-\mathbf{p}_{-}}\right)\pm\frac{e^{2}\tilde{\Phi}\mathbf{q}}{4\pi
a\varepsilon}\pm e\phi^{\rm ext}.$

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)):

$\displaystyle\mathbf{q}=-\frac{2\Lambda\sinh\left[\frac{l_{B}}{a}\tilde{\Phi}\mathbf{q}+e\beta\phi^{\rm
ext}\right]}{1+2\Lambda\cosh\left[\frac{l_{B}}{a}\tilde{\Phi}\mathbf{q}+e\beta\phi^{\rm
ext}\right]},$ (36)

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.

Figure 8: Charge distribution around a central ion and the resulting activity
coefficient, for different concentrations. (a-f) Charge distribution around a
central ion for six different coupling strengths ($\chi$). For weak coupling ,
the central ion is screened only by oppositely charged ions (e-f). As the
coupling increases, over screening and oscillations are observed. (g) Ionic
activity in a nanopore, based on the exact 1D solution for the mean field
equation (Eq. V.2, solid blue line) and the dilute limit approximation (Eq.
48, dashed red line).

###  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:

$\displaystyle H_{L/a}\approx\log\left(\frac{L}{a}\right)+\gamma_{\rm
Euler}+\frac{a}{2L}+\ldots,$ (37)

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:

$\displaystyle
q=-\frac{2\Lambda\sinh\left[2\frac{l_{B}}{a}H_{\frac{L}{a}}\left(q+Q\right)\right]}{1+2\Lambda\cosh\left[2\frac{l_{B}}{a}H_{\frac{L}{a}}\left(q+Q\right)\right]},$
(38)

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:

$\displaystyle q_{0}-q_{1}$ $\displaystyle=$
$\displaystyle-\frac{2\Lambda\sinh\left[-2q_{1}\frac{l_{B}}{a}\log(2)+2\frac{l_{B}}{a}H_{\frac{L}{a}}\left(q_{0}+Q\right)\right]}{1+2\Lambda\cosh\left[2q_{1}\frac{l_{B}}{a}\log(2)+2\frac{l_{B}}{a}H_{\frac{L}{a}}\left(q_{0}+Q\right)\right]}$
$\displaystyle q_{0}+q_{1}$ $\displaystyle=$
$\displaystyle-\frac{2\Lambda\sinh\left[2q_{1}\frac{l_{B}}{a}\log(2)+2\frac{l_{B}}{a}H_{\frac{L}{a}}\left(q_{0}+Q\right)\right]}{1+2\Lambda\cosh\left[-2q_{1}\frac{l_{B}}{a}\log(2)+2\frac{l_{B}}{a}H_{\frac{L}{a}}\left(q_{0}+Q\right)\right]},$

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:

$\displaystyle q_{1}$ $\displaystyle=$
$\displaystyle\frac{2\Lambda\sinh\left(2q_{1}\frac{l_{B}}{a}\log
2\right)}{1+2\Lambda\cosh\left(2q_{1}\frac{l_{B}}{a}\log 2\right)}$

The RHS of Eq. (V.1) is monotonically increasing, starting from $0$, 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:

$\Lambda_{\rm cr}=\frac{a}{4l_{B}\log 2-2a}.$ (41)

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:

$\displaystyle\mathbf{q}=-\frac{2\Lambda\sinh\left[\frac{l_{B}}{a}\left(\tilde{\Phi}\mathbf{q}+\alpha\delta_{i,0}\right)\right]}{1+2\Lambda\cosh\left[\frac{l_{B}}{a}\left(\tilde{\Phi}\mathbf{q}+\alpha\delta_{i,0}\right)\right]},$
(42)

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

$\displaystyle q_{i}$ $\displaystyle=$
$\displaystyle-\alpha\left(\tilde{\Phi}+\chi^{-1}\right)_{i,j}^{-1}\delta_{j,0},$
(43)

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:

$\displaystyle q_{i}$ $\displaystyle=$
$\displaystyle\frac{\alpha}{4\pi}\int_{2\pi}{\rm d}\omega\frac{\cos(\omega
i)}{2\log\left|2\sin\frac{\omega}{2}\right|-\chi^{-1}},$ (44)

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 $0$:

$\displaystyle\sum_{i=0}^{\infty}q_{i}=\alpha\lim_{\omega\rightarrow
0}\frac{1}{2\log\left|2\sin\frac{\omega}{2}\right|-\chi^{-1}}\rightarrow 0.$
(45)

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:

$\displaystyle\frac{1}{2\log\left|\Delta\omega\right|-\chi^{-1}}=\frac{1/2}{2\log
2-\chi^{-1}}\rightarrow\Delta\omega=4{\rm e}^{-\frac{a}{4\Lambda l_{B}}}.$

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:

$\displaystyle\beta\mu^{\rm
ex}=\frac{l_{B}}{a}\sum_{i}\frac{q_{i}}{i}=\frac{\alpha l_{{\rm B}}}{2\pi
a}\int_{2\pi}{\rm
d}\omega\frac{\log\left|2\sin\frac{\omega}{2}\right|}{2\log\left|2\sin\frac{\omega}{2}\right|-\chi^{-1}},$

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:

$\beta\mu^{\rm ex}\approx-\frac{2\pi^{2}l_{B}^{2}}{3a^{2}}\Lambda.$ (48)

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:

$\displaystyle\frac{\mu^{1D}}{\mu^{3D}}=\frac{\pi^{2}}{24}\frac{a}{\lambda_{D}},$
(49)

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.

Figure 9: Conductance of a KCl solution as a function of concentration, inside
CNT with varying surface charge and radii, fitted according to our model
(dashed lines). The experimental data (circles) was adapted from Secchi  _et
al._ (2016), where the surface charge was controlled by changing the pH . (a)
- A $3.5$nm wide pore, fitted barrier energy of $4.6k_{\mathrm{B}}T$ and 4
different surface charges, from bottom to top: $-3{\rm mC/m}^{2}$ (black),
$-1.6{\rm mC/m}^{2}$ (red), $3.4{\rm mC/m}^{2}$ (green), $5.2{\rm mC/m}^{2}$
(blue). (b) - A $10$nm wide pore, with only one surface charge: $12{\rm
mC/m}^{2}$. (c) - A $14$nm wide pore, with $3$ different surface charge,
fitted to: $63{\rm mC/m}^{2}$ (black), $213{\rm mC/m}^{2}$ (red) and $417{\rm
mC/m}^{2}$ (blue). (d)- A $34$nm wide pore, with fitted surface charge of:
$-7{\rm mC/m}^{2}$ (black), $55{\rm mC/m}^{2}$ (red) and $110{\rm mC/m}^{2}$
(blue).

##  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:

$\displaystyle G=e^{2}\upmu_{D}\frac{\pi R^{2}}{L}c,$ (50)

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.

Figure 10: Conductance of sub-nanometer pores as a function of ionic
concentration. Main figure: Conductance of a $3.5$Åwide CNT porin as a
function of KCl concentration (filled squares), adapted from Tunuguntla  _et
al._ (2017) and fitted according to our model (dashed lines). The bottom curve
(blue) has pH 3.5, which corresponds to zero charge, while the upper (black)
curve has pH 7.5, and is fitted to a surface charge of $6{\rm mC/m}^{2}$.
Inset: conductance in Gramicidin-A channel as a function of NaCl concentration
(filled squares, adapted from Finkelstein and Andersen (1981)) and its fit
(dashed lines). The fitted surface charge is $13{\rm mC/m}^{2}$.

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).

Figure 11: Fitted mobilities. The mobilities of $6$ datasets, Gramicidin-A,
CNT Porin, and 4 different CNT experiments, were fitted according to Eq. 50.
Experimental data and fits are shown if Figs. 7-8. Each pH was fitted
separately, and within each experiment the bins are ordered from low pH to
high. The bulk value for KCl mobility (last column) was taken from Lide (2004)

.

##  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.

## References

  * Hille (2001) B. Hille,  _Ionic channels in excitable membranes_ (sinauer associates Inc, Sunderland, MA, 2001).
  * Doyle  _et al._ (1998) D. A. Doyle, J. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L. Cohen, B. T. Chait, and R. MacKinnon, Science 280, 69 (1998).
  * Berneche and Roux (2001) S. Berneche and B. Roux, Nature 414, 73 (2001).
  * Hess and Tsien (1984) P. Hess and R. W. Tsien, Nature 309, 453 (1984).
  * Werner  _et al._ (2001) C. Werner, R. Zimmermann, and T. Kratzmüller, Colloids Surf. A 192, 205 (2001).
  * Vilatela and Marcilla (2015) J. J. Vilatela and R. Marcilla, Chem. Mater. 27, 6901 (2015).
  * Yang  _et al._ (2003) J. Yang, F. Lu, L. W. Kostiuk, and D. Y. Kwok, J. Micromechanics Microengineering 13, 963 (2003).
  * Kim  _et al._ (2010a) S. J. Kim, S. H. Ko, K. H. Kang, and J. Han, Nat. Nanotechnol. 5, 297 (2010a).
  * Kim  _et al._ (2010b) D.-K. Kim, C. Duan, Y.-F. Chen, and A. Majumdar, Microfluidics Nanofluidics 9, 1215 (2010b).
  * Mani and Bazant (2011) A. Mani and M. Z. Bazant, Phys. Rev. E 84, 061504 (2011).
  * Peters  _et al._ (2016) P. B. Peters, R. van Roij, M. Z. Bazant, and P. M. Biesheuvel, Phys. Rev. E 93, 053108 (2016).
  * Schmuck and Bazant (2015) M. Schmuck and M. Z. Bazant, SIAM J. Appl. Math 75, 1369 (2015).
  * Nonner and Eisenberg (1998) W. Nonner and B. Eisenberg, Biophys. J. 75, 1287 (1998).
  * Tedesco  _et al._ (2016) M. Tedesco, H. V. M. Hamelers, and P. M. Biesheuvel, J. Memb. Sci. 510, 370 (2016).
  * Lee  _et al._ (2010) C. Y. Lee, W. Choi, J.-H. Han, and M. S. Strano, Science 329, 1320 (2010).
  * Agrawal  _et al._ (2017) K. V. Agrawal, S. Shimizu, L. W. Drahushuk, D. Kilcoyne, and M. S. Strano, Nat. Nanotechnol. 12, 267 (2017).
  * Siria  _et al._ (2013) A. Siria, P. Poncharal, A.-L. Biance, R. Fulcrand, X. Blase, S. T. Purcell, and L. Bocquet, Nature 494, 455 (2013).
  * Guan and Reed (2012) W. Guan and M. A. Reed, Nano Lett. 12, 6441 (2012).
  * Daiguji  _et al._ (2004) H. Daiguji, P. Yang, A. J. Szeri, and A. Majumdar, Nano Lett. 4, 2315 (2004).
  * Faucher  _et al._ (2019) S. J. Faucher, N. R. Aluru, M. Z. Bazant, D. Blankschtein, A. H. Brozena, J. Cumings, J. P. de Souza, M. Elimelech, R. Epsztein, J. T. Fourkas, _et al._ , The Journal of Physical Chemistry C (2019).
  * Staverman (1952) A. J. Staverman, J. Chem. Soc. Faraday Trans. 48, 176 (1952).
  * Gross and Osterle (1968) R. J. Gross and J. Osterle, J. Chem. Phys. 49, 228 (1968).
  * Fair and Osterle (1971) J. C. Fair and J. F. Osterle, J. Chem. Phys. 54, 3307 (1971).
  * Sasidhar and Ruckenstein (1981) V. Sasidhar and E. Ruckenstein, J. Colloid Interface Sci. 82, 439 (1981).
  * Yao and Santiago (2003) S. Yao and J. G. Santiago, J. Colloid Interface Sci. 268, 133 (2003).
  * Petsev and Lopez (2006) D. N. Petsev and G. P. Lopez, J. Colloid Interface Sci. 294, 492 (2006).
  * Kim  _et al._ (2010c) S. J. Kim, Y.-A. Song, and J. Han, Chem. Soc. Rev. 39, 912 (2010c).
  * Chang and Yang (2011) C.-C. Chang and R.-J. Yang, Appl. Phys. Lett. 99, 083102 (2011).
  * Yan  _et al._ (2013) Y. Yan, Q. Sheng, C. Wang, J. Xue, and H.-C. Chang, J. Phys. Chem. C 117, 8050 (2013).
  * Zwolinski  _et al._ (1949) B. J. Zwolinski, H. Eyring, and C. E. Reese, J. Phys. Chem. 53, 1426 (1949).
  * Almers and McCleskey (1984) W. Almers and E. McCleskey, J. Physiol. 353, 585 (1984).
  * Sather and McCleskey (2003) W. A. Sather and E. W. McCleskey, Annu. Rev. Physiol. 65, 133 (2003).
  * Secchi  _et al._ (2016) E. Secchi, A. Niguès, L. Jubin, A. Siria, and L. Bocquet, Phys. Rev. Lett. 116, 154501 (2016).
  * Biesheuvel and Bazant (2016) P. M. Biesheuvel and M. Z. Bazant, Phys. Rev. E 94, 050601(R) (2016).
  * Tunuguntla  _et al._ (2017) R. H. Tunuguntla, R. Y. Henley, Y.-C. Yao, T. A. Pham, M. Wanunu, and A. Noy, Science 357, 792 (2017).
  * Sørensen and Sloth (1992) T. S. Sørensen and P. Sloth, J. Chem. Soc. Faraday Trans. 88, 571 (1992).
  * Rivera and Sørensen (1994) S. R. Rivera and T. S. Sørensen, Mol. Simul. 13, 115 (1994).
  * Lo and Chan (1995) W. Y. Lo and K.-Y. Chan, Mol. Phys. 86, 745 (1995).
  * Lee and Chan (1997) M. Lee and K.-Y. Chan, Chem. Phys. Lett. 275, 56 (1997).
  * Lee  _et al._ (1999) M. Lee, K.-Y. Chan, D. Nicholson, and S. Zara, Chem. Phys. Lett. 307, 89 (1999).
  * Boda  _et al._ (2000) D. Boda, D. D. Busath, D. Henderson, and S. Sokołowski, J. Phys. Chem. B 104, 8903 (2000).
  * Boda  _et al._ (2002) D. Boda, D. D. Busath, B. Eisenberg, D. Henderson, and W. Nonner, Phys. Chem. Chem. Phys. 4, 5154 (2002).
  * Boda  _et al._ (2006) D. Boda, M. Valiskó, B. Eisenberg, W. Nonner, D. Henderson, and D. Gillespie, J. Chem. Phys. 125, 034901 (2006).
  * Luo  _et al._ (2015) Z.-X. Luo, Y.-Z. Xing, Y.-C. Ling, A. Kleinhammes, and Y. Wu, Nat. Commun. 6, 6358 (2015).
  * Lozada-Cassou  _et al._ (1996a) M. Lozada-Cassou, W. Olivares, and B. Sulbarán, Phys. Rev. E 53, 522 (1996a).
  * Lozada-Cassou  _et al._ (1996b) M. Lozada-Cassou, W. Olivares, B. Sulbarán, and Y. Jiang, Physica A 231, 197 (1996b).
  * Colla  _et al._ (2016) T. Colla, M. Girotto, A. P. dos Santos, and Y. Levin, J. Chem. Phys. 145, 094704 (2016).
  * Jordan  _et al._ (1989) P. Jordan, R. Bacquet, J. McCammon, and P. Tran, Biophys. J. 55, 1041 (1989).
  * McLaren-Young-Sommerville (1958) D. McLaren-Young-Sommerville, (1958).
  * Schlögl (1955) R. Schlögl, Z. physik. Chem 3, 73 (1955).
  * Sonin (1976) A. A. Sonin, in  _Charged Gels and Membranes_ (Springer, 1976) pp. 255–265.
  * Cwirko and Carbonell (1989) E. H. Cwirko and R. Carbonell, J. Colloid Interface Sci. 129, 513 (1989).
  * Verbrugge and Hill (1990) M. W. Verbrugge and R. F. Hill, J. Electrochem. Soc. 137, 886 (1990).
  * Bowen and Welfoot (2002) W. R. Bowen and J. S. Welfoot, Chem. Eng. Sci. 57, 1121 (2002).
  * Cowan  _et al._ (1937) S. Cowan, J. Schulman, T. Teorell, W. Wilbrandt, A. Hill, J. Bikerman, R. Dean, and K. H. Meyer, J. Chem. Soc. Faraday Trans. 33, 1051 (1937).
  * Teorell (1953) T. Teorell, Prog. Biophys. Biophys. Chem. 3, 305 (1953).
  * Meyer and Sievers (1936) K. H. Meyer and J.-F. Sievers, Helv. Chim. Acta 19, 649 (1936).
  * Yaroshchuk (2012) A. Yaroshchuk, Microfluid Nanofluidics 12, 615 (2012).
  * Dydek and Bazant (2013) E. V. Dydek and M. Z. Bazant, AIChE J 59, 3539 (2013).
  * Epsztein  _et al._ (2019) R. Epsztein, E. Shaulsky, M. Qin, and M. Elimelech, J. Memb. Sci. (2019).
  * Sahu and Zwolak (2017) S. Sahu and M. Zwolak, Nanoscale 9, 11424 (2017).
  * Kamenev  _et al._ (2006) A. Kamenev, J. Zhang, A. Larkin, and B. I. Shklovskii, Physica A 359, 129 (2006).
  * Zhang  _et al._ (2005) J. Zhang, A. Kamenev, and B. I. Shklovskii, Phys. Rev. Lett. 95, 148101 (2005).
  * Zhang  _et al._ (2006) J. Zhang, A. Kamenev, and B. I. Shklovskii, Phys. Rev. E 73, 051205 (2006).
  * Levin (2006) Y. Levin, Europhys. Lett. 76, 163 (2006).
  * Bordin  _et al._ (2012) J. R. Bordin, A. Diehl, M. C. Barbosa, and Y. Levin, Phys. Rev. E 85, 031914 (2012).
  * Edwards and Lenard (1962) S. Edwards and A. Lenard, J. Math. Phys. 3, 778 (1962).
  * Lenard (1961) A. Lenard, J. Math. Phys. 2, 682 (1961).
  * Démery  _et al._ (2012) V. Démery, D. S. Dean, T. C. Hammant, R. R. Horgan, and R. Podgornik, Europhys. Lett. 97, 28004 (2012).
  * Bikerman (1942) J. J. Bikerman, Philos. Mag. 33, 384 (1942).
  * Borukhov  _et al._ (1997) I. Borukhov, D. Andelman, and H. Orland, Phys. Rev. Lett. 79, 435 (1997).
  * Kornyshev (2007) A. A. Kornyshev, J. Phys. Chem. B 111, 5545 (2007).
  * Kilic  _et al._ (2007) M. S. Kilic, M. Z. Bazant, and A. Ajdari, Phys. Rev. E 75, 021502 (2007).
  * Bazant  _et al._ (2009) M. Z. Bazant, M. S. Kilic, B. D. Storey, and A. Ajdari, Adv. Colloid. Interface. Sci. 152, 48 (2009).
  * Dukhin  _et al._ (2005) S. S. Dukhin, R. Zimmermann, and C. Werner, J. Colloid Interface Sci. 286, 761 (2005).
  * Lide (2004) D. R. Lide,  _CRC handbook of chemistry and physics_ , Vol. 85 (CRC press, 2004).
  * Holt  _et al._ (2006) J. K. Holt, H. G. Park, Y. Wang, M. Stadermann, A. B. Artyukhin, C. P. Grigoropoulos, A. Noy, and O. Bakajin, Science 312, 1034 (2006).
  * Whitby  _et al._ (2008) M. Whitby, L. Cagnon, M. Thanou, and N. Quirke, Nano Lett. 8, 2632 (2008).
  * Finkelstein and Andersen (1981) A. Finkelstein and O. S. Andersen, J. Membr. Biol. 59, 155 (1981).
  * Donev  _et al._ (2019) A. Donev, A. J. Nonaka, C. Kim, A. L. Garcia, and J. B. Bell, Phys. Rev. Fluids 4, 043701 (2019).
  * Balu and Khair (2018) B. Balu and A. S. Khair, Soft Matter 14, 8267 (2018).
  * Evans III  _et al._ (1961) R. Evans III, G. Watson, and E. Mason, J. Chem. Phys. 35, 2076 (1961).
  * Fu  _et al._ (2015) Y. Fu, Y. Jiang, S. Poizeau, A. Dutta, A. Mohanram, J. D. Pietras, and M. Z. Bazant, J. Electrochem. Soc. 162, F613 (2015).
  * Boda  _et al._ (2007a) D. Boda, W. Nonner, M. Valiskó, D. Henderson, B. Eisenberg, and D. Gillespie, Biophys. J. 93, 1960 (2007a).
  * Boda  _et al._ (2007b) D. Boda, M. Valiskó, B. Eisenberg, W. Nonner, D. Henderson, and D. Gillespie, Phys. Rev. Lett. 98, 168102 (2007b).
  * Esfandiar  _et al._ (2017) A. Esfandiar, B. Radha, F. Wang, Q. Yang, S. Hu, S. Garaj, R. Nair, A. Geim, and K. Gopinadhan, Science 358, 511 (2017).
  * Duan and Majumdar (2010) C. Duan and A. Majumdar, Nat. Nanotechnol. 5, 848 (2010).
  * Green  _et al._ (2014) Y. Green, S. Shloush, and G. Yossifon, Phys. Rev. E 89, 043015 (2014).
  * Green  _et al._ (2016) Y. Green, R. Eshel, S. Park, and G. Yossifon, Nano Lett. 16, 2744 (2016).
  * Corry (2011) B. Corry, Energy Environ. Sci. 4, 751 (2011).
  * Chan  _et al._ (2013) W.-F. Chan, H.-y. Chen, A. Surapathi, M. G. Taylor, X. Shao, E. Marand, and J. K. Johnson, ACS Nano 7, 5308 (2013).
  * Gravelle  _et al._ (2013) S. Gravelle, L. Joly, F. Detcheverry, C. Ybert, C. Cottin-Bizonne, and L. Bocquet, Proc. Natl. Acad. Sci. U.S.A. 110, 16367 (2013).
  * Barcilon (1992) V. Barcilon, SIAM J. Appl. Math 52, 1391 (1992).

##  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:

$\displaystyle\begin{cases}\nabla^{2}\phi_{\rm
in}(\mathbf{r})=\lambda_{D}^{-2}\phi_{\rm in}(\mathbf{r})&r<R\\\
\nabla^{2}\phi_{\rm out}(\mathbf{r})=0&r>R.\end{cases}$ (51)

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

$\displaystyle\left[\phi_{\rm out}(R,z)-\phi_{\rm in}(R,z)\right]$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle\left[\varepsilon_{\rm
out}\partial_{r}\phi_{\rm out}(R,z)-\varepsilon_{\rm in}\partial_{r}\phi_{\rm
in}(R,z)\right]$ $\displaystyle=$ $\displaystyle-q_{s}(0<z<L)$
$\displaystyle\left.\phi_{\rm out}(\mathbf{r})\right|_{r\rightarrow\infty}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle\phi_{\rm
in}(R,z=0)=\phi_{\rm out}(R,z=0)$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\phi_{\rm in}(R,z=L)=\phi_{\rm out}(R,z=L)$ $\displaystyle=$
$\displaystyle 0.$ (52)

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:

$\displaystyle\phi_{\rm in}(r,z)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}A_{n}\sin(\omega_{n}z)I_{0}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}r)$
$\displaystyle\phi_{\rm out}(r,z)$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}B_{n}\sin(\omega_{n}z)K_{0}(\omega_{n}r),$
(53)

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:

$\left[\varepsilon_{\rm out}\partial_{r}\phi_{\rm out}(R,z)-\varepsilon_{\rm
in}\partial_{r}\phi_{\rm in}(R,z)\right]=-q_{s}\sum_{n\,{\rm
odd}}\frac{4\sin(\omega_{n}z)}{\omega_{n}L}.$ (54)

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

$\displaystyle A_{n}I_{0}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}R)$
$\displaystyle=$ $\displaystyle B_{n}K_{0}(\omega_{n}R),$
$\displaystyle\varepsilon_{\rm
in}\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}I_{1}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}R)$
$\displaystyle=$ $\displaystyle-B_{n}\varepsilon_{\rm
out}\omega_{n}K_{1}(\omega_{n}R)+\frac{4q_{s}}{\omega_{n}L}.$

Finally, we get:

$\displaystyle\phi_{\rm in}(r,z)$ $\displaystyle=$
$\displaystyle\frac{2q_{s}}{L}\sum_{n\,{\rm
odd}}\frac{1}{\omega_{n}}\frac{K_{0}(\omega_{n}R)\sin(\omega_{n}z)I_{0}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}r)}{\varepsilon_{\rm
out}\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}K_{0}(\omega_{n}R)I_{1}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}R)+\varepsilon_{\rm
out}\omega_{n}I_{0}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}R)K_{1}(\omega_{n}R)}.$
$\displaystyle\phi_{\rm out}(r,z)$ $\displaystyle=$
$\displaystyle\frac{2q_{s}}{L}\sum_{n\,{\rm
odd}}\frac{1}{\omega_{n}}\frac{K_{0}(\omega_{n}r)\sin(\omega_{n}z)I_{0}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}R)}{\varepsilon_{\rm
out}\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}K_{0}(\omega_{n}R)I_{1}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}R)+\varepsilon_{\rm
out}\omega_{n}I_{0}(\sqrt{\omega_{n}^{2}+\kappa_{D}^{2}}R)K_{1}(\omega_{n}R)}.$

##  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:

$\phi_{\rm
out}(R,z)=\sum_{n=1}^{\infty}B_{n}\sin\left(n\pi\frac{z}{L}\right)K_{0}\left(n\pi\frac{R}{L}\right),$
(57)

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:

$\displaystyle\phi_{\rm out}(R,z)$ $\displaystyle=$
$\displaystyle-\sum_{n=1}^{\infty}B_{n}\sin\left(n\pi\frac{z}{L}\right)\left(\log\left(\frac{n\pi
R}{2L}\right)+\gamma_{\rm Euler}\right)$ (58) $\displaystyle+$ $\displaystyle
O\left(\left(\frac{R}{L}\right)^{2}\log\left(\frac{R}{L}\right)\right),$

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

$\log\left(\frac{n\pi R}{2L}\right)=\log\left(\frac{\pi
R}{2L}\right)+\log(n)\approx\log\left(\frac{\pi R}{2L}\right).$ (59)

which leads to the following approximation for the electrostatic potential:

$\displaystyle\phi_{\rm out}(R)$ $\displaystyle\approx$ $\displaystyle
M_{L/R}\sum_{n=1}^{\infty}B_{n}\sin(\omega_{n}z),$ (60)

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

$M_{L/R}=\log\left(\frac{2L}{\pi R}\right)-\gamma.$ (61)

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:

$\partial_{r}\phi_{\rm
out}(R,z)=-\sum_{n=1}^{\infty}B_{n}\omega_{n}K_{1}(\omega_{n}R)\sin(\omega_{n}z).$
(62)

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

$\partial_{r}\phi_{\rm
out}(R,z)\approx-\sum_{n=1}^{\infty}B_{n}\frac{\sin(\omega_{n}z)}{R}=-\frac{\phi_{\rm
out}(R,z)}{M_{L/R}R}.$ (63)

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

$\partial_{r}\phi_{\rm in}(R)=\frac{q_{s}}{\varepsilon_{\rm
in}}-\frac{\varepsilon_{\rm out}}{\varepsilon_{\rm in}}\frac{\phi_{\rm
in}(R)}{RM_{L/R}}.$ (64)