Current rectification by nanoparticles in bipolar nanopores

Two bipolar nanopore geometries are examined in what follows, namely a pore in which the two oppositely charged sections have the same length and wall charge density, and another in which the section with positive wall charges is four times shorter and has a wall charge density three times larger than the section with negative wall charges. In the first nanopore, referred to as C1, the section with positively charged walls has a charge density of $\sigma_{1}=+0.1\penalty\ C/{\rm m}^{2}$. The section with negatively charged walls has a wall charge density of $\sigma_{2}=-0.1\penalty\ C/{\rm m}^{2}$. The wall charge density is $\sigma_{\neq 1,2}=0\penalty\ C/{\rm m}^{2}$ for the rest of the walls. In the second channel, which is referred to as C2, the region with positively charged walls has a wall charge density of $\sigma_{1}=+0.24\penalty\ C/{\rm m}^{2}$, and the section with negatively charged walls has a wall charge density of $\sigma_{2}=-0.08\penalty\ C/{\rm m}^{2}$. In C2 the wall charge density is also $\sigma_{\neq 1,2}=0\penalty\ C/{\rm m}^{2}$ for the rest of the walls. The magnitudes of the surfaces charges adopted here are comparable to those found in experimental studies with Janus membranes Gao et al. (2014); Zhang et al. (2015); Zhang et al. (2017); Zhu et al. (2018). Both geometries feature changes in the height as a function of channel length. In both channels, the section with positively charged walls has a height of $17$ nm and the section with negatively charged walls has a height of $10$ nm. When nanoparticles are added to the nanopores, the ionic strength of the solution ($\approx 0.1$ M) is kept constant. The Debye length of the electrolytic solution is $\lambda_{\rm D}\approx 0.95\penalty\ \text{nm}$. Which is an order of magnitude smaller than the height in the narrowest section of the nanochannels considered here. Plots of the ionic concentrations as a function of the channel’s height in the wide and narrow sections of the nanopore are given in Figure S1 in the Supporting Information (SI).

We consider the addition of 10 nanoparticles in both channels. Moreover the effect of changing the surface charge of the nanoparticles from positive to negative is also studied here. Each nanoparticle has a hydrodynamic radius of about 1.14 nm and a surface charge density of about $\pm 0.33\penalty\ \text{C}/\text{m}^{2}$. See the Methods section below for more details on the nanoparticle model used here. The contributions to the ionic strength from the cations, anions and nanoparticles for the different cases considered here are summarized in Table 1.

Figure 1 shows steady state orthographic projections of the two bipolar nanopores as predicted by the coarse-grained molecular dynamics (MD) simulations. Sodium ions are represented as red spheres and Chlorine ions as blue spheres and the black spheres represent the interaction sites that make up the channel walls and the nanoparticles. The cases shown in Figure 1 correspond to nanopores where positively charged nanoparticles have been added. The same channels but with negatively charged nanoparticles can be found in Figure S2 of the SI. Figures 1 A and D correspond to C1 and C2, respectively, when no voltage bias is applied. It can be observed that there is an obvious excess of anions in the wider sections of the nanopores, where the wall has a positive charge. In the narrow section of the nanopores, where the wall charge is negative, a excess of cations is evident. In C1 the nanoparticles are distributed more or less evenly between the wide and narrow sections of the pore. In C2 the nanoparticles mostly remain in the long narrow section. With respect to the height of the channels, the location of the nanoparticles does not exhibit major differences between C1 and C2. The positively charged nanoparticles reside near the wall in the section of the pores with negative wall charge. In the section of the pores with positive wall charge the positively charged nanoparticles tend to stay near the center of the channels. As expected, the negatively charged nanoparticles do the opposite.

Figures 1 B and E show C1 and C2, respectively, when a negative bias, $-2\penalty\ \text{V}$, is applied. The details on how the voltage bias between the two ends of the nanopores is applied in the simulations are given in the Methods section below. Also note that the magnitude of the voltage bias used here is what has been commonly used in experiments with Janus membranes Gao et al. (2014); Zhang et al. (2015); Zhang et al. (2017); Zhu et al. (2018). In this case, a pronounced depletion zone, where the amount of both anions and cations drops, is observed around the transition from the narrow to the wide section of the nanopores. For this negative bias the general location of the nanoparticles inside the channels is not significantly changed compared to the equilibrium cases. However, in C1 it is more apparent that, similarly to what the ions do, the nanoparticles also stay away from the aforementioned depletion zone. Figures 1 C and F show C1 and C2, respectively, when a positive bias $+2\penalty\ \text{V}$, is applied. Here, contrary to what happens in the negative bias case, anions and cations appear to concentrate around the transition zone between the narrow and wide sections of the channels. For the positive bias, the positively charged nanoparticles move towards the wide section of the nanopores. They tend to concentrate around the transition zone between the narrow and wide sections of the cahnnels. As can be seen in Figure S3 in the SI, the negatively charged nanoparticles move towards the narrow section when the positive bias is applied.

The top row of Figure 2 shows concentration profiles as a function of the channel’s length (i.e. averaged over the channel’s width and depth). for (A) cations and (B) anions in C1 without nanoparticles (lines) and in C1 with positively charged nanoparticles (symbols). The bottom row of Figure 2 shows concentration profiles as a function of the channel’s length for (A) cations and (B) anions in C2 without nanoparticles and in C2 with positively charged nanoparticles (symbols). The same plots but for the cases where the nanoparticles have a negative surface charge are provided in the SI. As expected, the concentration of cations is larger in the narrow section of the slits, where the walls are negatively charged. In the wide section of the slits, where the wall is positively charged, the concentration of cations is smaller. A depletion zone is observed where the concentration of cations drops nearly to zero in the transition region between the narrow and wide sections of the nanopores. This depletion zone is larger for the $-2\penalty\ \text{V}$ bias case, and is also present in the equilibrium case. However, it is almost non-existent in the $+2\penalty\ \text{V}$ bias case. It has been previously reported that both cations and anions exhibit depletion zones in bipolar nanochannels, and that the depth of these depletion zones is sensitive to sign of the external voltage Fertig et al. (2020). The presence of either the negatively or the positively charged nanoparticles increases the length and depth of this depletion zone for the cations in both channels. It is also evident from Figure 2 that the cation depletion zone is longer and deeper in C1. Furthermore, it is the longest and deepest for C1 with a $-2\penalty\ \text{V}$ applied bias. For the $+2\penalty\ \text{V}$ applied bias no depletion zone is observed for the cations, and the maximum concentration occurs in the transition region between the narrow section and the wide section of the nanopores.

As can be seen in Figures 2 B and D, the concentration of anions is larger in the wide section of the slits, where the walls are positively charged. In the narrow section of the slits, where the wall is negatively charged, the concentration of anions is smaller. A depletion zone is also observed where the concentration of anions drops nearly to zero in the transition region between the narrow section and the wide section of the pore. This depletion zone again is observed at equilibrium, but is most pronounced in the $-2\penalty\ \text{V}$ case. As before, nanoparticles make the anion depletion zone significantly longer and deeper, specially in the C1 nanopore. The highest concentrations of anions occur for a $+2\penalty\ \text{V}$ bias applied, and they concentrate near the transition area between the narrow and wide sections of the slit. Note that they stay primarily on the side where the wall is positively charged.

From a practical perspective, the purpose of computer simulations is to design bipolar nanopores that can produce a usable current for a given applied voltage bias. Several theoretical and computational works have proposed current-voltage relationships for bipolar nanopores Mafé et al. (1990); Green et al. (2015); Moy et al. (2000); Corry et al. (2000); Long et al. (2018); Valiskó et al. (2019); Fertig et al. (2019); Fertig et al. (2020). Figure 3 shows the current-voltage relationship for the two slit pores considered here, as predicted by our coarse-grained MD calculations. Details on how the electric current is calculated are given in the Methods section below. In particular, Figure 3A shows the effect of adding nanoparticles (positive or negative) on the current-voltage relationship for the C1 nanopore. It can be observed that, for the $-2\penalty\ \text{V}$ applied bias, the electric current is suppressed by the presence of positive and negative charged nanoparticles. For the $+2\penalty\ \text{V}$ bias case, the effect of nanoparticles on the current is minor, and falls within the uncertainty of the calculations. The electric current rectification for C1 with positively charged nanoparticles is 22.3, which is about 13 times larger than without nanoparticles. The electric current rectification for C1 with negatively charged nanoparticles is 12.3, which is about 7 times larger than without nanoparticles.

Figure 3B shows the effect of adding nanoparticles, positively or negatively charged on the current-voltage relationship for the C2 nanopore. First, note that without nanoparticles the rectification ratios in C2, $|I_{\text{+2V}}^{\text{C2}}|/|I_{\text{-2V}}^{\text{C2}}|=1.82$, and C1 $|I_{\text{+2V}}^{\text{C1}}|/|I_{\text{-2V}}^{\text{C1}}|=1.78$ are very similar. However, when nanoparticles (either positively or negatively charged) are added to C2, the rectification ratio increases only by a factor of about 1.2. This represents a much smaller increase than the factor of 13 that is observed in C1 when positively charged nanoparticles are added. A summary of the rectification factors is given in Table 2.

Figure 3 also shows comparisons of the current-voltage relation from the MD simulations used here with predictions obtained using the analytic model for bipolar nanopores proposed by Green et al. (2015). The analytic model is based on the Poisson-Nernst-Planck (PNP) equations. It’s main assumptions are local electroneutrality in the reservoirs, cross-sectional local electroneutrality in the permselective regions and equal concentration of the electrolyte in the two reservoirs. The analytic predictions of that model have been verified by solving the fully coupled PNP equations using the finite elements numerical method. The analytic expression for the electric current has three distinct terms, one term accounts for the reservoirs’ Ohmic resistors, which also account for field focusing into the permselective interfaces. Another term describes the electric potential drop over the permselective regions. A third term describes the Donnan potential jumps at all three interfaces. The only adjustable parameter of the model is a diffusion coefficient for the ions.

For the predictions shown in Figure 3 the value obtained from fitting the analytic model to the electric currents obtained in the simulations is $\mathcal{D}=1.0\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$. This value can be contrasted with the diffusion coefficients of the ions inside the nanopores measured in the simulations. These are $\mathcal{D}_{\rm Na^{+}}\approx 9.05\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$ and $\mathcal{D}_{\rm Cl^{-}}\approx 5.85\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$. See the Methods section below for details on how these diffusion coefficients are calculated. This means that to predict similar currents the diffusion coefficient in the analytic model has to be about 5 to 9 times smaller than the diffusion coefficients in the MD simulations. The difference is reasonable given the different levels of description of the two models; the less detailed level of the continuum analytical model and the more detailed level of the MD simulations. Note That the analytical model is only applicable in the cases where there are no nanoparticles added. But in those cases the electric current shows agreement with the results obtained with the MD simulations employed here.

We have also examined the effect of the nanoparticle’s friction coefficient on the electric current. This was done by performing additional replicas of the simulations but using a friction coefficient for the nanoparticles about 8 times smaller than the one used for the results shown in Figure 3. See the Methods section below for details. Figure S10 and Table S3 in the SI show comparisons of the electric currents and the rectification ratios obtained in the different replicas performed for nanopore C1 with positively charged nanoparticles. To within the statistical uncertainty of the simulations all the replicas give the same electric currents and rectification ratios.

In the simulations presented here, the ionic strength of the solutions has been kept constant when adding the nanoparticles. Therefore one can immediately rule out that the observed increases in current rectification are due to an increase or decrease in the concentration of the total charges in the systems. Also, the concentration of particles has been kept very low, at about $0.7$ mM, and their mobility is about ten times smaller than the mobility of the ions (see the Methods section below). Therefore the contribution of the nanoparticles to the total electric current flowing through the nanopores is always less than $1\%$ (see Tables S4–S6 in the SI for the individual contributions of each charged species to the total electric currents). Moreover, changing the friction coefficient of the nanoparticles by a factor of 8 does not have a significant effect on the observed rectification ratios. Nevertheless, the nanoparticles can have a significant effect on the charge distribution inside the nanopores and therefore can produce significant changes in the nanopore’s electric potential.

Several theoretical studies based on the PNP equations Mauro (1962); Coster (1965); Sonin and Grossman (1972); Green et al. (2014); Chinaryan et al. (2014); Green et al. (2015); Green and Yossifon (2015), have shown that in a bipolar nanopore, the asymmetric and diode-like behavior of the current is due to the system layout. When the voltage is applied in a direction that coincides with the inherent electric field of the nanopore, then the current is enhanced. When the voltage is applied in the opposite direction, a so-called reverse bias, then the applied electric field operates against the inherent electric field, leading to a reduction in the resultant current Green et al. (2015). The effects of adding positively charged nanoparticles or negatively charged nanoparticles on the electric field are shown in Figures 4 A and B, respectively. Results are shown for the electric potential in C1 as a function of channel length, $\psi(x)$. The details on how the electric field of the nanopores as a function of length, $\psi(x)$, is calculated are given in the Methods section below. Results for C2 are given in Figure S6 of the SI.

One can clearly observe that the nanoparticles have a large effect on the electric potential. More specifically, in C1 negatively charged nanoparticles increase the depth of the minimum from $-1.6\penalty\ \text{V}$ to $-1.7\penalty\ \text{V}$, and positively charged nanoparticles increase the depth of the minimum from $-1.6\penalty\ \text{V}$ to $-3.3\penalty\ \text{V}$. For the $-2\penalty\ \text{V}$ applied bias case, negatively charged nanoparticles increase the depth of the minimum from $-4.2\penalty\ \text{V}$ to $-4.8\penalty\ \text{V}$. For this negative bias case, positively charged nanoparticles increase the depth of the minimum from $-4.2\penalty\ \text{V}$ to $-6.5\penalty\ \text{V}$. For the $+2\penalty\ \text{V}$ applied bias case, negatively charged nanoparticles decrease the depth of minimum in the electric potential from $-1.5\penalty\ \text{V}$ to $-1.3\penalty\ \text{V}$. For this positive bias case, positively charged nanoparticles increase the depth of the minimum from $-1.5\penalty\ \text{V}$ to $-2.6\penalty\ \text{V}$. In tables S4–S6 of the SI it can be observed that the deeper minimum in the electric potential observed for $-2\penalty\ \text{V}$ bias in the presence of nanoparticles significantly reduces the fluxes of anions and cations in that case. The rectification ratio is increased more by the positively charged nanoparticles because these produce a much deeper minimum in the electric potential for the $-2\penalty\ \text{V}$ bias, thereby causing the current to be diminished for that bias while maintaining the current enhancement for the $+2\penalty\ \text{V}$ bias.

In C2, nanoparticles do not have a significant effect on the inherent electric field of the nanopore. Overall, they tend to slightly increase the depth of the minimum in the electric potential. This leads to a slight overall reduction in the electric current that is observed in that pore upon addition of nanoparticles.

The coarse-grained, CG, parameters of the ion-ion effective potentials were obtained using relative entropy coarse-graining, RE-CG, Shell (2008), which provides a systematic approach to obtaining effective potentials for use in CG simulations, given a mapping function from all-atom, AA, to CG coordinates Hinckley and de Pablo (2015). The CG non-bonded effective potentials between ions are represented as the sum of a Coulombic interaction and a correction term, $U_{\rm corr}$, as follows,

Here $q_{i}$ and $q_{j}$ are the charges of the $i^{\rm th}$ and $j^{\rm th}$ ion. The dielectric permittivity of vacuum is given by $\epsilon_{0}$. Here $\epsilon(T)=249.4-0.788T/{\rm K}+7.20\times 10^{-4}(T/{\rm K})^{2}$ is the solution dielectric, $T$ is the absolute temperature and $r_{ij}$ is the inter-ion separation. $U_{\rm corr}$ is represented using cubic splines, and corrects the Coulombic potential to account for the effects of hydration and stearic overlap. The use of cubic splines permits the reproduction of subtle effects that could not be easily resolved with analytical expressions. The effective potentials recreate essentially all features of the radial distribution functions from AA simulations Hinckley and de Pablo (2015). The effective soft-core repulsion (i.e. radius) that arises after systematic coarse-graining is $0.2494\penalty\ \text{nm}$ for Sodium (Na⁺) ions and $0.4478\penalty\ \text{nm}$ for Chlorine (Cl^-) ions.

The nanopores considered here are channels with a slit geometry that is periodic in the $x$ and $y$ coordinates. The total size of the simulation box is set to $70\penalty\ \text{nm}\times\penalty\ 20\penalty\ \text{nm}\times 30.1\penalty\ \text{nm}$. The height of the channels and the wall charge density as a function of length, $\left[h(x),\sigma(x)\right]$, is given by,

Two types of bipolar nanochannels were simulated. One in which the two oppositely charged sections have the same length and another one in which the positively charged section is about four times shorter than the negatively charged section. The first channel is referred to as C1. In it $\ell=5\penalty\ \text{nm}$ in eq. (5.2). In this channel for $-30\penalty\ {\rm nm}\leq x\leq-5\penalty\ {\rm nm}$ the wall charge density is $\sigma_{1}=+0.1\penalty\ C/{\rm m}^{2}$ and the surface area of the wall in this section is $S_{1}=1000\penalty\ \text{nm}^{2}$. For $-5\penalty\ {\rm nm}\leq x\leq 20\penalty\ {\rm nm}$ the wall charge density is $\sigma_{2}=-0.1\penalty\ C/{\rm m}^{2}$, and the surface area of the wall in this section is also $S_{2}=1000\penalty\ \text{nm}^{2}$. The wall charge density is $\sigma_{\neq 1,2}=0\penalty\ C/{\rm m}^{2}$ for the rest of the walls, including the walls between sections of different height and walls that are perpendicular to the $x$ axis. The sections at the end of the channel, where the wall charge is zero, are intended to act as reservoirs. The total volume of C1 is $V_{C1}=25500\penalty\ \text{nm}^{3}$.

The second channel is denoted C2 and in it the negatively charged section is four times longer than the positively charged section. Therefore $\ell=20\penalty\ \text{nm}$ in eq. (5.2). And for $-30\penalty\ {\rm nm}\leq x\leq-20\penalty\ {\rm nm}$ the wall charge density is $\sigma_{1}=+0.24\penalty\ C/{\rm m}^{2}$ and the surface area of the wall in this section is $S_{1}=400\penalty\ \text{nm}^{2}$. For $-20\penalty\ {\rm nm}\leq x\leq 20\penalty\ {\rm nm}$ the wall charge density is $\sigma_{2}=-0.08\penalty\ C/{\rm m}^{2}$, and the surface area of the wall in this section is $S_{2}=1600\penalty\ \text{nm}^{2}$. The wall charge density is $\sigma_{\neq 1,2}=0\penalty\ C/{\rm m}^{2}$ for the rest of the walls, including the walls between sections of different height and walls that are perpendicular to the $x$ axis. The sections at the end of the channel, where the wall charge is zero, are intended to act as reservoirs. The total volume of C2 is $V_{C2}=23400\penalty\ \text{nm}^{3}$.

The systems are always globally neutral, i.e., $S_{1}\sigma_{1}+S_{2}\sigma_{2}=z_{c}en_{c}+z_{a}en_{a}+z_{p}en_{p}$. Where $n_{a}$ is the total number of anions (chloride) in the simulation box and $ez_{a}=-e$ is the charge of each anion. $e\approx 1.6\times 10^{-19}\penalty\ \text{C}$ is the elementary charge. $n_{c}$ is the total number of cations (Sodium) and $ez_{c}=+e$ is the charge of each cation. $n_{p}$ is the total number of nanoparticles in the simulation box and $ez_{p}$ is the charge of each nanoparticle. The Debye length which is used to estimate the size of the electrical double layer is $\lambda_{\rm D}=\sqrt{\dfrac{\epsilon_{0}\epsilon Vk_{\rm B}T}{e^{2}\sum_{i=1}^{N}z_{i}^{2}}}$. Where $V$ is the volume of the pore and $N$ is the total number of charged particles in the solution that flows through the pore. In all the simulations presented in this work $\lambda_{\rm D}\approx 0.95\penalty\ \text{nm}$.

In our coarse grained MD simulations the channel walls are modeled with a combination of neutral and charged sites that are kept fixed in space at their initial position. The charged wall sites are placed at random positions along the walls. The proportion of charged to neutral sites in a (charged) wall is determined by the wall charge density, the surface area of the wall, and by the charge of the wall charged sites. Each charged wall site is assigned a charge of $\pm 0.25e$ and interacts with the ions through Coulomb interactions.

Additionally, all wall sites, charged and uncharged, interact with the ions in solution through a repulsive Lenard-Jones potential,

where $\epsilon=1.68\penalty\ k_{\rm B}T$, $\sigma=0.4$ nm. A cutoff of $r_{\rm cut}=2^{1/6}\sigma$ is applied, and $r_{ij}$ is the distance between the ion and the wall site.

To apply the electric potential bias between the two ends of the bipolar nanopores we use a method commonly used in MD simulations for this purpose Gumbart et al. (2012). The method consists of applying a uniform electric field throughout the entire simulated periodic cell containing the nanopore. The external electric field $\bm{E}(\bm{r})$ is related to the external electric potential, $\psi_{E}(\bm{r})$, by $\bm{E}(\bm{r})=-\nabla\psi_{E}(\bm{r})$. In our simulations the external electric field is applied in the $x$ direction, $\bm{E}(\bm{r})=E_{x}\bm{\delta}_{x}$, where $\bm{\delta}_{x}$ is the unit vector in the $x$ direction and $E_{x}=\pm 0.0286\penalty\ \text{V}/\text{nm}$. The total length of the pores in the $x$ direction is $L=70\penalty\ \text{nm}$. The voltage bias is defined as $\Delta V=\psi_{E}(x=30\penalty\ \text{nm})-\psi_{E}(x=-40\penalty\ \text{nm})$. Therefore when $E_{x}=-0.0286\penalty\ \text{V}/\text{nm}$, $\Delta V=-E_{x}/L=+2\penalty\ \text{V}$. And when $E_{x}=+0.0286\penalty\ \text{V}/\text{nm}$, $\Delta V=-E_{x}/L=-2\penalty\ \text{V}$.

The coarse grained MD simulations were performed using LAMMPS Plimpton (1995) in the NVT ensemble using the Gronbench-Jensen Farago (G-JF) Langevin thermostat Grønbech-Jensen and Farago (2013); Grønbech-Jensen et al. (2014), with a damping factor of $\tau_{\rm m}=10\times 10^{-15}\penalty\ \text{s}$. The friction coefficient of ion $i$, $\zeta_{i}$, is related to the damping factor $\tau_{\rm m}$ by $\zeta_{i}=\dfrac{m_{i}}{\tau_{\rm m}}$, where $m_{i}$ is the mass of each ion. For Sodium $m_{{\rm Na}^{+}}=3.82\times 10^{-23}\penalty\ \text{g}$ (i.e. $22.9898\penalty\ \text{g/mol}$) and for Chlorine $m_{{\rm Cl}^{-}}=5.889\times 10^{-23}\penalty\ \text{g}$ (i.e. $35.453\penalty\ \text{g/mol}$). Therefore $\zeta_{{\rm Na}^{+}}=3.82\times 10^{-9}\penalty\ \text{N s}/\text{m}$ and $\zeta_{{\rm Cl}^{-}}=5.89\times 10^{-9}\penalty\ \text{N s}/\text{m}$. A time step $\Delta t=2\times 10^{-15}\penalty\ \text{s}$ and a temperature of $297\penalty\ \text{K}$ were used in all simulations. We relied on the particle mesh Ewald solver to calculate electrostatic interactions, with a real space cutoff of 1.2 nm. We used the modification of the three-dimensional Ewald summation technique for calculation of long-range Coulombic forces in systems with a slab geometry that are periodic in two dimensions and have a finite length in the third dimensionYeh and Berkowitz (1999).

Each nanoparticle is modeled by placing eight charged sites in the corners of a cube with sides of length $0.8\penalty\ \text{nm}$. Each site in the nanoparticle is assigned a charge of $\pm e$, therefore the total charge of each nanoparticle is $ez_{p}=-8e$. In our simulations, the relative distance between these sites is held approximately constant by using stiff harmonic springs. More specifically, adjacent sites are bonded by a spring with an equilibrium length of $0.8\penalty\ \text{nm}$ and a spring constant of $1.2\times 10^{5}\penalty\ k_{\rm B}T/{\rm nm}^{2}$. Sites in opposite corners of the cube are bonded by a harmonic potential with equilibrium length of $1.39\penalty\ {\rm nm}$ and an energy barrier of $1.2\times 10^{5}\penalty\ k_{\rm B}T/{\rm nm}^{2}$. Triplets of adjacent sites in the nanoparticles are also constrained with an angle potential with equilibrium value of $90^{\circ}$ and an energy barrier of $5.1\times 10^{2}\penalty\ k_{\rm B}T$. To estimate the charge density of the nanoparticles we consider a cube of side length $0.8$ nm. Therefore each nanoparticle has a surface area of $3.84\penalty\ \text{nm}^{2}$. Since each nanoparticle has eight charged sites with a charge of $\pm e$ the surface charge density is $\pm 0.33\penalty\ \text{C}/\text{m}^{2}$. The total mass of each nanoparticle is $m_{p}=305.6\times 10^{-23}\penalty\ \text{g}$ (i.e. $1839\penalty\ \text{g/mol}$). The same damping factor used for the ions is employed for the nanoparticles in the Langevin thermostat. Therefore the friction coefficient of each nanoparticles is $\zeta_{p}=305.6\times 10^{-9}\penalty\ \text{N s}/\text{m}$. The nanoparticle sites as well as the charged wall sites interact with the ions and with each other through Coulomb’s law. All wall sites, charged and uncharged, as well as the nanoparticle sites also interact with the ions and with each other through the repulsive Lenard-Jones potential given in eq. 5.2.

The diffusion constant of each ionic species in the bulk can be calculated a priori from the friction coefficients specified above using the Einstein relation, i.e. $\mathcal{D}=k_{\rm B}T/\zeta$. The values that result from that calculation are, $\mathcal{D}_{\rm Na^{+}}^{\rm bulk}\approx 1.1\times 10^{-12}\penalty\ \text{m}^{2}/\text{s}$, $\mathcal{D}_{\rm Cl^{-}}^{\rm bulk}\approx 6.9\times 10^{-13}\penalty\ \text{m}^{2}/\text{s}$ and $\mathcal{D}_{p}^{\rm bulk}\approx 1.3\times 10^{-14}\penalty\ \text{m}^{2}/\text{s}$. However the diffusion coefficients of charged species inside charged nanopores can be significantly different from diffusion coefficients in the bulk Valiskó et al. (2019). Therefore we calculate the diffusion coefficient of anions, cations and nanoparticles inside the confined environment of the nanopores from the mean-squared displacement, $\langle\Delta r^{2}(t)\rangle_{\rm eq}$, of each species inside the nanopores. This calculation is done without a voltage bias applied and therefore $\langle...\rangle_{\rm eq}$ indicates an average at equilibrium. Specifically the diffusion constants are obtained from the relation $\mathcal{D}=\langle\Delta r^{2}(t)\rangle_{\rm eq}/(6t)$. See Figure S7 of the SI for plots of the mean-squared displacements and see Table S2 also in the SI with the diffusion coefficients for each species in each system. The diffusion coefficient for each species does not change significantly from one system to another. The average values are, $\mathcal{D}_{\rm Na^{+}}^{\rm eff}\approx 9.05\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$, $\mathcal{D}_{\rm Cl^{-}}^{\rm eff}\approx 5.85\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$ and $\mathcal{D}_{p}^{\rm eff}\approx 0.24\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$.

The effective diffusion coefficients inside the nanopores observed in the simulations can be compared with the values obtained from a simple hydrodynamic calculation. The simple hydrodynamic calculation assumes that ions and noanoparticles can be treated as spheres. The diffusion coefficients are calculated from the relation $\mathcal{D}=k_{\rm B}T/\zeta$, with $\zeta=6\pi R\eta$ where $\eta$ is the viscosity of the solvent and $R$ is the hydrodynamic radius of each species. For the nanoparticles we consider that the eight sites that compose it are at a distance $0.69$ nm from its center. Also since each site interacts with other species through eq. (5.2) we add $0.4\times 2^{1/6}=0.44$ nm to the previous value. Therefore $R_{p}\approx 1.14$ nm. The soft-core repulsion radii given above in the Ions model subsection are used for the ions. Then, using $\eta\approx 1$ mPas, the diffusion coefficients given by this calculation are $\mathcal{D}_{\rm Na^{+}}^{\rm hyd}\approx 8.7\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$ and $\mathcal{D}_{\rm Cl^{-}}^{\rm hyd}\approx 4.8\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$ and $\mathcal{D}_{p}^{\rm hyd}\approx 1.9\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$. The effective diffusion coefficients for the ions inside the nanopore agree well with the simple hydrodynamic calculation. For the nanoparticles the effective diffusivity in the simulations is about 7.9 times smaller than the one obtained in the simple hydrodynamic calculation. Moreover, below we show that the current rectification results can also be replicated when a smaller friction coefficient is used for the nanoparticles in the simulations.

The electric current that flows trough the nanopore as a function of time, is calculated using a commonly employed method Crozier et al. (2001); Aksimentiev and Schulten (2005); Di Muccio et al. (2021), $I(t)=-\displaystyle\sum_{j=1}^{N}\dfrac{\Delta x_{j}(t)ez_{j}}{L\Delta t}$, where the sum is over all the ions and nanoparticles, $N$, and $ez_{j}$ is the charge of each species. $\Delta x_{j}(t)$ is the displacement in the $x$ direction of charged species $j$ on the time interval $\Delta t$ and $L$ is the total length of the channel. Plots of $I(t)$ during $24\penalty\ \text{ns}$ of our production runs are given in Figure S5 of the SI for all the different systems considered in here. Moreover, in Figure 3 we have reported the time average of $I(t)$ over those $24\penalty\ \text{ns}$ of our production runs for all the systems considered in this work. The statistical uncertainties for the electric currents plotted in Figure 3 and also given in Table S1 in the SI were calculated using the blocking transformation method of Flyvbjerg and Petersen Flyvbjerg and Petersen (1989), specifically Eq. (28) in that reference. The blocking transformation method is applicable to estimating statistical error in averages of time series data and properly accounts for the correlation in the data.

To calculate the total electric potential of the nanopores as a function of length, $\psi(x)$, shown in Figure (4) we use the Poisson equation integrated over the $y$ and $z$ coordinates, i.e. $\dfrac{\psi^{2}(x)}{dx^{2}}=-\dfrac{\rho_{e}(x)}{\epsilon_{0}\epsilon}$. Where $\rho_{e}(x)=N_{\rm A}e\left[z_{c}c_{c}(x)+z_{a}c_{a}(x)+z_{p}c_{p}(x)\right]+\dfrac{\sigma(x)}{h(x)}$ is the charge density as a function of the channel’s length (i.e. averaged over the channel’s width and depth). Where $N_{\rm A}$ is Avogadro’s number, $\sigma(x)$ and $h(x)$ were defined in eq. (5.2). $c_{c}(x),c_{a}(x)$ and $c_{p}(x)$ are molar concentration profiles of cations, anions and nanoparticles respectively. Each concentration profile is obtained by averaging $10^{6}$ simulation frames taken during $24\penalty\ \text{ns}$ of simulation after the systems have reached steady-state. In each frame the simulation box is divided into cubic bins of $(0.8\times 20\times 0.2)\penalty\ \text{nm}^{3}$. The number of particles of each species inside each bin are counted to calculate the concentrations. The steady-state molar concentration profiles are given in Figure 2 and Figures S3 and S4 in the SI. The Poisson equation is solved with boundary conditions, $\psi(x=-40\penalty\ \text{nm})=0$ and $\psi(x=30\penalty\ \text{nm})=-E_{x}/L$.

To check for reproducibility of the rectification results presented here we have performed additional replicas of the simulations for the system that shows the largest rectification ratio. Moreover we have also performed those replicas using a friction coefficient for the nanoparticles about 8 times smaller than the specified in the Nanoparticel model subsection above. More precisely, a friction coefficient of $\zeta_{p}=38.7\times 10^{-9}\penalty\ \text{N s}/\text{m}$ was used for the nanoparticles in the additional replicas. The friction coefficients of the ions were kept the same as those given in the Dynamics subsection above. Snapshots of the othographic projections of the initial conditions and steady-state configurations of these replicas are shown in Figure S8 of the SI. The mean-squared displacement at equilibrium, $\langle\Delta r^{2}(t)\rangle_{\rm eq}$, of each species inside the nanopores for each replica are given in Figure S9. The average values of the diffusion coefficients obtained from these $\langle\Delta r^{2}(t)\rangle_{\rm eq}$ are, $\mathcal{D}_{\rm Na^{+}}^{\rm eff}\approx 9.1\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$, $\mathcal{D}_{\rm Cl^{-}}^{\rm eff}\approx 6.1\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$ and $\mathcal{D}_{p}^{\rm eff}\approx 1.0\times 10^{-10}\penalty\ \text{m}^{2}/\text{s}$. Note that for the nanoparticles this diffusion coefficient is about 4 times larger than the one obtained in the simulations with the larger friction coefficient. Figure S10 and Table S3 in the SI show comparisons of the electric currents and the rectification ratios obtained in the different replicas performed.