Fluids at the Nanoscale: from continuum to sub-continuum transport

The two-centuries-old Navier-Stokes equation is remarkably robust at describing fluid flow down to the smallest scale, typically $\ell_{c}=1~{}\rm nm$ for water in normal pressure and temperature conditions (Bocquet & Charlaix (2010)). This length scale is essentially a lower bound for defining a fluid viscosity $\eta$. Indeed, in macroscopic fluid mechanics, the kinematic viscosity $\nu=\eta/\rho$, where $\rho$ is the mass density, plays the role of a diffusion coefficient for the fluid momentum. For such a diffusion coefficient to be defined, the time required for momentum to diffuse across the system, $\ell_{c}^{2}/\nu$, must be larger than the timescale of molecular motion, which is the microscopic origin of diffusion. A water molecule at a thermal agitation speed of $300~{}\rm m\cdot s^{-1}$ moves by its own size in $\tau_{c}=10^{-12}~{}\rm s$, which defines a molecular time scale. Therefore, viscosity may be defined down to a system size

Below this length scale, water structuring due to surfaces, memory effects and other sub-continuum phenomena come into play: these will be discussed in section 4. For water flow at $10~{}\rm nm$ length scales, the Reynolds number remains smaller than 0.1 up to fluid velocities of $10~{}\rm m\cdot s^{-1}$. Hence, in nanofluidic systems, inertial effects may be safely neglected, and the fluid flow is described by the Stokes equation:

where $p$ is the pressure and $f$ a body force, which may be due, for example, to the application of an electric field (see section 3.3).

Stokes flow is often solved with no-slip boundary conditions: the velocity of the liquid is assumed to vanish at a solid-liquid interface. This is, however, a limiting case of the more general Navier partial slip boundary condition, which enforces that the viscous stress at the interface should be balancing the solid liquid friction force. Within linear response theory, the friction force is proportional to the liquid velocity. For a fluid flowing in the direction $x$ along a surface of normal $z$, the force balance per unit area writes $\sigma_{xz}=\lambda v_{x},$ with $[\sigma]$ the stress tensor and $\lambda$ the friction coefficient per unit area (expressed in $\rm N\cdot s\cdot m^{-3}$). For a Newtonian fluid, $\sigma_{xz}=\eta\partial_{z}v_{x}$, which allows to rewrite the Navier boundary condition as

introducing the slip length $b=\eta/\lambda$. The slip length can be geometrically interpreted as the depth inside the solid where the linearly extrapolated fluid velocity profile vanishes. Accordingly, the no-slip boundary condition corresponds to $\lambda\to\infty$ or $b\to 0$. The effect of the partial slip condition is to simply shift the no-slip velocity profile by the slip velocity, which is not negligible roughly within a slip length from the wall. Since slip lengths up to tens of nanometers have been measured on atomically flat (and hydrophobic) surfaces, slippage is expected to play a crucial role in nanofluidics, and some of its effects will be discussed in the following sections. In the smallest channels, of size $R\ll b$, a perfect slip boundary condition may even be appropriate: the flow is then controlled by entrance effects.

Experimentally, the flow profile inside a nanofluidic channel can hardly be resolved, and one typically measures the total flow rate $Q$. Under a pressure drop $\Delta P$ and no-slip boundary conditions, the flow rate through a cylindrical channel of radius $R$ is given by the Hagen-Poiseuille formula:

This formula assumes a channel length $L\gg R$, and thereby neglects the effect of channel mouths on the flow rate. But the transition from a macroscopic reservoir to a nanoscale channel is a source of viscous dissipation, as the streamlines need to be bent in order for the fluid to enter the channel. These entrance effects may be examined by considering the flow through an infinitely thin nanopore, which is of interest in itself, given the geometry of certain nanofluidic devices (see section 2.1). This problem was addressed by Sampson (1891). For a nanopore of radius $R$ (and vanishing length) under pressure drop $\Delta P$, Sampson obtained the expression of the flow rate as

The scaling in Sampson’s formula naturally emerges from a Stokes equation where the only lengthscale is $R$: $\eta\Delta v=\nabla p\Rightarrow\eta v/R^{2}\sim\Delta P/R$, and the typical fluid velocity is $v\sim Q/R^{2}$. In order to estimate the flow rate through a channel taking into account entrance effects, one may simply add the hydrodynamic resistances of the pore (${\cal{R}}_{p}$) and the channel (${\cal{R}}_{c}$). If one writes $Q_{c}=(\Delta P)_{c}/{\cal{R}}_{c}$ and $Q_{p}=(\Delta P)_{p}/{\cal{R}}_{p}$, then the entrance-corrected flow rate $Q_{pc}$ is obtained by imposing $Q_{p}=Q_{c}=Q_{pc}$ and $\Delta P=(\Delta P)_{c}+(\Delta P)_{p}$, so that

An exact computation (Dagan et al. (1982)) shows that the error made by this a priori crude approximation is less than 1%. Equation 6 makes a continuous transition between the nanopore and nanochannel regimes, and shows that entrance effects are apparently negligible for channel lengths that exceed a few channel radii.

However, the above discussion has crucially not taken into account slippage, which, as we have highlighted in the previous section, is a strong effect at the nanoscale. Introducing a non-zero slip length $b$, the flow rate though a channel becomes

while the flow rate through a pore is not significantly affected (Gravelle et al. (2013)), since the source of dissipation in that case is mostly geometric. A full expression can be obtained by gathering previous results, but in the limit where $b\gg R$, the entrance-corrected flow rate becomes

Thus, the hydrodynamic resistance is actually dominated by entrance effects as long as the channel is shorter than the slip length, rather than the channel radius (see figure 3a). In the presence of significant slippage, one should check whether the low Reynolds number assumption still holds. The average velocity through a channel of radius $R=5~{}\rm nm$ and length $L=1~{}\rm\mu m$, with slip length $b=30~{}\rm nm$, under a pressure drop $\Delta P=1~{}\rm bar$ is $v=8~{}\rm mm\cdot s^{-1}$, which is 25 times faster than the no-slip result, but still well below the $1~{}\rm m\cdot s^{-1}$ threshold established in 3.1.1.

Entrance-effect-dominated transport is particularly relevant for biological nanochannels, due to their relatively small aspect ratio. A striking example is aquaporin, which was recently studied from a hydrodynamic point of view by Gravelle et al. (2013). Aquaporins are channel proteins that selectively transport water across the cell membrane. A simplified geometrical model for the aquaporin consists of two conical vestibules connected by a subnanometric channel where water flows in single file (figure 3b). The single-file transport is expected to be nearly frictionless, and therefore the limit to the aquaporin’s permeability is set by the entrance effects in the conical vestibules. Moreover, since the channel in question is less than 1 nm wide and slip lengths of the order of 10 nm are expected, perfect slip boundary conditions are relevant for the flow, and the viscous dissipation has a purely geometric origin in the curvature of the stream lines. With the notations of figure 3b, Gravelle et al. express the total hydrodynamic resistance of the aquaporin as

Both terms are the analogue of Sampson’s formula (eq. 5). $C_{\infty}=3.75$ replaces the factor 3 in the case of perfect slip boundary conditions, and the factor $\sin\alpha$ appears because at the cone-cylinder transition, the stream lines turn by an angle $\alpha$, as opposed to $\pi/2$ when the transition is from an infinite reservoir to a pore. With that, the permeability $K=R_{\rm AQP}^{-1}$ turns out to be a non-monotonous function on $\alpha$ (figure 3c), with the maximal permeability reached for cone angles in the range measured on aquaporin structures obtained by X-ray crystallography. It appears, therefore, that the geometry of the aquaporin is an optimum for hydrodynamic entrance effects, highlighting their particular relevance for nanoscale fluid transport.

Consider an aqueous solution of monovalent salt. Let $\rho_{+},\rho_{-}$ be the densities of positive and negative ions, respectively; $D$ the diffusion coefficient, here assumed to be the same for ions of either sign, and $\phi$ the electrostatic potential. In a mean-field treatment, the convective-diffusive dynamics of ions are described by a Smoluchowski equation:

where $e$ is the unit charge, $\phi$ the electrostatic potential and $v$ is the fluid velocity field. The mean-field assumption implies in particular that correlations between the ions can be neglected: the potential importance of such correlations in nanofluidics will be discussed in section 4. Until then, we proceed by specifying the electrostatic potential through Poisson’s equation,

where $\epsilon$ is the dielectric permittivity of water. For now we assume it to be isotropic, though this assumption may break down for nano-confined water, as we discuss in section 4. Lastly, we specify the flow velocity through the Stokes equation, which now includes an electrostatic term:

We now apply these three coupled equations to a specific geometry, though the discussion that follows could be generalised to channels of any shpae. For simplicity, we consider a slit-like channel of height $h$, width $w$ and length $L$, with $w,L\gg h$, connecting two reservoirs of salt solution at concentration $\rho_{s}$, extending along the direction $x$, between $z=0$ and $z=h$. When considering ion transport, it is important to note that most surfaces are charged in water, due either to the dissociation of surface groups or to the adsorption of ions (Perram et al. (1973), Grosjean et al. (2019), Mouhat et al. (2020)). We hence assume the channel wall carries a surface charge density $-\Sigma e$ ($\Sigma$ is expressed in elementary charges per unit area, and we assume here the surface charge to be negative).

We first neglect the coupling of ion transport to water transport, and consider the electrophoretic (EP) contribution to the ionic current under an applied electric field $E$: this means that we start by setting the fluid velocity $v=0$. In the steady state, the Smoluchowski equation (15) reduces to the so-called Nernst-Planck equations for the constant ionic fluxes (along the $x$ direction):

Together with the Poisson equation (16), these constitute the widely used Poisson-Nernst-Planck (PNP) framework. In our geometry, the condition $L\gg h$ ensures that in the middle of the channel $\nabla_{x}\rho_{\pm}=0$; moreover to first order in $E$, $\nabla_{x}\phi=-E$ and the densities reduce to their equilibrium values. Hence, the EP contribution to the ionic current writes

In order to compute $I_{\rm ep}$, one needs to find the equilibrium solution of the coupled PNP equations for $\rho_{+}$ and $\rho_{-}$. At equilibrium $j_{\pm}=0$ and the Nernst-Planck equations (18) can be integrated, imposing that in the reservoirs $\phi=0$ and $\rho_{\pm}=\rho_{s}$. This yields a Boltzmann distribution for the ions in the electrostatic potential

introducing a dimensionless potential $\psi$. Combining this with the Poisson equation (16) yields the Poisson-Boltzmann (PB) equation:

which introduces the Debye length $\lambda_{D}=(8\pi\rho_{s}\ell_{B})^{-1/2}$, with $\ell_{B}=e^{2}/(4\pi\epsilon k_{B}T)$ the Bjerrum length. In our geometry, the PB equation has an implicit solution in terms of an elliptic integral (Levine et al. (1975), Andelman (1995)). We will not exploit it here, however, and we will instead recover the relevant limiting behaviours from qualitative considerations.

It is well known that, roughly speaking, the Debye length sets the extension of the diffuse layer of counterions next to a charged surface (Israelachvili (2011)). Hence, if the channel height $h\gg\lambda_{D}$, its two opposing walls do not ’see’ each other. We expect the conductance to be the sum of a bulk term, and a surface term originating in the two Debye layers:

where $\rho_{s}h$ and $\Sigma$ account for the number of charge carriers in the bulk and at surfaces, respectively. In the opposite limit where $h\ll\lambda_{D}$, there is no more distinction between surface and bulk. All the quantities may be considered uniform across the channel: this is called the Debye overlap regime. However, one may not assume that the channel contains only counterions, and one should go back to the thermodynamic equilibrium with the reservoirs, which in this case bears the name of Donnan equilibrium (Bocquet & Charlaix (2010)). One has $\rho_{\pm}=\rho_{s}e^{\mp\psi}$, which implies a chemical equilibrium $\rho_{+}\rho_{-}=\rho_{s}^{2}$ in the channel. Going further, in the limit of long channel length, there should be local electroneutrality: $h(\rho_{+}-\rho_{-})=2\Sigma$. This yields

and the current-voltage relation in the Debye overlap is

Equation 24 displays the first peculiarity of small channels: one may not simply add the surface and bulk contributions. Table 1 (see appendix) lists the values of Debye length for different electrolyte concentrations, showing that the Debye overlap regime is indeed relevant for experimentally accessible nanofluidic systems (see section 2). Qualitatively, eqs. (22) and (24) both predict saturation of the conductance at low salt concentrations at a value determined by the surface charge. The saturation occurs when $\rho_{s}\sim\Sigma/h$, which can be recast in the form $h\sim\rho_{s}/\Sigma\equiv\ell_{Du}$. $\ell_{Du}$ is called the Dukhin length and quantifies the competition between bulk and surface contributions to the conductance. For a channel narrower than $\ell_{Du}$, surface contributions dominate, and vice versa. The Dukhin length is going to be important in our upcoming discussion of entrance effects.

[] \entryBjerrum length$\ell_{B}=\frac{e^{2}}{4\pi\epsilon_{0}k_{B}T}$. Distance between two unit charges at which their interaction energy is $k_{B}T$. \entryDebye length$\lambda_{D}=(8\pi\rho_{s}\ell_{B})^{-1/2}$. Thickness of the diffuse layer of counterions next to a charged surface. \entryDukhin length$\ell_{Du}=\Sigma/\rho_{s}$. Channel width below which surface conductance dominates over bulk conductance. \entryGouy-Chapman length$\ell_{GC}=(2\pi\Sigma\ell_{B})^{-1}$. Distance a unit charge must travel from a charged surface so that its electrostatic energy is reduced by $k_{B}T$.

At this point, a remark should be made concerning the range of validity of equations (22) and (24). Indeed, they have been derived from qualitative considerations, without reference to the exact solution of the PB equation. Now, from eqs. (19) and (20), one obtains more generally

where $\mathcal{E}=(\epsilon/2)\int_{0}^{h}(\partial_{z}\phi)^{2}\mathrm{d}z$ is the electrostatic energy per unit area. The electric double layer can be pictured as a capacitor with charge $\Sigma$, hence one would expect its electrostatic energy to scale as $\Sigma^{2}$; moreover, this is the prediction of the linearised PB equation, i.e. eq. (21) with the approximation $\sinh\psi\approx\psi$. This is in contrast to, eq. (22), which predicts a linear scaling of the conductance with $\Sigma$: this scaling must therefore come from the non-linearities of the PB equation. The PB equation may be linearised if the potential varies by less than $k_{B}T$ across the Debye layer (or across the channel if there is Debye overlap). This is the case in the high concentration/low surface charge limit and specifically when the Debye length is smaller than the so-called Gouy-Chapman length: $\ell_{GC}=(2\pi\Sigma\ell_{B})^{-1}$. On the other hand, Eq. (24), valid for the Debye overlap regime, is safe from a condition on $\ell_{GC}$, as it predicts both quadratic and linear scalings depending on the value of $\Sigma$, as long as there is Debye overlap.

Similarly to what we have discussed for liquid and gas transport, we may now ask, for ion transport, what is the additional electrical resistance due to the channel - reservoir interface. Equivalently, we may want to estimate the ionic conductance of a nanopore of small aspect ratio, say of radius $R$ and length $L\sim R$. The problem has first been considered in the context of biological channels by Hille (1968) and Hall (1975). Hall solved the electrostatic problem with an electrode at infinity and an equipotential disk accounting for the entrance of the pore. Translating the solution into scaling arguments, the current through the pore entrance is $I\sim\pi R^{2}\kappa_{b}\Delta V_{\rm out}/R$, where $\Delta V_{\rm out}$ is the voltage drop at the entrance of the pore, which is expected to occur over a distance $R$, and not over the macroscopic distance between the electrodes. This defines an electrical access resistance as the ratio $\Delta V_{\rm out}/I$. We now specialise to the thin Debye layer regime $\lambda_{D}\ll R$. If one simply sums the access resistance and the channel resistance as given by eq. (22) (adapted to cylindrical geometry), the current-voltage relation becomes

with $\alpha$ a geometric factor which is 2 in Hall’s computation, and introducing the bulk conductivity $\kappa_{b}=2e^{2}D\rho_{s}/(k_{B}T)$. The above equation predicts vanishing conductance as $\rho_{s}\to 0$, since the access resistance becomes infinite in this limit; however, this is not what is observed experimentally in short nanopores (Lee et al. (2012), Feng et al. (2016a)). This inconsistency arises because, for a surface-charged pore, the access current has a surface contribution, in addition to the bulk contribution. This surface contribution may be evaluated starting from charge conservation at the surface, which imposes a relation between the axial and radial components of the electric field, as pointed out by Khair & Squires (2008):

with $\kappa_{s}$ the surface conductance and $\theta$ the Heaviside function, accounting for the discontinuity of surface charge at the pore boundary, which leads to subtle consequences. Indeed, it reveals that the Dukhin length, $\ell_{Du}=\Sigma/\rho_{s}\sim\kappa_{s}/\kappa_{b}$, is the relevant lengthscale for the surface contribution to the electric field outside the pore, instead of the channel radius or Debye length. The Dukhin length appears here as an electrostatic healing length: feeding the surface conduction at the nanopore mouth disturbs the electric field lines in the bulk over a length $\ell_{Du}$. This interpretation is supported by the numerical results of Lee et al. (2012) as shown in figure 5a. The bending of the electric field lines can be estimated by a perturbative approach (Lee et al. (2012)) that leads to an analytical expression for the corrected current-voltage relation:

The surface-charged pore therefore appears, from the perspective of entrance effects, as an uncharged pore of effective size $R+\ell_{Du}$, since the geometrical prefactor $\beta\approx 2$. In the limit of vanishing salt concentration, $\rho_{s}\to 0$, the conductance indeed saturates, and, as expected, the entrance correction disappears in the limit of large aspect ratio pores. Lee et al. successfully compared the prediction of eq. (28) to experimental measurements in SiN nanopores in the range $100-500~{}\rm nm$. More recently, it has been used to describe conductance measurements in $\rm MoS_{2}$ nanopores down to 2 nm in diameter (Feng et al. (2016a)).

So far, we have neglected any coupling of ion transport to fluid transport. However, in the presence of charged surfaces, an external electric field exerts a net force on the charged Debye layer, which sets the fluid in motion. This interfacially-driven flow is termed electro-osmotic (EO) flow. The EO flow in turn drags along the ions in the Debye layer, which makes a supplementary contribution to the ionic current, that we denote $I_{\rm eo}$. This contribution has a convenient exact expression in the slit-like geometry considered in section 3.3.2. In analogy to eq. (19),

From the Poisson equation (16), we may replace $e(\rho_{+}-\rho_{-})=-\epsilon\partial^{2}_{z}\phi$. Then, integrating by parts,

Now, we may use the partial slip boundary condition $v(0)=b\partial_{z}v|_{0}$, as well as the electrostatic boundary condition $\partial_{z}\phi|_{0}=e\Sigma/\epsilon$. This stems from the Gauss theorem applied to the surface, and the assumption that the medium outside the channel has much lower dielectric permittivity than water (Andelman (1995)). Moreover, an integration of the Stokes equation (17) allows the replacement $\partial_{z}v=-(\epsilon E/\eta)\partial_{z}\phi$. Altogether, we obtain

where $\mathcal{E}=(\epsilon/2)\int_{0}^{h}(\partial_{z}\phi)^{2}\mathrm{d}z$ is the electrostatic energy per unit surface introduced in 3.3.1. It appears here that slippage has an additive contribution to the EO current, which strikingly does not depend on any electrolyte properties.

We now consider the result in two limiting cases. First, in the thin Debye layer limit, by comparing with eq. (25), one may identify

where $I_{\rm ep}^{\rm surf}$ is the surface contribution to the electrophoretic current. To quantify the importance of the EO contribution, one may compute the analogue of a Peclet number:

We have introduced here the ion diameter $d_{i}$ by using the Einstein relation $D=k_{B}T/(3\pi\eta d_{i})$ and we recall that $\ell_{GC}=(2\pi\Sigma\ell_{B})^{-1}$. In water at room temperature, $d_{i}\sim\ell_{B}\sim 0.7~{}\rm nm$, therefore, in the absence of slippage, the EO contribution is of the same order as the surface EP contribution, $Pe\sim 1$. Table in the supplemental appendix lists, for reference, values of the Gouy-Chapman length for typical surface charge values. As these are generally in the nanometre range, even in the case of moderate slippage ($b\sim 10~{}\rm nm$), there is a strong enhancement of the EO contribution. The threshold confinement below which the resulting surface contribution dominates over the bulk contribution is given by a rescaled Dukhin length: $\ell_{Du}^{*}=(b/\ell_{GC})\ell_{Du}$.

Second, in the Debye overlap regime, the EO current is readily determined from eq. (29), since the ion densities may then be considered uniform across the channel and are given by eq. (23). Uniform ionic densities also imply that the flow has no longer a surface, but rather a volume driving. It is simply a Poiseuille flow, with the pressure gradient $\Delta P/L$ replaced by the electric driving force $e(\rho_{+}-\rho_{-})E$. Altogether one obtains

Since in Debye overlap it makes no more sense to distinguish a surface and a bulk contribution, we define the Peclet number as the ratio of the EO current to the total EP current:

Figure 5d illustrates the dependence of this Peclet number on $h$ and $\ell_{GC}$, with $\rho_{s}=10^{-2}~{}\rm M$ so that there is Debye overlap. In the absence of slippage, the EO contribution dominates only for high surface charges. However, if a small slip length $b=5~{}\rm nm$ is introduced, the Peclet number exceeds 1 for all reasonable surface charges. Indeed, the Peclet number is then essentially determined by the ratio $b/\ell_{GC}$. This highlights that in the Debye overlap regime, ionic conduction should be mostly driven by electro-osmosis.

If strong EO flows are expected, then the reciprocal effect, streaming current, is expected to be large as well. The streaming current results from the application of a pressure gradient $\Delta P/L$: in Debye overlap, the charge density $2\Sigma/h$ is simply dragged along by the Poiseuille flow. The current reads

where we have phenomenologically defined the zeta potential ($\zeta$) via the streaming mobility. Independent measurements of the surface charge through the voltage-driven current and of the zeta potential from the pressure-driven current may allow in principle to estimate the slip length $b$. Such an estimate does not replace a direct measurement, however, in particular because the slip length may directly depend on surface charge (Xie et al. (2020)); and additional charge-surface coupling effects may also occur, as we discuss below.

A remark is in order at this point concerning the nature of surface charge, that we have so far assumed to be constant, whatever the conditions. Usually, surface charge is considered to result from the acid base reactivity of the surface when dipped into water, of the type $\rm[AH]_{s}\to[A^{-}]_{s}+H^{+}$, where the negative group $\rm[A^{-}]_{s}$ remains fixed on the surface, while the proton diffuses in solution. In the case of the air-water interface, surface charge may result from the adsorption of hydronium ions (Mamatkulov et al. (2017)), or charged impurities. The adsorption of surfactant impurities was also proposed as a charging mechanism for general hydrophobic surfaces (Uematsu et al. (2019)). Graphite and hBN, which are of particular relevance in nanofluidics, have a priori no obvious acid-base reactivity in water, but they may develop a surface charge, through chemisorption or physisorption of hydroxyde ions, as inferred from experiments (Secchi et al. (2016b), Siria et al. (2013)), and recently confirmed by ab initio simulations (Grosjean et al. (2019)). Such a picture of surface charge implies that it may actually depend on electrolyte concentration, which is an example of charge regulation. Qualitatively, the salt concentration affects the surface potential, and therefore the concentration of $\rm H^{+}$ and $\rm OH^{-}$ ions at the surface, which in turn affects the chemical (or the adsorption) equilibrium which governs the surface charge. Such charge regulation has been invoked to explain several experiments on carbon nanotubes, where a scaling $I\propto\rho_{s}^{1/2}$ or $I\propto\rho_{s}^{1/3}$ was observed (Secchi et al. (2016b), Liu et al. (2010), Pang et al. (2011)). Various models have been developed (Secchi et al. (2016b), Biesheuvel & Bazant (2016), Manghi et al. (2018), Uematsu et al. (2018)), predicting a rich panorama of different sublinear scalings of conductance with salt concentration, depending on the conditions.

A further step that may be required to accurately describe the surface charge is to take into account its mobility. Surfaces charges may be mobile when embedded in lipid bilayers, or with more relevance to nanofluidics, when resulting from adsorbed ions. For instance, the simulations of Grosjean et al. showed that physisorbed hydroxyde ions on graphene surfaces retain a high lateral mobility (Grosjean et al. (2019)). The effect of a mobile surface charge on ion and fluid transport coefficients (in the thin Debye layer regime) was the subject of several investigations (Maduar et al. (2015), Silkina et al. (2019), Mouterde & Bocquet (2018)). Using the framework of Mouterde & Bocquet (2018), one may introduce two new friction coefficients: $\lambda_{-}$ between the adsorbed (negative) ions and the wall, and $\xi_{-}$ between the adsorbed ions and the fluid, in addition to the water-wall friction coefficient $\lambda$. The force balance on the interfacial fluid layer results in a modified partial slip boundary condition involving the tangential electric field at the surface:

where the effective slip length is

The slip length is therefore reduced by the surface charge mobility, and the boundary condition involves an extra electric term. This has in general a moderating effect on the transport coefficients. For example, in the case of streaming current, the fluid brings along surfaces charges in addition to the counterions, so that the total current is reduced. Notable exceptions are diffusio-osmotic mobility and conductance: the overall ionic conductance increases when surface charges are able to move in response to the electric field.

We have just introduced ion-wall and ion-water friction for adsorbed ions in the interfacial layer. Now, for the smallest accessible channels, which are comparable to the ion size (Radha et al. (2016), Lee et al. (2010)), such friction needs to be taken into account for all the ions, and should appear at the level of the transport equations. This was proposed, for instance, to rationalise the non-linear voltage-pressure couplings observed in angstrom-scale slits by Mouterde et al. (2019). They observed that the streaming mobility $\mu$, defined by $I_{\rm str}=hw\mu\Delta P/L$, depends on the applied voltage in a qualitatively different way for graphite and BN slits. In the proposed ’surface PNP’ model, the ion-wall ($\lambda_{\pm}$) and ion-water ($\xi_{\pm}$) friction coefficients are introduced in the same way as above, except they no longer apply to adsorbed ’surface charge’ ions, but to all the regular salt ions. The equivalent of the Smoluchowski and Stokes equations need then to be rederived by considering the force balance on a single ion and on an element of fluid, respectively, as shown in the Appendix.

There is a rich phenomenology to be harnessed from the new couplings that appear when all components of the nanofluidic system interact with the surfaces. Nevertheless, such surface PNP theory is still in its infancy, and it requires input from more microscopic modelling in order to estimate the various friction coefficients it introduces. It even raises mathematical difficulties in the writing of a one-dimensional Poisson equation, whose Green’s function has infinite range. Finally, the relevance of this phenomenological theory should be assessed in light of the correlation and structuring effects – specific to the smallest scales – that will be described in the next section.

The relative particle number fluctuations in an open system of $N$ particles scale as $1/\sqrt{N}$. In a nanopore of radius $1~{}\rm nm$ and length $10~{}\rm nm$, there are on average about 1000 water molecules and these fluctuations are of a few percent and already not negligible. Fundamentally, the description of such fluctuations requires a statistical mechanics framework, such as density functional theory (Barrat & Hansen (2003)). However, these approaches remain quite formal and difficult to implement, particularly for finite length pores which are non translationally-invariant systems.

A simpler coarse-grained approach to fluctuations – the so-called fluctuating hydrodynamics (FH) – has been introduced by Landau & Lifshitz (1977). In fluctuating hydrodynamics, a random stress tensor $S$ is added to the Navier-Stokes equation:

where $S$ satisfies a fluctuation-dissipation theorem: $\langle S_{ij}(\mathbf{r},t)S_{kl}(\mathbf{r}^{\prime},t^{\prime})\rangle=2\eta k_{B}T(\delta_{ij}\delta_{kl}+\delta_{il}\delta_{jk})\delta(\mathbf{r}-\mathbf{r}^{\prime})\delta(t-t^{\prime})$, similarly to the fluctuating force in a Langevin equation; $\rho$ is the mass density of the fluid. This description introduces thermal fluctuations in a setting where, otherwise, a continuum hydrodynamic description holds. FH has been extensively discussed in the literature (Hauge & Martin-Lof (1973), Fox & Uhlenbeck (1970), Mashiyama & Mori (1978)) and we only consider here its basic implications for fluid transport in a nanopore.

In standard hydrodynamics, the fluid inside a nanopore is at rest when subject to no external force. The effect of fluctuations is to induce a stochastic flow through the nanopore, which can be pictured as a stochastic motion of the fluid centre of mass (CM). Detcheverry & Bocquet (2012) solved the FH equations for the velocity correlation function, and found that the CM motion was described by a non-Markovian Langevin equation of the type

where the random force $F(t)$ and the memory kernel $\xi$ are related by $\langle F(t)F(0)\rangle=2k_{B}T\xi(t)$. The memory kernel is indeed non-trivial (i.e. not reduced to a $\delta$ function), as confirmed by molecular dynamics (MD) simulations (Detcheverry & Bocquet (2013)). Memory effects in the diffusion of fluid inside a nanopore are indeed expected from the analysis of relevant timescales. A nanopore of radius $R$ and length $L$ contains a mass $m=\rho\pi R^{2}L$ of fluid, and its friction coefficient on the wall can be calculated from the Navier-Stokes equation as $\xi_{s}\approx 8\pi\eta L$. The velocity of the fluid CM then relaxes on a timescale $\tau\sim m/\xi_{s}\sim R^{2}/\nu$. But $R^{2}/\nu$ is the time required for fluid momentum to diffuse across the nanopore, that is the time for the wall friction force to establish itself. Therefore, the friction force may not adapt instantaneously to the CM velocity, and memory effects are to be expected.

Properties that are determined by short timescale dynamics may be affected by those memory effects as it was shown to be the case for solute mobility in confinement (Daldrop et al. (2017)). More generally, memory effects play an important role in barrier-crossing processes (Kappler et al. (2018)). But other properties only depend on the long timescale dynamics. In particular, the fluid centre of mass diffusion coefficient, given by $\mathcal{D}=k_{B}T/(2\int_{0}^{\infty}\xi(t)\mathrm{d}t)$, reduces to the Einstein expression $\mathcal{D}=k_{B}T/\xi_{s}$, that would be expected if the diffusion was Markovian. Similarly, the hydrodynamic slip length $b$ is found to satisfy a Green-Kubo relation whatever the form of the memory kernel $\xi(t)$ (Bocquet & Barrat (2013)).

From a practical perspective, one may evaluate the diffusion coefficient $\cal D$. From the Einstein relation, and the assumption of Poiseuille flow with slip length $b$, one obtains

where, in the most general case, the friction coefficient $\xi$ has two contributions, one from the channel interior and one from the entrance effects (Detcheverry & Bocquet (2012)). Overall, the fluid CM diffusion $\mathcal{D}$ can be seen as a supplementary contribution to the diffusion coefficient of a particle inside the nanopore. If the size of the particle is comparable to that of the nanopore, its self-diffusion is strongly hindered (Renkin (1954)) and the fluid contribution may actually dominate the particle diffusion (Detcheverry & Bocquet (2012)).

Thermal fluctuations are at the origin of the noise in ionic current measurements through a nanopore. At zero applied voltage, the total noise amplitude satisfies the Nyquist relation $\langle I^{2}\rangle=2k_{B}TG$, where $G$ is the nanopore conductance. The electro-osmotic term in the conductance makes a supplementary contribution to the noise, which can actually be traced back to the hydrodynamic fluctuations of the fluid CM described in the previous paragraph (Detcheverry & Bocquet (2012)).

When resolved in frequency, numerous experiments on artificial nanochannels (Smeets et al. (2008), Siwy & Fuliński (2002), Hoogerheide et al. (2009), Tasserit et al. (2010), Secchi et al. (2016b)) and biological pores (Wohnsland & Benz (1997), Bezrukov & Winterhalter (2000)), have shown that the current spectrum exhibits ’pink noise’, that is noise that scales with frequency $f$ as $f^{-\alpha}$, with $\alpha$ in the range $0.5-1.5$. Such a spectrum is traditionally described by the empirical Hooge’s law (Hooge (1969)):

with $A$ a coefficient inversely proportional to the number of charge carriers. Notably, in experiments, the $f^{-\alpha}$ scaling is observed down to frequencies below 1 Hz, showing that some correlations in the system under study exist even at such low frequencies. Correlations between ions have been put forward as a possible origin for the pink noise (Zorkot et al. (2016)). On the other hand, Gravelle et al. (2019) demonstrated in MD simulations that pink noise persisted even with non-interacting ions, if reversible ion adsorption was allowed on the channel wall. For simple diffusive dynamics, the longest correlation time that may be expected for adsorption-desorption processes is $R^{2}/D$, with $R$ the pore radius and $D$ the ion diffusion coefficient. For a 10 nm pore, this would correspond to a frequency cutoff of 100 MHz, below which the noise spectrum should be flat. Such a high cutoff frequency is in clear contradiction with experiment. However, Gravelle et al. showed that in the presence of reservoirs, excursions of ions outside the pore coupled to adsorption-desorption processes on the pore surface result in ionic correlations on much longer time scales, hence predicting pink noise down to very low frequencies.

So far, we have only considered thermal fluctuations of the fluid itself. However, the confining walls might also be subject to fluctuations. Wall fluctuations are known to enhance fluid self-diffusion by inducing flows, through the Taylor-Aris mechanism (Taylor (1953), Aris (1956)). On the other hand, static constrictions slow down diffusion in a channel through Fick-Jacobs entropic trapping (Malgaretti et al. (2014), Reguera & Rubí (2001)). Recently, MD simulations have shown diffusion enhancement due to the propagation of phonons in carbon nanotubes (Ma et al. (2015), Cruz-Chú et al. (2017)).

On the theoretical side, a very general framework for evaluating the self-diffusion coefficient for a fluid confined in a fluctuating channel was developed by Marbach et al. (2018). In a slit like-channel of height $h$, the fluctuation-corrected diffusion coefficient $D$ is formally related to the spectrum of wall fluctuation $S(k,\omega)$ by

where $D_{0}$ is the bare diffusion coefficient. Equation (43) has two limiting regimes, that are governed by a ’fluctuation-related’ Peclet number, which may be defined as $Pe=\ell^{2}/(D_{0}\tau)$, where $\ell$ and $\tau$ are the characteristic length and time scales of the fluctuations. The diffusion coefficient is enhanced for $Pe>1$ ($D=D_{0}(1+3\langle\delta h^{2}\rangle/h^{2})$), and reduced for $Pe<1$ ($D=D_{0}(1-\langle\delta h^{2}\rangle/h^{2})$) and both tendencies are observed experimentally (Marbach et al. (2018)). Equation (43) thus bridges the limiting Taylor-Aris and Fick-Jacobs results.

The permeability of a fluctuating channel is also expected to be impacted by wall fluctuations. However, permeability is related to the diffusion of the fluid centre of mass (see eq. (41)), which is different from the fluid self-diffusion; and extending the above framework to centre of mass diffusion remains to be done.

While a fluid appears disordered at hydrodynamic lengthscales, its molecular nature manifests itself under confinement in the form of structuring effects. Even at a non-confined solid-liquid interface, the attractive interactions between the solid and the liquid result in molecular layering, as verified experimentally by x-ray spectroscopy (Cheng et al. (2001)). In between two surfaces, a liquid adopts a layered structure even in the absence of interactions, simply due to geometrical constraints (Israelachvili & Pashley (1983)). The onset of structuring clearly represents the transition from continuum to sub-continuum transport.

The threshold confinement at which this transition occurs has been studied for both planar and cylindrical geometries using MD simulations. For water in carbon nanotubes, Thomas & McGaughey (2009) have found that it retains a bulk-like disordered structure down to a tube diameter of 1.39 nm ((10,10) tube). In a (9,9) tube (1.25 nm) water was found to structure in stacked hexagons (see figure 7), and in a single file chain in a 0.83 nm (6,6) tube. For water between two graphene sheets, four distinct layers could be observed at a 1.35 nm distance between sheets (Neek-Amal et al. (2016)), and a single monolayer of water was observed for confinement below 0.8 nm (Neek-Amal et al. (2018)). In should further be noted that the observation of an exotic square ice phase for few-layer water between two graphene sheets has been claimed (Algara-Siller et al. (2015)).

In the sub continuum-regime, transport properties show strong qualitative deviations from bulk expectations. The MD simulations of Thomas & McGaughey (2009) showed non-monotonous permeability for carbon nanotubes as a function of diameter, and similar observations were made for the capillary filling velocity of carbon nanotubes, studied by Gravelle et al. (2016). On the experimental side, this behaviour can be put in perspective with the recent results by Radha et al. (2016) for capillary-driven flow through slit-like graphene channels where a peak in flow rate at around 1.3 nm confinement was observed.

Thermodynamically, the effect of molecular fluid structuring may be described by a supplementary contribution to the pressure, the so-called disjoining pressure. For a fluid confined between two walls of area $\mathcal{A}$ separated by a distance $h$, in equilibrium with a reservoir at pressure $P_{0}$, it is defined in full generality as

where $G$ is the Gibbs free energy of the whole system (Israelachvili (2011), Barrat & Hansen (2003)).

In principle, $G$ varies with $h$ because of distance-dependent interactions between the surfaces; the disjoining pressure can then be seen as the force the fluid needs to exert in order to oppose those interactions. These are generally of two types. On the one hand, there are van der Waals interactions that scale algebraically as $h^{-2}$, yielding a contribution $\Pi_{vdW}=-{A_{H}\over 6\pi h^{3}}$ to the disjoining pressure, with $A_{H}$ the Hamaker constant (although more complex dependencies may occur in specific geometries). On the other hand, electrostatic interactions lead to contributions that scale exponentially with the distance $h$ due to ionic screening. The combination of these two interaction terms is the basis for the DLVO (Derjaguin-Landau-Verwey-Overbeek) theory.

Now, if there is fluid structuring that depends on confinement, it makes a supplementary contribution to the variation of $G$. This contribution exists even for a hard sphere fluid due to excluded volume effects, and is purely entropic in nature. It quantifies essentially how much the structure of the fluid in confinement differs from the one in the bulk. For hard spheres of diameter $\sigma$, the disjoining pressure is usually described by the expression (Israelachvili (2011), Kralchevsky & Denkov (1995)):

where $\rho_{\infty}$ is the bulk density. This is an oscillating function of $h$ that decays to 0 at large $h$, where there is no more fluid structuring. It is bounded by $\Pi_{d}(0)=\rho_{\infty}k_{B}T$, which evaluates to about 1000 bar for water at room temperature. In practice, for fluid flow measurements where a pressure drop is imposed between two reservoirs, this ( huge) contribution to the pressure does not intervene directly. However, in a capillary flow geometry, which involves liquid-vapour interfaces, the disjoining pressure may be much larger than the bare capillary pressure and act as the main driving force. Indeed, the capillary pressure drop across a meniscus of radius $R=1~{}\rm nm$ is typically $2\gamma/R=140~{}\rm bar$ (for water with surface tension $\gamma=72~{}\rm mN\cdot m^{-1}$). Driving by disjoining pressure was observed in simulations of carbon nanotube capillary filling (Gravelle et al. (2016)), where the subtle dependence of the water molecular structure – and therefore of the disjoining pressure – on nanotube radius resulted in the imbibition velocity displaying oscillatory behaviour versus the confinement. Experimentally, a similar phenomenology was observed by Radha et al. (2016), who measured the capillary flow of water across angstrom scale slits and found a non-monotonous dependence of the evaporation rate on confinement. These measurements could be interpreted in terms of disjoining pressure, which was indeed found to be of the order of 1000 bar in MD simulations of the experimental system (Neek-Amal et al. (2018)).

A first understanding of sub-continuum transport may be provided by extending the continuum description with effective, confinement-dependent, values for density, viscosity and slip length. In a layered fluid, the density in each layer is generally higher than in the bulk, and a general trend is that viscosity increases (Schlaich et al. (2017)). More exotically, an oscillating viscosity as a function of confinement was observed in simulations of water between two graphite slabs (Neek-Amal et al. (2016)). Neek-Amal et al. (2018) could reproduce the experimental results of Radha et al. (2016) by using simulated values of disjoining pressure, along with a Poiseuille formula with effective density and viscosity; however, they assumed a viscosity-independent slip length. Now, in general, care should be taken in defining a viscosity for a structured fluid: as its density is non-uniform, a position-dependent viscosity should be introduced in order to describe the details of the flow profile. For confinement below five molecular diameters (in the case of a Lennard-Jones fluid) Travis et al. (1997) have shown that even a position-dependent viscosity is not sufficient and a non-local viscosity kernel should be used (Zhang et al. (2004)). Ultimately, for confinement below a few molecular layers (typically 1 nm for water), the notion of viscosity itself – which is intrinsically a continuum quantity – looses its fundamental meaning.

Nevertheless, a simple picture of sub-continuum fluid transport is possible in the case of large slippage, since surface friction then becomes the main mechanism resisting fluid transport. This appears indeed when looking at the limit $b\gg R$ in the Poiseuille formula (7), which yields for the average flow velocity:

where we have obtained the second equality by relating the slip length to the liquid-wall friction coefficient $\lambda$: $b=\eta/\lambda$; note that, for simplicity, we forget here about entrance. The viscosity does not enter the permeability, as expected for surface dominated friction. The fluid moves indeed in the channel as a single block, and dissipation only intervenes at the surface. The crucial parameter governing the transport is then the solid-liquid friction coefficient $\lambda$. A model for evaluating $\lambda$, tested against MD simulations of water in carbon nanotubes, was proposed by Falk et al. (2010). They obtain

where $q_{0}$ is the wave-vector corresponding to the solid lattice spacing, $f_{q_{0}}$ is the Fourier component of the fluid-solid interaction potential at that wave-vector, and $S(q)$ is the liquid structure factor; $\tau$ is the correlation time of the force between the solid and the liquid, which is found to not depend on confinement. In this way, $\lambda$ is found to be directly related to the liquid structure and to the liquid wall interaction. Such a friction-dominated approach was applied successfully, for example, to predict the permeability of ultra-confined alcanes in a nanoporous matrix, showing a scaling dependence of the permeance on the alcane length (Falk et al. (2015)). However, predicting the solid-liquid friction coefficient remains to large extents an open problem. As we discuss in the following (see 4.2.6), its determination in MD simulations still remains ambiguous quantitatively, and the role of underlying quantum effects should be assessed.

The most extreme type of structuring is single-file arrangement of fluid molecules. In the case of water transport, this is the realm of biological channels (Horner & Pohl (2018b)), where water is conducted through sub-nanometric openings in transmembrane proteins. Single file transport may also occur in artificial channels: it was observed in MD simulations for carbon (Hummer et al. (2001)) and boron nitride (Won & Aluru (2007)) nanotubes of 0.8 nm in diameter, and experimental observation in carbon nanotubes was recently claimed (Tunuguntla et al. (2017)). Single-file transport is a field of its own and there are dedicated reviews on the subject (Horner & Pohl (2018b), Köfinger et al. (2011)); here we will only sketch the main ideas.

The notion of single file transport and the underlying exclusion transport models are not restricted to the study of confined fluids. Indeed, a general feature of particles in one dimension that cannot cross each other is sub-diffusive behaviour: the mean squared displacement of a particle scales with time as $\langle\Delta x^{2}\rangle\propto t^{1/2}$, as opposed to linear scaling for normal diffusion (Levitt (1973)). Such anomalous diffusion is encountered in a variety of fields (Fouad & Gawlinski (2017)), and may be understood in terms of the normal diffusion of vacancies in the 1D chain; in the case of fluid transport, it was studied by Chou (Chou (1999, 1998)). However, this effect is not expected to play a decisive role in single-file transport through short channels, where vacancies are unlikely to occur (Kalra et al. (2003)).

In the case where no vacancies are expected, a simple model for single-file transport was proposed decades ago by Finkelstein (1987) for describing water transport in biological channels. It stems from a global force balance on the water chain, and is as such equivalent to eq. (46), which we have derived as a strongly confined limit of the Poiseuille formula. In Finkelstein’s model, the chain of $N$ water molecules (of total length $L$) moving at velocity $v$ is subject to a friction force $-N\xi v$, and to a pressure driving $Nv_{w}\Delta P/L$, where $v_{w}$ is the volume of one water molecule. The force balance then leads to

which is indeed eq. (46), granted the identification $\xi=(2v_{w}/R)\lambda$. In the biophysics literature, the preferred quantity for characterising a channel’s transport properties is the unitary channel permeability $p_{f}$: $p_{f}$ is the number of water molecules that crosses the channel per unit time, per unit osmolyte concentration difference applied across the channel, therefore expressed in $\rm m^{3}\cdot s^{-1}$. For an osmolyte concentration difference $\Delta\rho_{s}$, the osmotic pressure difference is $k_{B}T\Delta\rho_{s}$, and hence the unitary permeability is

where we have defined $a=L/N$ the average spacing between water molecules. The diffusion coefficient of the water chain centre of mass may be expressed through the Einstein relation as $D=k_{B}T/(N\xi)$, and therefore $p_{f}=v_{w}D/a^{2}$, hence the diffusion coefficient and the unitary permeability may be used interchangeably. A typical value for the permeability in biological channels such, as Aquaporin or Gramicidin A, is in the range of $10^{-14}-10^{-13}~{}\rm cm^{3}\cdot s^{-1}$, which corresponds to 1-10 water molecules crossing the channel per nanosecond under a 1 bar pressure drop.

In an alternative phenomenological description, transport through the channel is viewed as an activated process with activation energy $\Delta G^{\ddagger}$ (Horner & Pohl (2018b)). The permeability is then expressed in the framework of transition state theory as $p_{f}=v_{w}\,\nu_{0}\,\exp[-\Delta G^{\ddagger}/k_{B}T]$, with $\nu_{0}\sim 10^{13}$ s^-1 a molecular attempt frequency. For biological channels, this relation is well verified by independent measurements of $p_{f}$ and $\Delta G^{\ddagger}$, with $\Delta G^{\ddagger}$ of the order of $5~{}\rm kcal\cdot mol^{-1}$. However, values of $\Delta G^{\ddagger}$ remains debated for transport measurements in carbon nanotubes (Horner & Pohl (2018a), Tunuguntla et al. (2018)).

Although it provides a general guiding line, Finkelstein’s formula is challenged both in artificial and biological channels. In the case of protein channels, the assumption of uniformly smooth walls breaks down, since there are discrete hydrogen bonding sites (Horner & Pohl (2018b)). Hence $p_{f}$ is not found to be inversely proportional to channel length (Saparov et al. (2006)): it rather has an exponential dependence on the number of hydrogen bonding sites (Horner et al. (2015)). Such an exponential dependence suggests a collective transport mechanism of the water chain, with bursts requiring the breaking of multiple hydrogen bonds at once. There is, however, a notable disagreement between experiments and simulations, as the simulations of water transport through polyalanine channels showed no dependence of $p_{f}$ on the channel length (Portella et al. (2007)).

In simulations of carbon nanotube channels, all the single-file water molecules were clearly shown to move in a correlated fashion (Hummer et al. (2001)): water transport occurs when all the molecules simultaneously shift by one molecular diameter. These dynamics were successfully described by a continuous time random walk model (Berezhkovskii & Hummer (2002)), or equivalently by diffusion of a collective coordinate of the water molecules (Zhu et al. (2004)), which actually echoes our discussion of fluctuations in a more general setting (eq. 40 and Detcheverry & Bocquet (2012)).

Overall, there is still much to understand about the collective motions and subtle surface interactions involved in single-file transport, and it is an active field of research.

As confined water becomes structured, its response properties to an external electric field are accordingly modified. When an electric field $\bf E$ is applied in water, the individual molecules reorient and polarise, creating an electric field themselves. The total electric field $\bf E$ is then the sum of the polarization field of the water molecules and of the externally applied field $\mathbf{D}/\epsilon_{0}$, where $\bf D$ is called the electric displacement and $\epsilon_{0}$ is the vacuum permittivity. In bulk water, the dielectric response, that is the relation between $\bf D$ and $\bf E$, may be expressed through a single scalar quantity, the relative permittivity $\epsilon\approx 80$: $\mathbf{D}=\epsilon\epsilon_{0}\mathbf{E}$. However, the most general (static) linear response may be anisotropic, space-dependent and non-local:

The relative permittivity is then a tensor with components $\epsilon_{\alpha\beta}(\mathbf{r},\mathbf{r}^{\prime})$. While MD simulations show that the dielectric response in water may be considered local ($\epsilon_{\alpha\beta}(\mathbf{r},\mathbf{r}^{\prime})=\epsilon_{\alpha\beta}(\mathbf{r})\delta(\mathbf{r}-\mathbf{r^{\prime}})$), it becomes anisotropic and space-dependent in the vicinity of interfaces, as a consequence of the water layering (Bonthuis et al. (2012),Bonthuis et al. (2011)). Qualitatively, the orientations of the water dipoles are anti-correlated in the direction perpendicular to the interface, resulting in a reduced permittivity in that direction, while the permittivity is largely unaffected parallel to the interface. In planar confinement, this behaviour could be captured by an effective medium model, in which the water is described by a space-independent, but anisotropic permittivity $(\epsilon_{\parallel},\epsilon_{\perp})$ (the parallel direction is not confined). While $\epsilon_{\parallel}$ essentially retains its bulk value, $\epsilon_{\perp}$ is reduced by up to an order of magnitude for confinements below 1 nm (Schlaich et al. (2016), Zhang et al. (2013)). Such a reduction of the perpendicular dielectric response was recently observed experimentally for water confined between a graphite and a boron nitride crystal (Fumagalli et al. (2018)). A deviation from the bulk value was measured up to nearly 100 nm confinement (figure 7). The results were well described by assuming each interface carried a 7 Å thick layer of very low permittivity ($\epsilon=2.1$) ”electrically dead water”.

An effective medium model based on MD simulations was developed by Loche et al. (2019) for cylindrical confinement of water in carbon nanotubes. Similarly to the case planar case, the radial permittivity is found to be reduced by up to an order of magnitude for tube radii smaller than 3 nm (figure 7). However, for the smallest tubes (below 1 nm radius), the longitudinal permittivity is found to increase with respect to its bulk value, and it skyrockets 1 to 3 orders of magnitude in a 0.4 nm radius tube, where water is in a single file arrangement.

Knowledge of these modified dielectric properties is important when considering ion transport in strong confinement, as will be discussed in the following (see 4.3.1).

We will close this section on liquid structuring with a word of warning. Many of the results we have presented so far have been obtained using classical MD simulations, and like any type of simulation, these come with some underlying assumptions. In particular, as we are dealing with strongly confined fluids where all the molecules interact with the confining surface, the modelling of that surface plays an increasingly important role. Since the solid is represented by rigid balls interacting with classical force fields, no electron dynamics may be described by the simulation. However, some coupling between water and electron dynamics may be expected, based on a lengthscale argument. The maximum of the dielectric response of a solid indeed occurs at a lengthscale $k_{\rm TF}^{-1}$, where $k_{\rm TF}$ is the Thomas-Fermi wave-vector (Mahan (1990)). Typically, $k_{TF}^{-1}\approx 1~{}\rm nm$ in graphite, hence the electrons in graphite could expected to respond to the dynamics of individual water molecules.

There is growing evidence that such electron-water couplings indeed occur. A first consequence of taking into account electron dynamics is that a solid should appear polarisable: in simulations, the fluctuations of water next to a graphite surface were found to be strongly affected by the polarisability of carbon atoms (Misra & Blankschtein (2017)). On the experimental side, the induction of an electronic current by a water flow inside a carbon nanotube, as well as the reverse phenomenon, were observed (Ghosh et al. (2003), Rabinowitz et al. (2020)).

These results call for caution in the interpretation of MD simulations, which, as we have pointed out, lack electron dynamics. They also call for the development of new simulation methods that would take such dynamics into account, while not limited by the very small system size of ab initio simulations.

In bulk water, ions interact via a Coulomb potential $\phi(r)=e/(4\pi\epsilon_{0}\epsilon_{w}r)$, where $\epsilon_{w}$ is the water dielectric permittivity. But inside a nanochannel, ions are no longer surrounded by a homogeneous fluid, and their interaction potential may be affected by the dielectric properties of the confining medium. In order to assess the importance of this effect, we consider the simple situation where an infinitely long cylindrical channel of radius $R$ is filled with water having isotropic permittivity $\epsilon_{w}$, and the medium outside the channel is a homogeneous dielectric of permittivity $\epsilon_{s}$, with typically $\epsilon_{s}\ll\epsilon_{w}$. The electrostatic potential around an ion placed in the middle of the channel can then be determined by solving Poisson’s equation in the presence of the dielectric discontinuities, which has been done analytically by several authors (Parsegian (1969), Levin (2006), Teber (2005), Kavokine et al. (2019)).

The complicated analytical result may be interpreted in the following way (see figure 8). At short distances (much smaller than the channel radius), only the dielectric response of the water is visible and the potential is $\phi(x)\sim 1/(\epsilon_{w}x)$. At long distances (much larger than the channel radius), it is the dielectric response of the confining medium that matters and $\phi(x)\sim 1/(\epsilon_{s}x)$. At intermediate distances, the electric field lines are essentially parallel to the channel due to the dielectric contrast $\epsilon_{w}\gg\epsilon_{s}$, and the potential, which is well described by an exponential function, resembles a 1D Coulomb potential:

with $\alpha$ a numerical coefficient that depends on the ratio $\epsilon_{w}/\epsilon_{s}$; $\alpha=6.3$ for $\epsilon_{w}/\epsilon_{s}=40$ (Teber (2005)). It is notable that so far no assumption on the channel radius was made, so that formally this 1D regime exists for a channel of any size. However, it is only relevant if it leads to ion-ion interactions stronger than $k_{B}T$, that is if $e\phi(R)>k_{B}T$, for a monovalent ion. This defines a limiting channel radius, below which ionic interactions are affected by the confining medium: $R_{c}\sim 7.5~{}\rm nm$. These modified Coulomb interactions therefore may have some effect in any single-digit nanopore, and they have essentially two practical consequences.

First, because ions produce a stronger Coulomb potential in confinement than in the bulk, they have a larger Born self-energy. This phenomenon was pointed out decades ago by Parsegian (1969) in the case of ions crossing lipid bilayers, and has since then been commonly considered in the study of biological ion channels. For small channels where the 1D contribution to the potential dominates, the additional self-energy due to confinement is ${\cal E}_{s}=e\phi(0)/2$, with $\phi(x)$ given by eq. (51). An expression valid for any $R$ was obtained by Teber, and it is plotted in figure 8. It shows that ${\cal E}_{s}$ is larger than $k_{B}T$ up to $R\sim 4~{}\rm nm$. ${\cal E}_{s}$ represents a supplementary energy barrier for entering the channel. To illustrate the consequence of this Born self-energy, let us consider a neutral channel in the absence of correlation effects: a barrier ${\cal E}_{s}=4k_{B}T$ reduces the ionic concentration by a factor $e^{{\cal E}_{s}/(k_{B}T)}=100$.

Second, increased Coulomb interactions with respect to the bulk may actually result in correlation effects. To assess the potential importance of correlations, one may introduce a coupling parameter $\Gamma=e\phi((\pi R^{2}\rho)^{-1})/(k_{B}T)$, where $\rho$ is the ion concentration inside the channel, so that $(\pi R^{2}\rho)^{-1}$ is the average distance between ions (Bocquet & Charlaix (2010)). Let us take as an example a channel of radius 2 nm, where ${\cal E}_{s}\sim 2.5k_{B}T$ from eq. (51). An effective Coulomb potential of magnitude $\sim 5k_{B}T$ then extends around an ion over a distance $\alpha R\sim 10~{}\rm nm$. Therefore, for concentrations above $\rho=10^{-2}~{}\rm M$, $\Gamma>1$ and correlations are expected to become important.

Several remarks can be made here. A first point is that we have considered monovalent ions so far; for ions of valence $z$, the self-energy, for example, would be multiplied by a factor $z^{2}$. Then, we have gone through the discussion assuming for simplicity that water has an isotropic permittivity $\epsilon_{w}$. We have highlighted in the previous section that this assumption breaks down in cylindrical channels of radius smaller than 5 nm, and, though it does not affect the qualitative phenomenology, dielectric anisotropy, as predicted by Loche et al. (2019), should then be taken into account. It mainly affects the behaviour of the potential at short distances $x\leq R$. We should also mention that when entering a channel smaller than its hydrated radius, an ion pays an additional energy penalty due to the shedding of its hydration shell (Richards et al. (2012), Epsztein et al. (2019)). For chloride, the total hydration energy is as high as $155~{}k_{B}T$, but the hydrated radius is 0.4 nm, so that only partial dehydration arises and only for the smallest pores. Finally, we have assumed that the confining medium is a dielectric medium with uniform permittivity. This is a good model for a lipid membrane, but artificial confining materials may have a variety of electronic properties (see 4.2.6), and their influence on effective Coulomb interactions has not been thoroughly investigated. In the case of confinement by a perfect metal, the ionic interactions were found to be exponentially screened over a distance of order of the channel radius (Kondrat & Kornyshev (2011), Loche et al. (2019)). Such effective interactions were found to lead to like-charge attraction in ionic liquids confined in carbon nanopores (Futamura et al. (2017)).

We have highlighted in the previous paragraph that correlations may become important in ion transport through single digit nanopores. If it is the case, then mean-field theories, such as the PNP framework, may not be directly applied, and the determination of equilibrium properties, such as ion concentrations inside the nanopore, requires some form of exact statistical mechanics. In systems with high aspect ratio (typically nanotubes), the problem can be reduced to a one-dimensional model of the Ising or 1D Coulomb gas type, and the partition function may then be exactly computed through a (functional) transfer matrix formalism (Zhang et al. (2006), Lee et al. (2014a), Kavokine et al. (2019)). Outside of the 1D geometry, variational methods may be used (Buyukdagli et al. (2010)), and methods for incorporating ion pairs into Poisson-Boltzmann theory have been developed (Levin (2002)). A general feature of these calculations (Zhang et al. (2006), Buyukdagli et al. (2010)) is that they predict filling transitions: namely, a non-analytic behaviour of ion concentrations in the nanopore as a function of the salt concentration in the reservoir, which strongly deviates from the mean-field Donnan equilibrium result, and may break local electroneutrality (Zhang et al. (2006)). One should note, however, that breakdown of electroneutrality does not require ionic correlations and it has recently been shown to arise in a mean-field setting (Levy et al. (2020)).

If linear response theory applies, then transport properties such as ionic conductance may be determined directly from the concentrations of charge carriers at equilibrium. However, this is no longer the case when correlations are very strong, leading to the formation of tightly bound Bjerrum pairs of oppositely charged ions. This was first realised by Onsager, who showed that in a three-dimensional weak electrolyte – an electrolyte that forms Bjerrum pairs – there is a quadratic current-voltage relation, a phenomenon known as the second Wien effect (Onsager (1934), Kaiser et al. (2013)). Bjerrum pairing in one-dimensional confinement was studied by Kavokine et al. (2019), and stable ion pairs were shown to arise typically for confinement below 2 nm, similarly to what was observed in MD simulations (Nicholson & Quirke (2003)). The pairing resulted in some very non-linear behaviour governed by discrete particle effects, echoing the phenomenology of Coulomb blockade (Beenakker (1991), Kaufman et al. (2015)). As such, ion correlations produced quantised transport behaviour in a purely classical system.

The coupling of ion transport to fluid transport in the non-mean-field, non linear regime remains largely unexplored, though it promises a rich phenomenology. As an example, a recent simulation (Li et al. (2017)) found that pressure driven water flow through a nanochannel could be blocked by an ion tightly bound inside. Lastly, we should mention, at the frontier of sub-continuum ion transport, the Grotthus-like translocation of protons, which has been observed in a single file (Tunuguntla et al. (2017)) and a single plane (Gopinadhan et al. (2019)) of water, but still remains poorly understood in strongly confined systems.