# Scaling for rectification of bipolar nanopores as a function of a modified
Dukhin number: the case of 1:1 electrolytes

Dávid Fertig  Zsófia Sarkadi  Mónika Valiskó  Dezső Boda<EMAIL_ADDRESS>1Center for Natural Sciences, University of Pannonia, P. O. Box 158, H-8201
Veszprém, Hungary

###### Abstract

The scaling behavior for the rectification of bipolar nanopores is studied
using the Nernst-Planck equation coupled to the Local Equilibrium Monte Carlo
method. The bipolar nanopore’s wall carries $\sigma$ and $-\sigma$ surface
charge densities in its two half regions axially. Scaling means that the
device function (rectification) depends on the system parameters (pore length,
$H$, pore radius, $R$, concentration, $c$, voltage, $U$, and surface charge
density, $\sigma$) via a single scaling parameter that is a smooth analytical
function of the system parameters. Here, we suggest using a modified Dukhin
number,
$\mathrm{mDu}=|\sigma|l_{\mathrm{B}}^{*}\lambda_{\mathrm{D}}HU/(RU_{0})$,
where $l_{\mathrm{B}}^{*}=8\pi l_{\mathrm{B}}$, $l_{\mathrm{B}}$ is the
Bjerrum length, $\lambda_{\mathrm{D}}$ is the Debye length, and $U_{0}$ is a
reference voltage. We show how scaling depends on $H$, $U$, and $\sigma$ and
through what mechanisms these parameters influence the pore’s behavior.

##  I Introduction

Nanopores are channels in membranes connecting two baths that have radii
comparable to the screening length, $\lambda$, of the electrolyte whose ions
are transported through the pore. If the nanopore is asymmetric in the axial
dimension along the pore, the ionic current is larger at one sign of the
voltage (ON state) than at the opposite sign (OFF state), that is, the
nanopore rectifies. In this paper, we are interested in the case of
electrostatically asymmetric cylindrical nanopores, specifically, bipolar
nanopores daiguji_nl_2005 ; constantin_pre_2007 ; karnik_nl_2007 ;
vlassiouk_nl_2007 ; kalman_am_2008 ; vlassiouk_acsnanno_2008 ; yan_nl_2009 ;
cheng_acsnano_2009 ; nguyen_nt_2010 ; guo_acccr_2013 ; hato_pccp_2017 ;
matejczyk_jcp_2017  whose wall carries positive surface charge, $\sigma$, on
one side and negative surface charge, $-\sigma$, on the other side (Fig. 1B,
$\sigma$ is always positive in this work).

In a previous publication fertig_jpcc_2019 , we showed that, for fixed surface
charge density, $\sigma$, pore length, $H$, and voltage, $U$, rectification
defined as

$r=\frac{I^{\mathrm{\mathrm{ON}}}}{I^{\mathrm{OFF}}}$ (1)

scales with the parameter $\xi=R/(\lambda z_{\mathrm{if}})$ where $R$ is the
radius of the cylindrical nanopore, and $z_{\mathrm{if}}=\sqrt{z_{+}|z_{-}|}$
with $z_{+}$ and $z_{-}$ being the valences of the cations and the anions,
respectively. In Eq. 1, $I^{\mathrm{ON}}$ and $I^{\mathrm{OFF}}$ are the
magnitudes of the total currents at the forward and reverse biased voltages,
respectively.

By scaling, we mean that the device function (rectification, in this case)
depends on the input variables via a single parameter that is a well-defined
analytic function of the input variables:

$r=r\left[\xi(R,c,z_{+},z_{-},\dots)\right].$ (2)

The $\xi$ parameter may depend on concentration, $c$, via the screening length
for which the Debye length

$\lambda_{\mathrm{D}}=\left(\dfrac{ce^{2}}{\epsilon_{0}\epsilon
kT}\sum_{i}z_{i}^{2}\nu_{i}\right)^{-1/2},$ (3)

is an obvious choice. In Eq. 3, $k$ is Boltzmann’s constant, $T$ is the
absolute temperature ($298.15$ K in this work), $e$ is the elementary charge,
$\nu_{i}$ is the stoichiometric coefficient of ionic species $i$, $c$ is the
electrolyte concentration that is related to the ionic bath concentrations via
$c=c_{+}/|z_{-}|=c_{-}/z_{+}=c_{+}/\nu_{-}=c_{-}/\nu_{+}$, $\epsilon$ is the
dielectric constant of the solvent ($78.45$ in this work), and $\epsilon_{0}$
is the permittivity of vacuum. Another choice for the screening length is the
one provided by the Mean Spherical Approximation (MSA) blum_mp_1975 ;
blum_jcp_1977 ; nonner_bj_2000 , denoted by $\lambda_{\mathrm{MSA}}$.

The basic idea behind the scaling behavior is the following. Double layers
(DL) are formed at the charged wall of the nanopore and extend into the middle
of the pore in the radial dimension. They overlap in the centerline if
$\lambda_{\mathrm{D}}/R\gg 1$. If they overlap, coions are excluded from the
respective regions. The concept of coion is linked to the sign of the surface
charge, namely, in the positively charged region ($\sigma$, “p” region) the
cation is the coion, while in the negatively charged region ($-\sigma$, “n”
region) the anion is the coion. Exclusion of an ionic species results in
depletion zones for that ionic species. A depletion zone acts as a low-
resistance element for the given ionic species along the pore axis if we
imagine the consecutive segments of the pore as resistors connected in series.

Figure 1: The geometries of the nanopores under consideration. (A) The
uniformly charged nanopore’s device function is selectivity. The scaling
behavior of this pore was studied in our previous work. sarkadi2021arxiv  (B)
The bipolar nanopore’s device function is rectification. This is the system
studied here. The pores are cylindrical with radius $R=1$ nm and lengths
$H=2-24$ nm.

In the short and narrow pores we consider, the mechanism of rectification is
that the depletion zones of the coions are deeper in the OFF state of the
voltage than in the ON state. Rectification scales with
$\lambda_{\mathrm{D}}/R$ because $\lambda_{\mathrm{D}}/R$ tunes the degree of
overlap of the DLs, and, consequently, the deepness of the depletion zones.
fertig_jpcc_2019  The OFF voltage can make an already deep depletion zone even
deeper and can produce rectification. The importance of the
$\lambda_{\mathrm{D}}/R$ ratio was also emphasized in numerous earlier
studies. daiguji_nl_2005 ; vlassiouk_nl_2007 ; yan_nl_2009 ;
albrecht_chapter_2013 ; Abgrall_2008 ; bocquet_csr_2010 ; daiguji_csr_2010 ;
eijkel_csr_2010 ; zangle_csr_2010 ; madai_pccp_2018 ; dal_cengio_jcp_2019

While $\sigma$ was constant in our study for the scaling of rectification
fertig_jpcc_2019 , its effect was examined separately in a study,
fertig_pccp_2020  where a non-monotonic $\sigma$-dependence of rectification
was observed for electrolytes containing multivalent ions (2:2, 2:1, 3:1).
This behavior is caused by strong ionic correlations resulting in charge
inversion, overcharging, and anion leakage. These phenomena were reproduced by
a hybrid simulation method that unifies the advantages of a simple transport
equation and a particle simulation method, particularly, the Nernst-Planck
(NP) equation and the Local Equilibrium Monte Carlo (LEMC) method that is a
generalization of the Grand Canonical Monte Carlo (GCMC) method to non-
equilibrium situations. The methodology will be described in detail in section
II.

The non-monotonic behavior was not reproduced by the Poisson-Nernst-Planck
(PNP) theory that is a mean-field continuum theory using the Poisson-Boltzmann
(PB) theory with the NP equation. For 1:1 electrolytes, where ionic
correlations are weak, the curves were monotonic and the NP+LEMC and PNP
results agreed well.

In a recent paper, sarkadi2021arxiv  we considered a cylindrical nanopore with
a uniform surface charge on its wall (Fig. 1A). For this nanopore, the obvious
device function is selectivity that can be defined as
$S_{+}=I_{+}/(I_{+}+I_{-})$, where $I_{i}$ is the ionic current carried by
ionic species $i$. We simulated nanopores of different radii, lengths, and
surface charge values, and also varied the concentration of the electrolyte
and voltage. We restricted ourselves to 1:1 electrolytes (as we do in this
study) because we were curious when scaling works instead of when it does not
work.

We found that identifying a “universal” scaling parameter that includes all
the relevant parameters ($R$, $H$, $\sigma$, $c$, and $U$) is a real
challenge. A solid result that we reported is that the appropriate scaling
parameter for the infinitely long ($H/R\rightarrow\infty$) pore (nanotube
limit) is the Dukhin number defined as bazant_pre_2004 ; chu_pre_2006 ;
bocquet_chemsocrev_2010

$\mathrm{Du}=\frac{|\sigma|}{eRc},$ (4)

where $e$ is the elementary charge. Du was originally introduced by Bikerman
bikerman_1940  to characterize the ratio of the surface and volume
conductances considering electrokinetic phenomena. Later, Dukhin adopted the
idea (see Ref. dukhin_advcollsci_1993  and references therein) in studying
electrophoretic phenomena. His name was linked to the variable thanks to
Lyklema who suggested the Du notation to salute Dukhin lyklema_book_1995 ,
although the name ‘Bikerman number’ is also present in the literature. The
Dukhin number proved to be useful in nanopore studies by replacing the colloid
particle’s radius with the pore radius to characterize a geometrical
constraint in relation to an electrolyte property, the screening length.
bazant_pre_2004 ; chu_pre_2006 ; khair_jfm_2008 ; das_langmuir_2010 ;
bocquet_chemsocrev_2010 ; zangle_csr_2010 ; lee_nanolett_2012 ; yeh_ijc_2014 ;
ma_acssens_2017 ; xiong_scc_2019 ; poggioli_jpcb_2019 ; dalcengio_jcp_2019 ;
kavokine_annualrev_2020 ; noh_acsnano_2020

We proposed using an alternative formula for Du that, for a 1:1 electrolyte,
reads as

$\mathrm{Du}=\frac{|\sigma|l_{\mathrm{B}}^{*}\lambda_{\mathrm{D}}^{2}}{eR},$
(5)

where $l_{\mathrm{B}}^{*}=8\pi l_{\mathrm{B}}$ and
$l_{\mathrm{B}}=e^{2}/4\pi\epsilon_{0}\epsilon kT$ is the Bjerrum length. It
is based on the fact that the Debye length can be written as
$\lambda_{\mathrm{D}}^{2}=1/(8\pi l_{\mathrm{B}}c)=1/(l_{\mathrm{B}}^{*}c)$.
This form of Du makes it possible to relate quantities with the dimensions of
distance to each other. Also, we can use $\lambda_{\mathrm{MSA}}$ instead of
$\lambda_{\mathrm{D}}$ in Eq. 5 if we feel appropriate. A screening length
that depends on confinement levy_jcis_2020  is also possible.

For the nanotube limit ($H/R\rightarrow\infty$), the variables $H$ and $U$ are
absent; we used equilibrium calculations (PB and GCMC) to obtain selectivity.
For a finite nanopore, these variables become increasingly important because
the DLs that appear at the entrances of the pore near the membrane on the two
sides influence the ionic distributions inside the pore. In the nanohole limit
($H/R\rightarrow 0$), we found that a modified Dukhin number

$\mathrm{mDu}=\mathrm{Du}\frac{H}{\lambda_{\mathrm{D}}}=\frac{|\sigma|l_{\mathrm{B}}^{*}\lambda_{\mathrm{D}}H}{eR}$
(6)

is a much better scaling parameter.

Both the nanotube and nanohole limits are relevant experimentally. The
nanotube limit is obvious in cases, where a channel is fabricated in a thick
membrane as in PET nanopores Siwy_2004 . The nanohole limit is relevant in the
case of thin membranes, for example, graphene, MoS2, or WS2 garaj_n_2010 ;
garaj_pnas_2013 ; OHern_nl_2014 ; rollings_nc_2016 ; thiruraman_nl_2018 ;
thiruraman_acsnano_2020 .

If we plot $S_{+}$ as a function of $\mathrm{mDu}$ with a logarithmic scaling
on the abscissa, e.g., $S_{+}$ vs. $\lg(\mathrm{mDu})$, the resulting curve is
a sigmoid whose inflection point has a basic importance. First, the inflection
point separates state points for which surface conduction dominates from state
points where volume (bulk) conduction dominates. $S_{+}\sim 0.5$ corresponds
to the non-selective cases, when volume conduction dominates; both cations and
anions contribute to the current. $S_{+}\sim 1$ corresponds to the selective
cases, when surface conduction dominates; the DL extends over the pore in the
$r$ dimension and only the counterions contribute to the current. The
inflection point (where $S_{+}\approx 0.75$) offers itself to be a transition
point between these limiting cases.

When we talk about scaling, however, we are interested in the cases in between
where $S_{+}$ is intermediate. The inflection point is a useful asset because
if the inflection points of various scaling curves are at the same
$\lg(\mathrm{mDu})$ values, we can say that the curves fall onto each other.
Such a role of the inflection point was not found for bipolar nanopores.

The limiting behavior of the voltage dependence for $U\rightarrow 0$ is
different in selective and bipolar nanopores. In the selective case (uniform
$\sigma$, Fig. 1A), the $U\rightarrow 0$ limit provides a limiting value for
selectivity,

$\lim_{U\rightarrow 0}S_{+}(U)=S_{+}^{0},$ (7)

the slope conductance limit. In Ref. sarkadi2021arxiv  we reported a quadratic
dependence of the inflection point on voltage. Voltage, however, is absent in
Eq. 6. Including $U$ in the scaling parameter requires more data.

In the case of the rectifying case (bipolar $\sigma$, Fig. 1B), rectification
is based on the fact that the voltage has different effects in the ON and OFF
cases. Larger voltage produces larger rectification because larger voltage
produces larger differences between the concentration profiles in the ON and
OFF cases. Accordingly, rectification vanishes in the $U\rightarrow 0$ limit:

$\lim_{U\rightarrow 0}r(U)=1.$ (8)

We will show in this paper that $\sigma$ and $U$ have similar roles in the
bipolar nanopore: they both tune rectification. They do it with different
mechanisms, but from the point of view of scaling, with similar results.
Therefore, we redefine mDu for bipolar nanopores as

$\mathrm{mDu}=\frac{|\sigma|l_{\mathrm{B}}^{*}\lambda_{\mathrm{D}}H}{eR}\frac{U}{U_{0}},$
(9)

where $U_{0}$ is an arbitrary voltage to make mDu dimensionless ($U_{0}=200$
mV in this work). For $U=200$ mV (a value that we used in most of our
simulations), Eqs. 6 and 9 provide the same value.

The length of the nanopore, $H$, also influences the device behavior
differently in the selective and bipolar cases. Because bipolar nanopores are
necessarily finite due to their asymmetry, the nanotube limit
($H\rightarrow\infty$) does not make sense in this case. Indeed, the original
Dukhin number is not an appropriate scaling parameter for the bipolar
nanopore. In this paper, we examine whether mDu is any better for this
purpose.

The length of the pore influences the selective pore’s behavior in the middle
of the pore because the DLs formed at the entrances extend in the pore in the
axial dimension; the degree of this extension is characterized by
$\lambda_{\mathrm{D}}/H$. The behavior of the bipolar nanopore is also
sensitive to $H$. The depletion zones are formed at the junction of the two
differently charged regions (“n” and “p”) and they are deeper in the OFF state
if the nanopore is longer.

In this paper, we examine the applicability of a scaling parameter
($\mathrm{mDu}$, Eq. 9) that includes all relevant system parameters, $R$, $H$
$c$, $\sigma$, and $U$. Even if scaling is not perfect, showing our results in
terms of such a parameter is a compressed way of analyzing the results. This
is useful even if we express the device function in a subset of the parameter
space, e.g., changing only a few (one or two) of all the parameters, while
keeping the rest of them fixed.

The device function is rectification (Eq. 1) in this study. It proved to be
useful in our previous work sarkadi2021arxiv  that selectivity was defined as
a bounded function ($0.5<S_{+}<1$) that produced a sigmoid curve as a function
of $\lg(\mathrm{mDu})$. Following the practice of that paper, we pursue
sigmoid curves in this work as well, and introduce the following definition
for rectification:

$r^{\prime}=\frac{r-1}{r+1}=\frac{I^{\mathrm{ON}}-I^{\mathrm{OFF}}}{I^{\mathrm{ON}}+I^{\mathrm{OFF}}}.$
(10)

This function is $1$ for perfect rectification ($r\rightarrow\infty$) and $0$
for no rectification ($r=1$).

##  II Model and method

The model of the nanopore/membrane/baths system has been constructed of two
baths separated by a membrane and connected with a cylindrical pore through
the membrane. The results presented in this paper will be expressed as
functions of cylindrical coordinates ($z$ and $r$) due to the rotational
symmetry of the system around the axis of the pore. The simulation domain is a
cylinder with a $15-30$ nm length and $8-9$ nm radius (depending on $c$ and
$H$), while the pore’s length, $H$, is varied between $2$ and $24$ nm. The
membrane and the pore are confined by hard walls with a surface charge placed
on the wall of the pore (Fig. 1). The radius of the pore is constant in this
study ($R=1$ nm).

The pore has two regions: a positively charged region with surface charge
density $\sigma$ and a negatively charged region with surface charge density
$-\sigma$; the length of each is $H/2$. The surface charge density on the
cylinder is modeled as a collection of fractional point charges on a
$0.2\times 0.2$ nm grid ($0.3\times 0.3$ for longer pores above $H=15$ nm.

The electrolyte is modeled in the implicit solvent framework. This means that
the effect of the water molecules on the ions is taken into account implicitly
by response functions. One response function is the dielectric response
(screening) represented by the dielectric constant of the electrolyte. In our
study, it is constant system-wide and appears in the denominator of the
Coulomb potential. The interaction potential between ions can be expressed as
the sum of the hard-sphere potential and the Coulomb potential:

$u_{ij}(r)=\left\\{\begin{array}[]{ll}\infty&\quad\mathrm{if}\quad
r<R_{i}+R_{j}\\\
\dfrac{1}{4\pi\epsilon_{0}\epsilon}\dfrac{z_{i}z_{j}e^{2}}{r}&\quad\mathrm{if}\quad
r\geq R_{i}+R_{j}\\\ \end{array}\right.$ (11)

where $R_{i}$ and $R_{j}$ are the radii ($R_{+}=R_{-}=0.15$ nm in this study),
$z_{i}$ and $z_{j}$ are the valences of the different ionic species, and $r$
is the distance between ions $i$ and $j$.

The other response function is dynamic in nature; the water molecules collide
with ions and impede their diffusion. We take that effect into account by a
diffusion coefficient function, $D_{i}(\mathbf{r})$, in the NP equation that
reads as

$\mathbf{j}_{i}(\mathbf{r})=-\frac{1}{kT}D_{i}(\mathbf{r})c_{i}(\mathbf{r})\nabla\mu_{i}(\mathbf{r}),$
(12)

where $\mathbf{j}_{i}(\mathbf{r})$, $D_{i}(\mathbf{r})$, $c_{i}(\mathbf{r})$,
and $\mu_{i}(\mathbf{r})$ are the flux density, the diffusion coefficient
profile, the concentration profile, and the electrochemical potential profile
of ionic species $i$, respectively.

For the diffusion coefficient profile, $D_{i}(\mathbf{r})$, we use a piecewise
constant function, where the value in the baths is $1.334\times 10^{-9}$ m2s-1
for both ionic species, while it is the tenth of that inside the pore,
$D_{i}^{\mathrm{pore}}$, as in our earlier works. matejczyk_jcp_2017 ;
madai_jcp_2017 ; madai_pccp_2018 ; fertig_jpcc_2019 ; fertig_pccp_2020 ;
sarkadi2021arxiv  These particular choices do not qualitatively affect our
conclusions.

The electrochemical potential profile can be expressed as

$\mu_{i}(\mathbf{r})=\mu_{i}^{0}+kT\ln
c_{i}(\mathbf{r})+\mu_{i}^{\mathrm{BMF}}(\mathbf{r})+z_{i}e\Phi(\mathbf{r}),$
(13)

where $\mu_{i}^{0}$ is a temperature dependent reference term (it carries the
information about the unit of $c_{i}(\mathbf{r})$; strictly speaking, $\ln
c_{i}$ should be written as $\ln(c_{i}/c_{0})$, where $c_{0}$ is a reference
concentration also contained by $\mu_{i}^{0}$), $\Phi(\mathbf{r})$ is the mean
(ensemble or time averaged) electrical potential, and
$\mu_{i}^{\mathrm{BMF}}(\mathbf{r})$ is an excess term beyond the mean field
(BMF) approximation.

The chemical and electrical contributions,
$\mu_{i}^{\mathrm{CHEM}}(\mathbf{r})$ and $\mu_{i}^{\mathrm{EL}}(\mathbf{r})$,
can be defined formally even if they are difficult to separate in experiments.
The chemical term is

$\mu_{i}^{\mathrm{CHEM}}(\mathbf{r})=\mu_{i}^{0}+kT\ln
c_{i}(\mathbf{r})+\mu_{i}^{\mathrm{BMF}}(\mathbf{r}),$ (14)

while the electrical term is

$\mu_{i}^{\mathrm{EL}}(\mathbf{r})=z_{i}e\Phi(\mathbf{r}).$ (15)

The BMF and EL terms together constitute the excess term.

The ion transport is steady-state; the profiles in the above equation do not
depend on time. The transport is maintained by constant (time-independent)
boundary conditions for the concentration and the electrical potential on the
two sides of the membrane on the two half cylinders confining the simulation
cell from left and right.

Fixing the concentration ($c_{i}^{\mathrm{L}}$ and $c_{i}^{\mathrm{R}}$) and
the electrical potential ($\Phi^{\mathrm{L}}$ and $\Phi^{\mathrm{R}}$) on the
left (L) and right (R) boundaries means fixing the electrochemical potentials
that can be written as

$\mu_{i}^{\mathrm{L}}=\mu_{i}^{0}+kT\ln
c_{i}^{\mathrm{L}}+\mu_{i}^{\mathrm{BMF,L}}+z_{i}e\Phi^{\mathrm{L}}$ (16)

and

$\mu_{i}^{\mathrm{R}}=\mu_{i}^{0}+kT\ln
c_{i}^{\mathrm{R}}+\mu_{i}^{\mathrm{BMF,R}}+z_{i}e\Phi^{\mathrm{R}},$ (17)

where $\mu_{i}^{\mathrm{BMF,L}}$ and $\mu_{i}^{\mathrm{BMF,R}}$ are the excess
chemical potentials in the absence of an external field determined by the
Adaptive GCMC method of Malasics et al. malasics_jcp_2010 . The
$z_{i}e\Phi^{\mathrm{L}}$ and $z_{i}e\Phi^{\mathrm{R}}$ terms are the
interactions with the applied electrical potentials in the two baths.
Prescribing $\Phi^{\mathrm{L}}$ and $\Phi^{\mathrm{R}}$ on the system’s
boundary means that we use an electrostatic Dirichlet boundary condition.
Voltage is defined as $U=\Phi^{\mathrm{R}}-\Phi^{\mathrm{L}}$ in our study.

To make use of the NP equation, we need a relation between the concentration
profile, $c_{i}(\mathbf{r})$, and the electrochemical potential profile,
$\mu_{i}(\mathbf{r})$. Most of our results reported in this work have been
obtained by the LEMC technique developed by Boda and Gillespie boda_jctc_2012
. It is a variant of the GCMC method that is generally applied in global
equilibrium ($\mu_{i}$ is constant in space), while only local equilibrium is
assumed in LEMC in subvolumes of the simulation cell. Inside these subvolumes,
$\mu_{i}(\mathbf{r})$ is constant, but it varies between subvolumes. The
resulting varying $\mu_{i}(\mathbf{r})$ profile provides a driving force for
the ionic transport. Coupled with NP (denoted by NP+LEMC), an iterative
procedure is obtained that provides a solution where $c_{i}(\mathbf{r})$ is
statistical-mechanically consistent with $\mu_{i}(\mathbf{r})$ and also they
together produce a flux via the NP equation that satisfies conservation of
mass ($\nabla\cdot\mathbf{j}(\mathbf{r})=0$).

The LEMC simulations use the interaction potentials given in Eq. 11 also
called charged hard spheres or the primitive model of electrolytes. Overlap of
hard-sphere ions with each other and with the hard walls in the system (wall
of the pore, the membrane, and the system’s outer boundary) is forbidden.

The other method uses the PB theory to relate $c_{i}(\mathbf{r})$ to
$\mu_{i}(\mathbf{r})$. In this case, explicit particles are not present in the
calculation; ions are described by the $c_{i}(\mathbf{r})$ functions. The
hard-sphere cores are absent, so the ions are treated as point charges. When
we couple the PB theory to the NP equation, the PNP theory is obtained that is
a mean-field continuum theory applied extensively to bipolar nanopores
daiguji_nl_2005 ; constantin_pre_2007 ; vlassiouk_nl_2007 ; karnik_nl_2007 ;
vlassiouk_acsnanno_2008 ; kalman_am_2008 ; yan_nl_2009 ; nguyen_nt_2010 ;
szymczyk_jpcb_2010 ; singh_jap_2011 ; singh_jpcb_2011 ; singh_apl_2011 ;
oeffelen_oneplos_2015 ; tajparast_bba_2015 . Note that the NP equation has
also been coupled to a density functional theory by Gillespie.
gillespie_jpcm_2002 ; gillespie_pre_2003 ; gillespie_jpcb_2005 ;
gillespie_bj_2008_energetics

Being mean field means that the ions interact with other ions only via the
mean electrical potential produced by them. Formally, this means that the BMF
term in Eq. 13 is zero. The resulting electrochemical potential is that of an
ideal solution. The PNP theory is quite standard in nanopore studies and is
described elsewhere in detail. constantin_pre_2007 ; burger_nonlin_2012 ;
eisenberg_pnp_2019  Our procedure has been detailed in Ref. matejczyk_jcp_2017
.

Because the NP+LEMC method is less standard, we provide a brief description in
Appendix A, while more details can be found in our earlier papers.
boda_jctc_2012 ; boda_arcc_2014 ; boda_jml_2014 ; fertig_hjic_2017

##  III Results and Discussion

###  III.1 Mechanism of rectification

Let us introduce our discussion with axial concentration profiles to pin down
the mechanism by which a bipolar nanopore rectifies ionic current. The cross-
section-averaged concentration profile is defined as

$c_{i}(z)=\dfrac{1}{(R(z)-R_{i})^{2}\pi}\int\limits_{0}^{R(z)-R_{i}}c_{i}(z,r)2\pi
r\mathrm{d}r,$ (18)

where $R(z)$ is the radius of the pore for $-H/2<z<H/2$ and the radius of the
simulation cell outside the pore. Fig. 2 shows that changing the sign of the
voltage form ON to OFF results in reduced concentrations inside the pore.
Regions where concentrations are reduced are generally termed depletion zones
in the semiconductor, nanopore, and ion channel literature. Depletion zones of
coions are formed as results of the surface charges that repulse the coions.
In the “p” zone cations have depletion zones (they are the coions), while in
the “n” zone the anions have depletion zones (they are the coions). The basis
of the mechanism of rectification is that the OFF voltage makes the depletion
zones deeper compared to the ON state.

Figure 2: Axial concentration profiles for a fixed concentration, $c=0.1$ M,
pore radius, $R=1$ nm, and pore length, $H=6$ nm. The surface charge,
$\sigma$, and voltage, $U$, are chosen in the two panels in such a way that
their product is constant; $\sigma=0.1$ $e$/nm2 and $U=100$ mV in the top
panel, while $\sigma=0.5$ $e$/nm2 and $U=20$ mV in the bottom panel. Blue and
red lines refer to cation and anion profiles, respectively. Red and light red
(blue and light blue) profiles refer to OFF and ON states, respectively. In
general, blue and red colors refer to cations and anions, respectively, here
and later figures.

This indicates that in the following discussions we need to pay extra
attention to the depletion zones in the OFF state because tuning certain
parameters ($\sigma$, $c$, or $H$, for example) changes rectification via
modifying the depletion zones. We need to state right at the beginning,
however, that influencing the ON-state profiles is also a legal way to tune
rectification. Indeed, as we will show, increasing voltage have a stronger
effect on increasing the ON-state concentration profiles than the OFF-state
concentration profiles.

Fig. 2 shows that the $c_{i}(z)$ profiles behave differently in the top and
the bottom panels, still, rectification is similar in the two cases. The two
panels refer to the same $\mathrm{mDu}$, but different pairs of $\sigma$ and
$U$. The agreement of rectifications indicates that mDu may be an appropriate
scaling parameter, but we need to examine the scaling curves to show this.

###  III.2 Scaling curves

Next, let us turn our attention to the device-level behavior, namely, let us
show the device function (rectification) as a function of the chosen scaling
parameter, mDu. As it was stated before, the scaling parameter does not need
to be a perfect one. Our main premise here is that even an imperfect scaling
parameter is useful in explaining the behavior of the device and in analyzing
the effects of the various parameters.

Since we have many system parameters ($R$, $H$, $c$, $U$, and $\sigma$ and
there would be even more for multivalent electrolytes such as $z_{+}$ and
$z_{-}$), we found it unfeasible to find the Holy Grail of scaling parameters,
a parameter as a function of which all the rectification curves fall onto each
other perfectly. Such a “perfect” scaling was found for the infinitely long
selective pore in our previous work, sarkadi2021arxiv  but deviations were
found when we decreased the length of the pore. An interplay between radial
and axial effects took place that also occurs in the bipolar pore considered
here.

Figure 3: Rectification, $r^{\prime}$, as a function of
$\mathrm{mDu}=|\sigma|l_{\mathrm{B}}^{*}\lambda_{\mathrm{D}}HU/RU_{0}$. In
panel (A) pore radius and pore length are kept fixed at $R=1$ nm and $H=6$ nm,
while concentrations are $c=0.1$ M (top panel) and $c=1$ M (bottom panel). The
various symbols/curves have been obtained by either changing $\sigma$ at fixed
$U$ or by changing $U$ at fixed $\sigma$ as indicated in the legends. The
inset in the bottom panel for $c=1$ M was obtained by using
$\lambda_{\mathrm{MSA}}=0.415$ nm instead of $\lambda_{\mathrm{D}}=0.304$ nm
in the case of the NP+LEMC data. The curve (open blue symbols) was obtained by
merging all the NP+LEMC data and fitting a sigmoid. In panel (B) pore radius
and voltage are kept fixed at $R=1$ nm and $U=200$ mV, while pore lengths are
$H=6$ (top panel) and $H=24$ nm (bottom panel). The various symbols/curves
refer to different concentrations and have been obtained by changing $\sigma$.

So our purpose here is to show that changing $\sigma$, $U$, and $H$ have
similar effects on the behavior of the nanopore. They are all in the numerator
of mDu, so if scaling is valid, doubling any of these parameters will double
rectification. Although such an accurate scaling is not present, we can put
forward a weaker statement that increasing any of these parameters will
increase rectification.

Concentration, for a given electrolyte, determines the screening length, so it
determines the degree of overlap of the DLs formed at the nanopore’s wall.
That way, the relation of the screening length and pore radius,
$\lambda_{\mathrm{D}}/R$ , determines rectification if everything else is
constant as it was shown in our earlier publication. fertig_jpcc_2019  Because
this question was analyzed in detail in that paper, we pay less attention to
it here; instead, we fix the pore radius at $R=1$ nm and study the device
behavior at various fixed concentrations, $c=0.01$, $0.1$, and $1$ M, with
$\sigma$, $H$, and $U$ being the main variables.

Fig. 3 shows scaling curves under different circumstances. Fig. 3A shows that
we obtain similar behavior either by changing $\sigma$ or $U$. As we will show
through the behavior of the ionic distributions later in this work, increasing
any of these two parameters increases rectification. This figure, furthermore,
implies linearity: it does not matter if we change $\sigma$ or $U$ to the same
degree, we get the same effect although via different mechanisms as implied by
Fig. 2. Fig. 3A suggests that mDu is a fare scaling parameter as far as the
relation of $\sigma$ and $U$ is concerned. A close examination of the
differences between the two panels, however, reveals that we have different
scaling behaviors for $c=0.1$ M and $1$ M.

This difference is better shown in Fig. 3B that shows scaling curves for
different concentrations inside one panel for a fixed voltage, $U=200$ mV.
Indeed, the different curves for $c=0.01$, $0.1$, and $1$ M in the top panel
for $H=6$ nm are relatively far apart. The $c=0.1$ and $1$ M curves are
closer, while the $c=0.01$ M curve deviates more. This is the result of the
fact that the screening lengths at the studied concentrations are comparable
in size to the length and radius of the pore; the Debye lengths are
$\lambda_{\mathrm{D}}=3.042$, $0.962$ , and $0.304$ nm for $c=0.01$, $0.1$,
and $1$ M, respectively. Such wide DLs couple the axial and radial behaviors
of the ionic distributions.

If the insight about the effect of the axial interaction between the “p” and
“n” regions in short pores is right, we should observe better scaling in
longer pores. Indeed, the bottom panel of Fig. 3B for $H=24$ nm shows that the
curves for $c=0.1$ and $1$ M scale appropriately. The pore is long enough
compared to the screening lengths of the $c=0.1$ and $1$ M electrolytes. The
$c=0.01$ M curve is still off, but to a lesser extent than in the shorter
pore.

The main message of Fig. 3B is that linear scaling as a function of $H$
($\mathrm{mDu}\sim H$) works quite well, though it is far from being perfect.
On the basis of the available results, we may predict that scaling would
improve with increasing $H$, but at this point of the research, this is only a
hypothesis. A further complication of the picture is that
$\lambda_{\mathrm{D}}$ is not necessarily a good screening length in
confinement as implied by the study of Levy et al. levy_jcis_2020 .

For simplicity, we use the Debye length for $\lambda$ in this study for both
NP+LEMC and PNP. In the inset of the bottom panel of Fig. 3A, however, we also
show the results for the choice $\lambda=\lambda_{\mathrm{MSA}}$ in the case
of NP+LEMC for the large concentration of $c=1$ M where $\lambda_{\mathrm{D}}$
and $\lambda_{\mathrm{MSA}}$ are quite different ($0.304$ and $0.415$,
respectively). As it was already pointed out in Ref. fertig_jpcc_2019 , using
$\lambda_{\mathrm{MSA}}$ for NP+LEMC and $\lambda=\lambda_{\mathrm{D}}$ for
PNP produces better agreement between the two models. This is also seen here:
the NP+LEMC and PNP curves agree if we use the appropriate screening lengths
for the two methods.

Fundamentally, $\lambda_{\mathrm{MSA}}$ suits LEMC better because MSA also
takes BMF correlations into account even if in an approximate manner. We use
$\lambda_{\mathrm{D}}$ in this work because it is easier to calculate, and
using it also makes our point. The equations for $\lambda_{\mathrm{MSA}}$ are
found in earlier publications. blum_mp_1975 ; blum_jcp_1977 ; nonner_bj_2000 ;
fertig_jpcc_2019

The PNP curve in Fig. 3A has also been computed by fixing $\sigma$ and
scanning $U$. The resulting curve practically coincides with the one shown in
the figure for fixed $U$ and scanned $\sigma$. This implies that scaling works
better in the mean-field approximation, while ionic correlations may cause
deviations in the NP+LEMC data. This is more apparent for multivalent
electrolytes where ionic correlations are strong (data not shown).
fertig_pccp_2020

Here, we are restricted to the question of why changing $\sigma$, $H$, or $U$
has similar effects on rectification for a given value of $c$ and $R$. To
answer this question, we need to examine more detailed results provided by the
simulations, e.g., concentration, electrical potential, and electrochemical
potential profiles.

###  III.3 Self consistent NP+LEMC system

Cation and anion concentration profiles are separated both in the radial and
axial dimensions. In the radial dimension, it is the surface charge that
separates the cation and anion profiles and produces the radial DL. The
overlap of this DL in the centerline of the pore produces the depletion zone
that is tuned by the applied field. Here, we do not show radial profiles; they
have been discussed in detail in our previous publications, for example, in
Figs. 2 and 5 of Ref. fertig_jpcc_2019  and Figs. 5, 6, and 9 of Ref.
fertig_pccp_2020 .

In the axial dimension, it is the applied field and the charge imbalance along
the $z$ axis that modifies the ionic distributions in a given region resulting
in axial DLs. If the screening length is large, the radial and axial DLs
overlap and their structures are mutually correlated.

An accurate screening, however, seems to be based on a decoupling of the
radial and axial effects. If the pore is short compared to
$\lambda_{\mathrm{D}}$ ($\lambda_{\mathrm{D}}/H$ is large), the “p” region has
a strong effect on the ionic distributions in the “n” region, and vice versa.
The neighboring region influences the radial behavior of the concentration
profiles, the degree of DL overlap, and, thus, the formation of depletion
zones. The relation of the pore shape and the screening length ($R$ vs. $H$
vs. $\lambda_{\mathrm{D}}$) will be elaborated on in a subsequent publication.

Fig. 4 explains how $c_{i}(z)$, $\Phi(z)$, and $\mu_{i}(z)$ cooperate to
provide a self consistent solution of the NP+LEMC system (and the PNP system
as well). The top and bottom rows show the ON and OFF states ($U=200$ and
$-200$ mV), respectively. The left panels show the cross-section-averaged (Eq.
18) mean electrical potential profiles (blue) and their two components:

$\Phi(z)=\Phi^{\mathrm{APPL}}(z)+\Phi^{\mathrm{ION}}(z),$ (19)

where $\Phi^{\mathrm{APPL}}(z)$ is the applied electrical potential (red)
obtained by solving Laplace’s equation (Eq. 31) with the prescribed Dirichlet
boundary conditions ($\Phi=0$ on the left end, while $\Phi=U$ on the right
end), while $\Phi^{\mathrm{ION}}(z)$ is the mean electrical potential produced
by the ions and the surface charges inside the pore (black).
$\Phi^{\mathrm{ION}}(z)$ has been computed by “on the fly” sampling in the
LEMC simulations, but solving Poisson’s equation leads to the same result. A
potential profile with positive and negative peaks in the “p” and “n” regions
is produced by the procedure.

Figure 4: Cross-section-averaged axial profiles for the mean electrical
potential ($\Phi(z)$, left) and the electrochemical potential of the cation
($\mu_{+}(z)$, right). Top and bottom rows refer to the ON and OFF voltages
($U=200$ and $-200$ mV), respectively. The mean electrical potential (blue) is
split into the terms of the applied potential (red) and the potential produced
by the ions and the pore charges (black) as defined in Eq. 19. The
electrochemical potential (green) is split into the chemical (black) and
electrical (blue) terms as defined in Eqs. 13-15. The blue curves on the right
are equal to the blue curves on the left multiplied by $e/kT$. The
$\mu_{i}^{\mathrm{CHEM}}(z)$ and $\mu_{i}(z)$ terms are shifted with
$\mu_{i}^{\mathrm{L}}$. The open green symbols represent the derivative of the
electrochemical potential, $\mathrm{d}\mu_{+}(z)/\mathrm{d}z$. The parameters
are $R=1$ nm, $H=12$ nm, $c=0.1$ M, and $\sigma=1$ $e$/nm2.

The total mean electrical potential, $\Phi(z)$, multiplied by $z_{i}e$ is
defined as the electrical part, $\mu_{i}^{\mathrm{EL}}(z)$, of the
electrochemical potential, $\mu_{i}(z)$ (see Eqs. 13 and 15), while the
difference between them is defined as the chemical term,
$\mu_{i}^{\mathrm{CHEM}}(z)$ (see Eq. 14). The $\mu_{i}^{\mathrm{CHEM}}(z)$
and $\mu_{i}^{\mathrm{EL}}(z)$ contributions are plotted by black and blue
lines, while $\mu_{i}(z)$ is plotted by green lines in the right panels of
Fig. 4. Because the gradient of $\mu_{i}$ is the driving force of the
transport, we show the derivative of $\mu_{i}(z)$ with open symbols in the
bottom-right panel to illustrate that $\mathrm{d}\mu_{i}(z)/\mathrm{d}z$ has a
minimum where $\mu_{i}^{\mathrm{CHEM}}(z)$ profile has a minimum.

The minimum of the $\mu_{i}^{\mathrm{CHEM}}(z)$ profile, on the other hand,
corresponds to a minimum in the $c_{i}(z)$ profile (a depletion zone) because
$\mu_{i}^{\mathrm{CHEM}}(z)\sim\ln c_{i}(z)$. This is because the BMF term in
Eq. 14 is small for the 1:1 system studied here.

The picture is complicated, but as a rule of thumb, we can keep in mind that
$\mu_{i}(z)$ (green curve) must be monotonic, so $\mu_{i}^{\mathrm{CHEM}}(z)$
must balance $\mu_{i}^{\mathrm{EL}}(z)$ accordingly. The result is that
depletion zones are formed in the OFF state at the junction (the interface of
the “p” and “n” region).

In summary, in the self consistent calculation of the NP+LEMC system, ionic
distribution profiles, $c_{i}(z,r)$, must be obtained that produce a mean
electrical potential, $\Phi(z,r)$, that, together with
$\mu_{i}^{\mathrm{CHEM}}(z,r)$, produce a $\mu_{i}(z,r)$ profile that,
together with $c_{i}(z,r)$, produce a flux density, $\mathbf{j}_{i}(z,r)$,
that satisfies conservation of mass, $\nabla\cdot\mathbf{j}_{i}(z,r)=0$.

###  III.4 Slope conductance approach

Plotting concentration profiles is not necessarily the best way to get the
appropriate insight in the case of nanopores, especially, if the nanopore’s
behavior is driven by depletion zones as it was illustrated by Fig. 2.

It has already been realized and discussed in numerous several publications of
ours gillespie_bj_2008_ca ; gillespie_bj_2008_nanopore ; he_jacs_2009 ;
boda_jgp_2009 ; malasics_bba_2010_trivalent ; valisko_jcp_2019 ;
fertig_pccp_2020 ; boda_entropy_2020  that the conductance of the pore is
better described by the reciprocal of the concentration than the concentration
itself as soon as the resistance of the pore is determined by the low-
resistance elements of the consecutive segments of the pore along the $z$-axis
that are considered as resistors connected in series. This approach works if
there are so deep depletion zones in the pore that the importance of the peaks
are dwarfed compared to them.

This is the case, for example, in ion channels that are naturally narrow pores
where coions are excluded from regions dominated by charged amino acids as in
the case of calcium channels gillespie_bj_2008_ca ; boda_jgp_2009 ;
malasics_bba_2010_trivalent  whose selectivity filter is lined by E or D amino
acids. Coions, however, can also be excluded from wide pores if the screening
length of the electrolyte is large enough compared to the radius of the pore.
madai_pccp_2018 ; fertig_jpcc_2019 ; valisko_jcp_2019 ; fertig_pccp_2020 ;
boda_entropy_2020 ; sarkadi2021arxiv ; levine_jcis_1975 ; balme_sr_2015 ;
uematsu_jpcb_2018 ; green_jcp_2021

An elegant quantification of this idea is the slope conductance approach
gillespie_bj_2008_ca ; gillespie_bj_2008_nanopore ; boda_jgp_2009  in which we
assume that the chemical potential is constant in the radial dimension
($\mu_{i}(z,r)\approx\mu_{i}(z)$) inside the pore. The current carried by an
ionic species $i$, $I_{i}$, is computed from the cross-sectional integral of
$\mathbf{j}_{i}(z,r)$, inside the pore. By substituting $j_{i}(z,r)$ from the
NP equation (Eq. 12) we obtain that

$\displaystyle I_{i}$ $\displaystyle\approx$
$\displaystyle-\dfrac{z_{i}eD_{i}^{\mathrm{pore}}}{kT}\dfrac{\mathrm{d}\mu_{i}(z)}{\mathrm{d}z}\int\limits_{0}^{R-R_{i}}c_{i}(z,r)2\pi
r\mathrm{d}r=$ (20) $\displaystyle=$
$\displaystyle-\dfrac{z_{i}eD_{i}^{\mathrm{pore}}}{kT}\dfrac{\mathrm{d}\mu_{i}(z)}{\mathrm{d}z}A_{i}c_{i}(z),$

where $A_{i}=(R-R_{i})^{2}\pi$ is the effective cross section (where
$c_{i}(z,r)$ is nonzero for $-H/2<z<H/2$) and $c_{i}(z)$ is the cross-section-
averaged axial concentration profile defined in Eq. 18. If we divide by
$c_{i}(z)$ and integrate, Eq. 20 can be expressed as

$g_{i}=\dfrac{I_{i}}{U}=-\dfrac{z_{i}^{2}e^{2}A_{i}D_{i}^{\mathrm{pore}}}{kT}\left(\int\limits_{-H/2}^{H/2}\dfrac{\mathrm{d}z}{c_{i}(z)}\right)^{-1},$
(21)

where $g_{i}$ is the conductance for ionic species $i$ and

$U=\frac{1}{z_{i}e}\int_{-H/2}^{H/2}\mathrm{d}\mu_{i}(z)=\frac{1}{z_{i}e}\left[\mu_{i}(H/2)-\mu_{i}(-H/2)\right]\approx\Phi(H/2)-\Phi(-H/2)$
(22)

is the potential drop across the pore. Eq. 22 assumed that the concentration
is the same on the two sides of the membrane, so only the
$\mu_{i}^{\mathrm{EL}}(z)=z_{i}e\Phi(z)$ term of $\mu_{i}$ is different at
$z=-H/2$ and $H/2$.

The conductance, therefore, is associated with the reciprocal of the integral
of $c_{i}^{-1}(z)$ along the pore. If $c_{i}(z)$ is very small,
$c_{i}^{-1}(z)$ is very large, its integral is large, and the conductance is
small (resistance large).

Because the radii and diffusion constants of the cations and anions are the
same in this study, their currents are similar, and we can write for the
rectification that

$r\sim\frac{I^{\mathrm{ON}}_{i}}{I^{\mathrm{OFF}}_{i}}\sim\dfrac{\left[\int\left(c_{i}^{\mathrm{ON}}(z)\right)^{-1}dz\right]^{-1}}{\left[\int\left(c_{i}^{\mathrm{OFF}}(z)\right)^{-1}dz\right]^{-1}}=\dfrac{\int\left(c_{i}^{\mathrm{OFF}}(z)\right)^{-1}dz}{\int\left(c_{i}^{\mathrm{ON}}(z)\right)^{-1}dz},$
(23)

that is, rectification is proportional to the ratio of the integrals of the
reciprocal concentration profiles in the OFF and ON states for either $i$.
This is shown in Fig. 5 for the case already depicted in Fig. 2.

Figure 5: Axial reciprocal concentration profiles are shown for the same state
point shown in Fig. 2 in the top and the middle panels. The surface charge and
the voltage $\sigma=0.1$ $e$/nm2 and $U=100$ mV in the top panel, while
$\sigma=0.5$ $e$/nm2 and $U=20$ mV in the bottom panel (their product is
fixed). Blue and red lines refer to cation and anion profiles, respectively.
Red and light red (blue and light blue) profiles refer to OFF and ON states,
respectively. The bottom panel shows the ratio of the ON and OFF concentration
profiles.

###  III.5 Concentration profiles and reciprocal concentration profiles

#### The effect of voltage and surface charge

In Fig. 5, peaks in the OFF profile correspond to depletion zones, so we can
visualize depletion zones better this way. The above analysis and this figure
explain why Fig. 2 alone is misleading if we want to understand the
mechanisms. There, the relation between the ON and OFF curves’ behaviors is so
different in the top and bottom panels that it is hard to deduct from that
figure why the pore rectifies similarly in the two cases.

Figure 6: Axial reciprocal concentration profiles for fixed pore radius, $R=1$
nm, and voltage, $U=200$ mV. The pore length is $H=6$ nm on the left hand
side, while it is $H=24$ nm on the right hand side. The surface charge,
$\sigma$, and concentration, $c$, are chosen such a way in the different
panels that the mechanism of scaling is illustrated through the areas under
the curves (see Eq. 21). In the middle panels, the base point is $\sigma=0.1$
$e$/nm2 and $c=0.1$ M. In the top-left panel, $H$ is decreased to $6$ nm,
while $\sigma$ is increased to $0.4$ $e$/nm2 so that mDu is the same as in the
middle-right panel. In the top-right panel, $H$ is increased to $24$ nm, while
$\sigma$ is decreased to $0.025$ $e$/nm2 so that mDu is the same as in the
middle-left panel. In the bottom panels, $c$ is increased to $1$ M, while
$\sigma$ is increased to $0.3162$ $e$/nm2 so that mDu is the same as in the
middle panels ($\sigma\lambda_{\mathrm{D}}$ is the same). Blue and red lines
refer to cation and anion profiles, respectively. Red and light red (blue and
light blue) profiles refer to OFF and ON states, respectively.

The explanation is that Fig. 2 overemphasizes the peaks, while it is the
depletion zones that are relevant for rectification. Looking at Fig. 5, the
contradiction can be resolved easily. The ratio of the integrals of the OFF
and the ON profiles can be similar even if the actual behaviors of these
profiles are quite different. Indeed, if we plot the ratio of the ON and OFF
concentration profiles (bottom panel), it is apparent that they are similar in
the two cases. Note that one must be careful with the
$c_{i}^{\mathrm{ON}}(z)/c_{i}^{\mathrm{OFF}}(z)$ profiles because generally
the integrals of these profiles do not have anything to do with the
rectification. Despite this restriction, the
$c_{i}^{\mathrm{ON}}(z)/c_{i}^{\mathrm{OFF}}(z)$ function is useful in drawing
conclusions for the general behavior of the device function. If the
$c_{i}^{\mathrm{ON}}(z)/c_{i}^{\mathrm{OFF}}(z)$ profile is similar in two
cases, rectification may also be similar.

This is what we can see in Fig. 5 by comparing the top and the middle panels.
The top panel refers to $\sigma=0.1$ $e$/nm2 and $U=100$ mV; a relatively
small surface charge produces smaller peaks and less deep depletion zones,
namely, counter- and coions are separated less. A larger voltage is needed to
attain a separation of the ON and OFF profiles so that we obtain a given
rectification $r^{\prime}=0.27$. The middle panel refers to the same $\sigma
U$ product, but for a larger surface charge, $\sigma=0.5$ $e$/nm2, and a
smaller voltage, $U=20$ mV. In this case, we obtain a similar rectification
$r^{\prime}=0.23$ by a larger $\sigma$ separating cation and anion profiles
more, but, in turn, a smaller voltage separating the ON and OFF profiles. The
bottom panel shows the $c_{i}^{\mathrm{ON}}(z)/c_{i}^{\mathrm{OFF}}(z)$
profiles.

The bottom line is that $\sigma$ and $U$ do different things, but with similar
effects. The surface charge separates cation and anion profiles; in each
region there is always a peak of counterions and a depletion zone of coions.
The voltage then modulates these profiles, but it modulates them in different
ways at its ON and OFF signs. The OFF voltage makes the depletion zones
deeper, while the ON voltage makes the peaks higher. Increasing either
$\sigma$ or $U$ increases the difference between the ON and OFF profiles, but
via different mechanisms.

Figure 7: Axial reciprocal concentration profiles (middle and top rows) and
the $c_{i}^{\mathrm{ON}}(z)/c_{i}^{\mathrm{OFF}}(z)$ profiles (bottom row). In
the respective columns only one parameter is changed: (A) surface charge, (B)
pore length, (C) voltage. The rest of the parameters are kept fixed and are
$R=1$ nm, $c=0.1$ M, $H=6$ nm, $U=200$ mV, and $\sigma=1$ $e$/nm2. Blue and
red lines refer to cation and anion profiles, respectively. The (C) column
shows the $c_{i}(z)$ profiles, while columns (A) and (B) shows the
$c_{i}^{-1}(z)$ profiles.

#### The effect of pore length and surface charge

Now let us turn our attention to the interplay of $\sigma$ and $H$. The
product of these two parameters determines the net pore charge in a half
region, $\pm\sigma H\pi R$, if $R$ is constant. It is plausible that more net
charge repels more coions from the pore thus producing a deeper depletion
zone.

Fig. 6 shows reciprocal concentration profiles to illuminate this. The middle
row shows our two base points, where $\sigma=0.1$ $e$/nm2 and $c=0.1$ M
(voltage is kept fixed in these calculations at $U=200$ mV). The left and
right panels show the results for $H=6$ and $24$ nm, respectively. The $H=24$
nm results refer to $4$ times larger $\mathrm{mDu}$, and, thus, a larger
rectification according to scaling in Fig. 3. Indeed, this resonates with the
profiles; the ratio of the areas under the OFF profiles and the ON profiles is
much larger for $H=24$ nm than for $H=6$ nm.

So, what can we do if we want to get a rectification that is similar to that
for the $\sigma=0.1$ $e$/nm2, $c=0.1$ M, $H=24$ nm case (middle right panel)?
One choice is to decrease $H$ to $6$ nm, but increase $\sigma$ to $0.4$
$e$/nm2 (top left panel). Indeed, the relation between the OFF and ON curves
is similar in this panel and the middle right panel (follow the green arrow).
The other choice is to increase $\sigma$ and decrease $\lambda_{\mathrm{D}}$
so that $\sigma\lambda_{\mathrm{D}}$ is the same value. The bottom right panel
shows the data for $c=1$ M that corresponds to a Debye length that is
$\sqrt{10}$ times smaller than in the $c=0.1$ case
($\lambda_{\mathrm{D}}=0.304$ and $0.962$ nm, respectively). The corresponding
surface charge, therefore, is $\sqrt{10}$ times larger, namely,
$\sigma=0.3162$ $e$/nm2. Again, the relation of the OFF and ON profiles are
similar to that in the panel above (follow the green arrow).

Inversely, what can we do if we want to get a rectification that is similar to
that for the $\sigma=0.1$ $e$/nm2, $c=0.1$ M, $H=6$ nm case (middle left
panel)? One choice is to increase $H$ to $24$ nm, but decrease $\sigma$ to
$0.025$ $e$/nm2 (top right panel, follow the green arrow). The other choice is
the same as in the above paragraph: to increase $\sigma$ and decrease
$\lambda_{\mathrm{D}}$ so that $\sigma\lambda_{\mathrm{D}}$ is the same value
(bottom left panel, follow the green arrow). The conclusions regarding the
relations of the OFF and ON profiles are the same as in the previous
paragraph.

#### Changing voltage, pore length, and surface charge individually

Figs. 2, 5, and 6 showed cases for the same mDu to reveal how and why scaling
works. Next, we show profiles to explain how and why changing a single
parameter ($\sigma$, $H$, or $U$) results in a change in device behavior that
is similar in the three cases. For this, we change only one parameter while
all the others are kept fixed. This is shown in Fig. 7.

Fig. 7A shows the reciprocal profiles for the ON state (top panel), the OFF
state (middle panel), and their ratio (bottom panel). Blue and red colors
refer to cations and anions, respectively, while different shades of blue and
red refer to different surface charge densities. As $\sigma$ increases, the
peaks in the OFF-state $c_{i}^{-1}(z)$ profiles increase. Increasing
$c_{i}^{-1}(z)$, on the other hand means deepening depletion zones. The
$c_{i}^{\mathrm{ON}}(z)/c_{i}^{\mathrm{OFF}}(z)$ profile is very similar to
the OFF-state $c_{i}^{-1}(z)$ profile. This means that rectification is
chiefly determined by the behavior in the OFF state.

Fig. 7B is similar except that now the different shades of blue and red refer
to different pore lengths, $H$, and that we plot the profiles as functions of
$z/H$. Plotting this way reveals that increasing $H$ has little effect on the
ON profiles, while it produces larger peaks in the OFF-state $c_{i}^{-1}(z)$
profiles, e.g., deeper depletion zones.

The result that rectification is rather determined by the OFF state when we
change $H$ is also supported by Fig. 8 that shows the nanopores’ specific
conductances defined as

$\kappa=\frac{IH}{UR^{2}\pi}$ (24)

in the ON and OFF states. Because $\kappa$ is plotted on a logarithmic scale,
the distance between the two curves corresponds to rectification.
Rectification increases because the rate of the decrease of the OFF
conductance is larger than the rate of the increase of the ON conductance.

Figure 8: The specific conductance of the nanopore as a function of $H$ for
the parameters $R=1$ nm, $c=0.1$ M, $\sigma=1$ $e$/nm2, and $U=200$ mV. The
two curves represent the ON and OFF states, while their difference corresponds
to rectification due to the logarithmic scale of the ordinate.

As opposed to Figs. 7A and B that showed reciprocal concentration profiles,
Fig. 7C shows the concentration profiles for different voltages. This is
because changing voltage rather changes the ON-state concentration profiles.
Comparing the $c_{i}^{\mathrm{ON}}(z)/c_{i}^{\mathrm{OFF}}(z)$ profiles to
either the ON or the OFF profiles, we observe that the trend showed by them
(increasing peaks with increasing voltage) agrees with the trend shown by the
ON-state concentration profile.

Figure 9: The conductance of the nanopore as a function of $U$ for the
parameters $R=1$ nm, $c=0.1$ M, $\sigma=1$ $e$/nm2, and $U=200$ mV. The two
curves represent the ON and OFF states, while their ratio corresponds to
rectification.

The result that rectification is rather determined by the ON state when we
change $U$ is also supported by Fig. 9 that shows the nanopores’ total
conductances in the ON and OFF states computed as $g=g_{+}+g_{-}$ ($g_{i}$ was
defined in Eq. 21). It is apparent that rectification increases because the
rate of the increase of the ON-state conductance is larger than the rate of
the increase of the ON conductance.

##  IV Summary

We showed that by a little modification of the mDu parameter (Eq. 9)
introduced in our previous study sarkadi2021arxiv , a useful scaling parameter
can be obtained for bipolar nanopores by including the voltage. The case of
bipolar nanopores (bipolar surface charge pattern) is more complicated from
several points of view than the case of selective nanopores (uniform surface
charge density). The device function is rectification obtained from two
separate simulation for $U$ and $-U$. The nanotube limit,
$H\rightarrow\infty$, for which we obtained an exact scaling in Ref.
sarkadi2021arxiv , makes sense only as a long-pore limit here. The two
oppositely charged half pores “communicate” with each other; the pore charge
and the ionic distribution of one influences the ionic distribution of the
other.

Despite this complex picture and the fact that mDu contains all the relevant
parameters ($R$, $H$, $c$, $\sigma$, and $U$), we found a good, although not
perfectly accurate, scaling behavior. This behavior is useful not only from a
fundamental point of view, but also from a practical point of view. If we
know, from measurements, the properties of the nanopore for a given set of
parameters, we can predict it for another set of parameters as soon as scaling
is valid.

Scaling, however, is not necessarily a universal feature. Strong ionic
correlations present for multivalent ions may cause deviations from the smooth
monotonic scaling observed for 1:1 electrolytes. Therefore, the scaling
behavior as reported here and in our earlier publications madai_pccp_2018 ;
fertig_jpcc_2019 ; sarkadi2021arxiv  can serve as a gold standard the
deviations from which indicate the strength of ionic correlations.

##  Appendix A Local Equilibrium Monte Carlo

Let us divide the simulation cell into subvolumes $\mathcal{B}^{\alpha}$,
characterize this subvolume with the electrochemical potential,
$\mu_{i}^{\alpha}$, and assume that this value is constant in
$\mathcal{B}^{\alpha}$. We designate the $\mu_{i}^{\alpha}$ and
$c_{i}^{\alpha}$ value to the mass center of the $\mathcal{B}^{\alpha}$ volume
element. We assume that the subvolumes are open systems with a constant volume
($V^{\alpha}$), temperature ($T$), and electrochemical potential
($\mu_{i}^{\alpha}$) and that they are in local equilibrium. In the spirit of
GCMC simulations, we apply particle insertion/deletions with the acceptance
probability $\min(1;p^{\alpha}_{i,\chi}(\mathbf{r}))$, where

$p^{\alpha}_{i,\chi}(\mathbf{r})=\dfrac{N^{\alpha}_{i}!(V^{\alpha})^{\chi}}{(N^{\alpha}_{i}+\chi)!}\exp\left(-\dfrac{\Delta
U(\mathbf{r})-\chi\mu^{\alpha}_{i}}{kT}\right).$ (25)

Here, $N^{\alpha}_{i}$ is the number of ions of type $i$ in subvolume
$\mathcal{B}^{\alpha}$ before insertion/deletion, $\Delta U(\mathbf{r})$ is
the change of the system’s potential energy during particle insertion to
position $\mathbf{r}$ (or deletion from there), $\chi=1$ for insertion, and
$\chi=-1$ for deletion. The difference between this method and equilibrium
GCMC is that $\mu_{i}$ is space-dependent and the acceptance criterion is
referred to a given subvolume instead of the whole simulation cell. The effect
of the surrounding of subvolume $\mathcal{B}^{\alpha}$, however, is taken into
account in the simulation through the energy change that includes all the
interactions from other subvolumes, not only interactions between ions in
$\mathcal{B}^{\alpha}$.

Although it is tempting to consider the subvolume $\mathcal{B}^{\alpha}$ as a
distinct thermodynamic system with its own ensemble of states and to include
the effect of other subvolumes as an external constraint, this is not the
case. The ensemble of states belongs to the whole system because ion
configurations in subvolume $\mathcal{B}^{\alpha}$ should be collected for
every possible ion configurations of all the other subvolumes. Therefore, the
independent variables of this ensemble are $T$ and
$\\{V^{\alpha},\mu_{i}^{\alpha}\\}$, where $\alpha$ and $i$ run over the
volume elements and particle species, respectively. For comparison, the
variables in global equilibrium are $T$, $V$, and $\mu_{i}$, where $V$ is the
total volume and $\mu_{i}$ does not depend on space.

We solve the NP+LEMC system iteratively. The electrochemical potential is
adjusted until conservation of mass
($\nabla\cdot\mathbf{j}_{i}(\mathbf{r})=0$) is satisfied. The procedure can be
summarized as

$\mu^{\alpha}_{i}[n]\,\,\xrightarrow{\mathrm{LEMC}}\,\,c^{\alpha}_{i}[n]\,\,\xrightarrow{\mathrm{NP}}\,\,\mathbf{j}^{\alpha}_{i}[n]\,\,\xrightarrow{\nabla\cdot\mathbf{j}=0}\,\,\mu^{\alpha}_{i}[n+1].$
(26)

The electrochemical potentials for the next iteration,
$\mu_{i}^{\alpha}[n+1]$, are computed from the results of the previous
iteration, $c_{i}^{\alpha}[n]$, on the basis of the divergence-theorem (also
known as Gauss-Ostrogradsky’s theorem). The continuity equation is converted
to a surface integral:

$0=\int_{\mathcal{B}^{\alpha}}\nabla\cdot\mathbf{j}_{i}(\mathbf{r})\,dV=\oint_{\mathcal{S}^{\alpha}}\mathbf{j}_{i}(\mathbf{r})\cdot\mathbf{n}(\mathbf{r})\,da,$
(27)

where volume $\mathcal{B}^{\alpha}$ is bounded by surface
$\mathcal{S}^{\alpha}$ and $\mathbf{n}(\mathbf{r})$ denotes the normal vector
pointing outward at position $\mathbf{r}$ of the surface.

The system is solved by digitizing the problem. Every $\mathcal{S}^{\alpha}$
surface is divided into $\mathcal{S}^{\alpha\beta}$ elements. Along with these
elements, $\mathcal{B}^{\alpha}$ and $\mathcal{B}^{\beta}$ are adjacent cells.
It is assumed that the concentration, the gradient of the electrochemical
potential, the flux density and the diffusion coefficient are constant on a
surface $\mathcal{S}^{\alpha\beta}$. They are denoted by hat:
$\hat{c}^{\alpha\beta}_{i}$, $\nabla\hat{\mu}_{i}^{\alpha\beta}$,
$\hat{\mathbf{j}}_{i}^{\alpha\beta}$, and $\hat{D}_{i}^{\alpha\beta}$. The
$\hat{c}^{\alpha\beta}_{i}$ values are obtained from the values
$c_{i}^{\alpha}$ and $c_{i}^{\beta}$ via linear interpolation. The
$\nabla\hat{\mu}_{i}^{\alpha\beta}$ values are also obtained from
$\mu_{i}^{\alpha}$ and $\mu_{i}^{\beta}$ assuming linearity.

Thus the integral in Eq. 27 for a given surface $\mathcal{S}^{\alpha}$ is
replaced by a sum over the surface elements that constitute
$\mathcal{S}^{\alpha}$:

$0=\sum_{\beta,\mathcal{S}^{\alpha\beta}\in\mathcal{S}^{\alpha}}\hat{\mathbf{j}}_{i}^{\alpha\beta}\cdot\mathbf{n}^{\alpha\beta}a^{\alpha\beta},$
(28)

where $a^{\alpha\beta}$ is the area of surface element
$\mathcal{S}^{\alpha\beta}$ and $\mathbf{n}^{\alpha\beta}$ is the outward
normal vector in the center of $\mathcal{S}^{\alpha\beta}$. The iteration
procedure is described by the following steps:

  1. 1.

An appropriately chosen initial set of electrochemical potentials is chosen
($\mu_{i}^{\alpha}[1]$; in general: $\mu_{i}^{\alpha}[n]$, where $[n]$ denotes
the $n$th iteration).

  2. 2.

Using these $\mu_{i}^{\alpha}[n]$ parameters as inputs, LEMC simulations are
performed. The resulting concentrations are denoted by $c_{i}^{\alpha}[n]$.

  3. 3.

The flux computed from the $\\{\mu_{i}^{\alpha}[n],c_{i}^{\alpha}[n]\\}$ pair
usually does not satisfy Eq. 28. The next set of electrochemical potential is
calculated by assuming that the
$\\{\mu_{i}^{\alpha}[n+1],c_{i}^{\alpha}[n]\\}$ pair does satisfy Eq. 28. If
we write the value of $\hat{\mathbf{j}}_{i}^{\alpha\beta}$ as given by the NP
equation into Eq. 28, we obtain

$0=\sum_{\beta}\hat{D}_{i}^{\alpha\beta}\hat{c}_{i}^{\alpha\beta}[n]\nabla\hat{\mu}_{i}^{\alpha\beta}[n+1]\cdot\mathbf{n}^{\alpha\beta}\,a^{\alpha\beta},$
(29)

where $\beta$ runs over all the surface elements $\mathcal{S}^{\alpha\beta}$
that constitute $\mathcal{S}^{\alpha}$. Eq. 29 is a system of linear
equations; the unknown variables are denoted by
$\mu_{i}^{\alpha,\mathrm{CAL}}[n+1]$, where $\mathrm{CAL}$ refers to the fact
that these values come from calculations by solving 29.

  4. 4.

To achive faster and more robust convergence in the case of large driving
forces, the electrochemical potential used in the $(n+1)$th iteration is mixed
from the values calculated in the $(n+1)$th iteration from Eq. 29 and the
values mixed in the $n$th iteration:

$\displaystyle\mu_{i}^{\alpha,\mathrm{MIX}}[n+1]$ $\displaystyle=$
$\displaystyle b_{i}\mu_{i}^{\alpha,\mathrm{CAL}}[n+1]$ (30) $\displaystyle+$
$\displaystyle(1-b_{i})\mu_{i}^{\alpha,\mathrm{MIX}}[n],$

where $b_{i}$ is a mixing parameter that determines the ratio of mixing. If
the parameter is close to 1, faster iteration can be achieved, however, it may
result in the system fluctuating between local minima. Smaller $b_{i}$ values
prevent this fluctuation at the price of making the convergence slower. In
this work we used values $b=0.5-0.7$.

  5. 5.

The input of the $(n+1)$th LEMC simulation is
$\mu_{i}^{\alpha,\mathrm{MIX}}[n+1]$.

The energy change $\Delta U$ contains not only the interactions between
particles (and interactions of particles with the pore), but also the
interaction with an external electrical potential,
$\Phi^{\mathrm{APPL}}(\mathbf{r})$. This applied potential is calculated by
solving Laplace’s equation

$\nabla^{2}\Phi^{\mathrm{APPL}}(\mathbf{r})=0$ (31)

for the system (inside the blue line) with the prescribed Dirichlet boundary
condition on the system’s boundaries ($\Phi^{\mathrm{L}}$ and
$\Phi^{\mathrm{R}}$ on the left and right blue line, respectively).

In the rotationally symmetric present geometry studied here, the elementary
cells are $\Delta z\times\Delta r$ rectangles in the $(z,r)$ plane. In three-
dimensional space, these correspond to concentric rings with the $z$-axis in
their centers.

Because of statistical uncertainties of the LEMC simulations, the iteration
does not converge to an exact value of the ionic flux density, but it
fluctuates around a limiting value. The final solution is obtained from a
running average over iterations. Longer LEMC simulations and more iterations
result in a more reliable outcome. In this study, we used $40-100$ iterations
sampling $30$million configurations in each. This resulted in currents with an
error smaller than $1$ %.

## Acknowledgement(s)

We gratefully acknowledge the financial support of the National Research,
Development and Innovation Office – NKFIH K124353. We acknowledge KIFÜ for
awarding us access to resource based in Hungary at Szeged.

## References

  * (1) H. Daiguji, Y. Oka, and K. Shirono.  Nanofluidic diode and bipolar transistor.  Nano Lett., 5(11):2274–2280, 2005.
  * (2) D. Constantin and Z. S. Siwy.  Poisson-Nernst-Planck model of ion current rectification through a nanofluidic diode.  Phys. Rev. E, 76(4):041202, 2007.
  * (3) R. Karnik, C. Duan, K. Castelino, H. Daiguji, and A. Majumdar.  Rectification of ionic current in a nanofluidic diode.  Nano Lett., 7(3):547–551, 2007.
  * (4) I. Vlassiouk and Z. S. Siwy.  Nanofluidic diode.  Nano Lett., 7(3):552–556, 2007.
  * (5) E. B. Kalman, I. Vlassiouk, and Z. S. Siwy.  Nanofluidic bipolar transistors.  Adv. Mater., 20(2):293–297, 2008.
  * (6) I. Vlassiouk, S. Smirnov, and Z. Siwy.  Nanofluidic ionic diodes. comparison of analytical and numericalsolutions.  ACS Nano, 2(8):1589–1602, 2008.
  * (7) R. Yan, W. Liang, R. Fan, and P. Yang.  Nanofluidic diodes based on nanotube heterojunctions.  Nano Lett., 9(11):3820–3825, 2009.
  * (8) L.-J. Cheng and L. J. Guo.  Ionic current rectification, breakdown, and switching in heterogeneous oxide nanofluidic devices.  ACS Nano, 3(3):575–584, 2009.
  * (9) G. Nguyen, I. Vlassiouk, and Z. S Siwy.  Comparison of bipolar and unipolar ionic diodes.  Nanotech., 21(26):265301, 2010.
  * (10) W. Guo, Y. Tian, and L. Jiang.  Asymmetric ion transport through ion-channel-mimetic solid-state nanopores.  Acc. Chem. Res., 46(12):2834–2846, 2013.
  * (11) Z. Ható, M. Valiskó, T. Kristóf, D. Gillespie, and D. Boda.  Multiscale modeling of a rectifying bipolar nanopore: explicit-water versus implicit-water simulations.  Phys. Chem. Chem. Phys., 19(27):17816–17826, 2017.
  * (12) B. Matejczyk, M. Valiskó, M.-T. Wolfram, J.-F. Pietschmann, and D. Boda.  Multiscale modeling of a rectifying bipolar nanopore: Comparing Poisson-Nernst-Planck to Monte Carlo.  J. Chem. Phys., 146(12):124125, 2017.
  * (13) D. Fertig, B. Matejczyk, M. Valiskó, D. Gillespie, and D. Boda.  Scaling behavior of bipolar nanopore rectification with multivalent ions.  J. Phys. Chem. C, 123(47):28985–28996, 2019.
  * (14) L. Blum.  Mean spherical model for asymmetric electrolytes.  Mol. Phys., 30(5):1529–1535, 1975.
  * (15) L. Blum and J. S. Hoeye.  Mean spherical model for asymmetric electrolytes. 2. Thermodynamic properties and the pair correlation function.  J. Phys. Chem., 81(13):1311–1316, 1977.
  * (16) W. Nonner, L. Catacuzzeno, and B. Eisenberg.  Binding and Selectivity in L-Type Calcium Channels: A Mean Spherical Approximation.  Biophys. J., 79(4):1976–1992, 2000.
  * (17) Zs. Sarkadi, D. Fertig, Z. Ható, M. Valiskó, and D Boda.  From nanotubes to nanoholes: scaling of selectivity in uniformly charged nanopores through the Dukhin number for 1:1 electrolytes, 2021.
  * (18) T. Albrecht, T. Gibb, and P. Nuttall.  Ion transport in nanopores.  In Engineered Nanopores for Bioanalytical Applications, pages 1–30. Elsevier BV, 2013.
  * (19) P. Abgrall and N. T. Nguyen.  Nanofluidic devices and their applications.  Anal. Chem., 80(7):2326–2341, 2008.
  * (20) L. Bocquet and E. Charlaix.  Nanofluidics, from bulk to interfaces.  Chem. Soc. Rev., 39(3):1073–1095, 2010.
  * (21) H. Daiguji.  Ion transport in nanofluidic channels.  Chem. Soc. Rev., 39(3):901–911, 2010.
  * (22) J. C. T. Eijkel and A. van den Berg.  Nanofluidics and the chemical potential applied to solvent and solute transport.  Chem. Soc. Rev., 39(3):957, 2010.
  * (23) T. A. Zangle, A. Mani, and J. G. Santiago.  Theory and experiments of concentration polarization and ion focusing at microchannel and nanochannel interfaces.  Chem. Soc. Rev., 39(3):1014, 2010.
  * (24) E. Mádai, B. Matejczyk, A. Dallos, M. Valiskó, and D. Boda.  Controlling ion transport through nanopores: modeling transistor behavior.  Phys. Chem. Chem. Phys., 20(37):24156–24167, 2018.
  * (25) S. Dal Cengio and I. Pagonabarraga.  Confinement-controlled rectification in a geometric nanofluidic diode.  J. Chem. Phys., 151(4):044707, 2019.
  * (26) D. Fertig, M. Valiskó, and D. Boda.  Rectification of bipolar nanopores in multivalent electrolytes: effect of charge inversion and strong ionic correlations.  Phys. Chem. Chem. Phys., 22(34):19033–19045, 2020.
  * (27) M. Z. Bazant, K. Thornton, and A. Ajdari.  Diffuse-charge dynamics in electrochemical systems.  Phys. Rev. E, 70(2):021506, 2004.
  * (28) K. T. Chu and M. Z. Bazant.  Nonlinear electrochemical relaxation around conductors.  Phys. Rev. E, 74(1):011501, 2006.
  * (29) L. Bocquet and E. Charlaix.  Nanofluidics, from bulk to interfaces.  Chem. Soc. Rev., 39(3):1073–1095, 2010.
  * (30) J. J. Bikerman.  Electrokinetic equations and surface conductance. a survey of the diffuse double layer theory of colloidal solutions.  Trans. Farad. Soc., 35:154, 1940.
  * (31) S.S. Dukhin.  Non-equilibrium electric surface phenomena.  Adv. Coll. Interf. Sci., 44:1–134, 1993.
  * (32) J. J. Lyklema, A. de Keizer, B.H. Bijsterbosch, G.J. Fleer, and M.A. Cohen Stuart (Eds.).  Solid-Liquid Interfaces.  Fundamentals of Interface and Colloid Science 2. Elsevier, Academic Press, 1995.
  * (33) A. S. Khair and T. M. Squires.  Surprising consequences of ion conservation in electro-osmosis over a surface charge discontinuity.  J. Fluid Mech., 615:323–334, 2008.
  * (34) S. Das and S. Chakraborty.  Effect of conductivity variations within the electric double layer on the streaming potential estimation in narrow fluidic confinements.  Langmuir, 26(13):11589–11596, 2010.
  * (35) C. Lee, L. Joly, A. Siria, A.-L. Biance, R. Fulcrand, and L. Bocquet.  Large apparent electric size of solid-state nanopores due to spatially extended surface conduction.  Nano Lett., 12(8):4037–4044, 2012.
  * (36) H.-C. Yeh, M. Wang, C.-C. Chang, and R.-J. Yang.  Fundamentals and modeling of electrokinetic transport in nanochannels.  Israel J. Chem., 54(11-12):1533–1555, 2014.
  * (37) Y. Ma, J. Guo, L. Jia, and Y. Xie.  Entrance effects induced rectified ionic transport in a nanopore/channel.  ACS Sensors, 3(1):167–173, 2017.
  * (38) T. Xiong, K. Zhang, Y. Jiang, P. Yu, and L. Mao.  Ion current rectification: from nanoscale to microscale.  Sci. China Chem., 62(10):1346–1359, 2019.
  * (39) A. R. Poggioli, A. Siria, and L. Bocquet.  Beyond the tradeoff: Dynamic selectivity in ionic transport and current rectification.  J. Phys. Chem. B, 123(5):1171–1185, 2019.
  * (40) S. Dal Cengio and I. Pagonabarraga.  Confinement-controlled rectification in a geometric nanofluidic diode.  J. Chem. Phys., 151(4):044707, 2019.
  * (41) N. Kavokine, R. R. Netz, and L. Bocquet.  Fluids at the nanoscale: From continuum to subcontinuum transport.  Annu. Rev. Fluid Mech., 53(1), 2020.
  * (42) Y. Noh and N. R. Aluru.  Ion transport in electrically imperfect nanopores.  ACS Nano, 14(8):10518–10526, 2020.
  * (43) A. Levy, J. P. de Souza, and M. Z. Bazant.  Breakdown of electroneutrality in nanopores.  J. Coll. Inter. Sci., 579:162–176, 2020.
  * (44) Z. Siwy, E. Heins, C. C. Harrell, P. Kohli, and C. R. Martin.  Conical-nanotube ion-current rectifiers: the role of surface charge.  J. Am. Chem. Soc., 126(35):10850–10851, 2004.
  * (45) S. Garaj, W. Hubbard, A. Reina, J. Kong, D. Branton, and J. A. Golovchenko.  Graphene as a subnanometre trans-electrode membrane.  Nature, 467(7312):190–193, 2010.
  * (46) S. Garaj, S. Liu, J. A. Golovchenko, and D. Branton.  Molecule-hugging graphene nanopores.  Proc. Nat. Acad. Sci., 110(30):12192–12196, 2013.
  * (47) S. C. O’Hern, M. S. H. Boutilier, J.-C. Idrobo, Y. Song, J. Kong, T. Laoui, M. Atieh, and R. Karnik.  Selective ionic transport through tunable subnanometer pores in single-layer graphene membranes.  Nano Lett., 14(3):1234–1241, 2014.
  * (48) R. C. Rollings, A. T. Kuan, and J. A. Golovchenko.  Ion selectivity of graphene nanopores.  Nat. Comm., 7(1), 2016.
  * (49) J. P. Thiruraman, K. Fujisawa, G. Danda, P. M. Das, T. Zhang, A. Bolotsky, N. Perea-López, A. Nicolaï, P. Senet, M. Terrones, and M. Drndić.  Angstrom-size defect creation and ionic transport through pores in single-layer MoS2.  Nano Lett., 18(3):1651–1659, 2018.
  * (50) J. P. Thiruraman, P. M. Das, and M. Drndić.  Stochastic ionic transport in single atomic zero-dimensional pores.  ACS Nano, 14(9):11831–11845, 2020.
  * (51) E. Mádai, M. Valiskó, A. Dallos, and D. Boda.  Simulation of a model nanopore sensor: Ion competition underlines device behavior.  J. Chem. Phys., 147(24):244702, 2017.
  * (52) A. Malasics and D. Boda.  An efficient iterative grand canonical monte carlo algorithm to determine individual ionic chemical potentials in electrolytes.  J. Chem. Phys., 132(24):244103, 2010.
  * (53) D. Boda and D. Gillespie.  Steady state electrodiffusion from the Nernst-Planck equation coupled to Local Equilibrium Monte Carlo simulations.  J. Chem. Theor. Comput., 8(3):824–829, 2012.
  * (54) A. Szymczyk, H. Zhu, and B. Balannec.  Ion rejection properties of nanopores with bipolar fixed charge distributions.  J. Phys. Chem. B, 114(31):10143–10150, aug 2010.
  * (55) K. P. Singh and M. Kumar.  Effect of surface charge density and electro-osmotic flow on ionic current in a bipolar nanopore fluidic diode.  J. Appl. Phys., 110(8):084322, 2011.
  * (56) K. P. Singh and M. Kumar.  Effect of nanochannel diameter and debye length on ion current rectification in a fluidic bipolar diode.  J. Phys. Chem. C, 115(46):22917–22924, 2011.
  * (57) K. P. Singh, K. Kumari, and M. Kumar.  Ion current rectification in a fluidic bipolar nanochannel with smooth junction.  Appl. Phys. Lett., 99(11), 2011.
  * (58) L. van Oeffelen, W. Van Roy, H. Idrissi, D. Charlier, L. Lagae, and G. Borghs.  Ion current rectification, limiting and overlimiting conductances in nanopores.  PLoS ONE, 10(5):e0124171, 2015.
  * (59) M. Tajparast, G. Virdi, and M. I. Glavinović.  Spatial profiles of potential, ion concentration and flux in short unipolar and bipolar nanopores.  Biochim. Biophys. Acta (BBA) - Biomem., 1848(10, Part A):2138–2153, 2015.
  * (60) D. Gillespie, W. Nonner, and R. S. Eisenberg.  Coupling Poisson-Nernst-Planck and density functional theory to calculate ion flux.  J. Phys. Condens. Matter, 14(46):12129–12145, 2002.
  * (61) D. Gillespie, W. Nonner, and R. S. Eisenberg.  Density functional theory of charged, hard-sphere fluids.  Phys. Rev. E, 68(3):031503, 2003.
  * (62) D. Gillespie, L. Xu, Y. Wang, and G. Meissner.  (De)constructing the Ryanodine Receptor: Modeling Ion Permeation and Selectivity of the Calcium Release Channel.  J. Phys. Chem. B, 109(32):15598–15610, 2005.
  * (63) D. Gillespie.  Energetics of divalent selectivity in a calcium channel: The Ryanodine Receptor case study.  Biophys. J., 94(4):1169–1184, 2008.
  * (64) M. Burger, B. Schlake, and M.-T. Wolfram.  Nonlinear Poisson–Nernst–Planck equations for ion flux through confined geometries.  Nonlinearity, 25(4):961–990, 2012.
  * (65) R. Eisenberg.  PNP What is in a name?  aug 2019.
  * (66) D. Boda.  In R. A. Wheeler, editor, Ann. Rep. Comp. Chem., volume 10, chapter 5 Monte Carlo Simulation of Electrolyte Solutions in Biology: In and Out of Equilibrium, pages 127–163. Elsevier, 2014.
  * (67) D. Boda, R. Kovács, D. Gillespie, and T. Kristóf.  Selective transport through a model calcium channel studied by Local Equilibrium Monte Carlo simulations coupled to the Nernst-Planck equation.  J. Mol. Liq., 189:100–112, 2014.
  * (68) D. Fertig, E. Mádai, M. Valiskó, and D. Boda.  Simulating ion transport with the NP+LEMC method. Applications to ion channels and nanopores.  Hung. J. Ind. Chem., 45(1):73–84, 2017.
  * (69) D. Gillespie and D. Boda.  The anomalous mole fraction effect in calcium channels: A measure of preferential selectivity.  Biophys. J., 95(6):2658–2672, 2008.
  * (70) D. Gillespie, D. Boda, Y. He, P. Apel, and Z.S. Siwy.  Synthetic nanopores as a test case for ion channel theories: The anomalous mole fraction effect without single filing.  Biophys. J., 95(2):609–619, 2008.
  * (71) Y. He, D. Gillespie, D. Boda, I. Vlassiouk, R. S. Eisenberg, and Z. S. Siwy.  Tuning transport properties of nanofluidic devices with local charge inversion.  JACS, 131(14):5194–5202, 2009.
  * (72) D. Boda, M. Valiskó, D. Henderson, B. Eisenberg, D. Gillespie, and W. Nonner.  Ion selectivity in L-type calcium channels by electrostatics and hard-core repulsion.  J. Gen. Physiol., 133(5):497–509, 2009.
  * (73) M. Malasics, D. Boda, M. Valiskó, D. Henderson, and D. Gillespie.  Simulations of calcium channel block by trivalent ions: Gd3+ competes with permeant ions for the selectivity filter.  Biochim. et Biophys. Acta - Biomembranes, 1798(11):2013–2021, 2010\.
  * (74) M. Valiskó, B. Matejczyk, Z. Ható, T. Kristóf, E. Mádai, D. Fertig, D. Gillespie, and D. Boda.  Multiscale analysis of the effect of surface charge pattern on a nanopore’s rectification and selectivity properties: from all-atom model to poisson-nernst-planck.  J. Chem. Phys., 150(14):144703, 2019.
  * (75) D. Boda, M. Valiskó, and D. Gillespie.  Modeling the device behavior of biological and synthetic nanopores with reduced models.  Entropy, 22(11):1259, 2020.
  * (76) S. Levine, J. R. Marriott, G. Neale, and N. Epstein.  Theory of electrokinetic flow in fine cylindrical capillaries at high zeta-potentials.  J. Coll. Interf. Sci., 52(1):136–149, 1975.
  * (77) S. Balme, F. Picaud, M. Manghi, J. Palmeri, M. Bechelany, S. Cabello-Aguilar, A. Abou-Chaaya, P. Miele, E. Balanzat, and J. M. Janot.  Ionic transport through sub-10 nm diameter hydrophobic high-aspect ratio nanopores: experiment, theory and simulation.  Sci. Rep., 5(1), 2015.
  * (78) Y. Uematsu, R. R. Netz, L. Bocquet, and D. J. Bonthuis.  Crossover of the power-law exponent for carbon nanotube conductivity as a function of salinity.  J. Phys. Chem. B, 122(11):2992–2997, 2018.
  * (79) Y. Green.  Ion transport in nanopores with highly overlapping electric double layers.  J. Chem. Phys., 154(8):084705, 2021.