Ion transport measurements are widely used as an indirect probe for various properties of confined electrolytes. It is generally assumed that the ion concentration in a nanoscale channel is equal to the ion concentration in the macroscopic reservoirs it connects to, with deviations arising only in the presence of surface charges on the channel walls. Here, we show that this assumption may break down even in a neutral channel due to electrostatic correlations between the ions arising in the regime of interaction confinement, where Coulomb interactions are reinforced due to the presence of the channel walls. We focus on a one-dimensional channel geometry, where an exact evaluation of the electrolyte’s partition function is possible with a transfer operator approach. Our exact solution reveals that in nanometer-scale channels, the ion concentration is generally lower than in reservoirs and depends continuously on the bulk salt concentration, in contrast to the conventional mean-field theory that predicts an abrupt filling transition. We develop a modified mean-field theory taking into account the presence of ion pairs that agrees quantitatively with the exact solution and provides predictions for experimentally relevant observables, such as the ionic conductivity. Our results will guide the interpretation of nanoscale ion transport measurements.

## I. INTRODUCTION

A channel connects two reservoirs filled with a salt solution at concentration *c*_{out}. What is the salt concentration *c*_{in} inside the channel? The straightforward answer *c*_{in} = *c*_{out} is challenged as soon as the channel’s dimensions are at the nanometer scale.^{1} A deviation typically occurs because of the presence of a surface charge density Σ on the channel walls, which results in an imbalance of the concentrations $cin\xb1$ of the positive and negative ions. In a cylindrical channel of radius *R* that is smaller than the electrolyte’s Debye length, under the assumption of local electroneutrality, the concentrations are given by the famous Donnan equilibrium result,^{2}

Equation (1) is widely used to infer a channel’s surface charge from measurements of its conductivity at different salt concentrations. For sufficiently small surface charges (2Σ/*R* ≪ *c*_{out}), Eq. (1) predicts *c*_{in} = *c*_{out} even at extreme nanoscales. Importantly, this prediction underlies the method for extracting confined ion mobilities from transport measurements, which has been applied down to 7-Å-wide two-dimensional channels.^{3} Yet, physically, *c*_{in} = *c*_{out} stems from the assumption that the electrolyte solutions, both in the reservoirs and in the channel, behave as ideal gases of non-interacting ions. While such a description is valid in the bulk reservoirs at reasonable salt concentrations,^{4} it must be challenged in the nanometer-scale channel that is subject to *interaction confinement*^{5}—a reinforcement of the effective Coulomb interactions between the ions due to the dielectric contrast between the solvent (water) and the channel wall.^{2,5–13}

Due to interaction confinement, ions face a *self-energy barrier* *E*_{s} when entering the channel.^{6,7} It was first noted by Parsegian^{6} that this should result in ion exclusion: the salt concentration within the channel is then given by an Arrhenius scaling $cin=coute\u2212Es/kBT$ under the assumption of non-interacting ions. However, the result becomes more subtle as the confinement-reinforced ionic interactions are taken into account. Within a mean-field description of a spherical nanopore, Dresner^{14} predicted an abrupt filling transition, where *c*_{in} was a discontinuous function of *c*_{out}. Later, Palmeri and co-workers^{15,16} recovered a similar transition using a three-dimensional model of a cylindrical channel, treated within the variational field theory formalism of Netz and Orland.^{17} While this approach could be applied to a realistic geometry, it took electrostatic correlations into account only approximately.

An exact treatment of electrostatic correlations is possible upon simplification of the geometry to a purely one-dimensional model, with the channel wall being taken into account by introducing an effective confined Coulomb interaction. The 1D electrolyte can then be mapped onto an Ising or 1D Coulomb-gas-type model; the transfer matrix solution of such models was used, for example, to discuss the capacitance of nanoporous systems.^{18–20} The lattice models may be taken to the continuum limit, and the resulting path integral solutions have been used to discuss various ion-exchange phase transitions that arise in the presence of fixed discrete charges inside the channel^{8,21,22} and the ionic Coulomb blockade phenomenon.^{12} Such models are particularly rich theoretically, as they support a mapping to non-Hermitian quantum mechanics.^{23} Nevertheless, to our knowledge, the fundamental problem of ion filling in an uncharged channel has not been tackled within this framework.

In this paper, we treat the ion filling problem in the interaction confinement regime using an exactly solvable one-dimensional model. We find that the value of *c*_{in} is strongly affected by the formation of Bjerrum pairs—pairs of oppositely charged ions—within the channel, which precludes the occurrence of an abrupt filling transition. This is in contrast to the prediction of Palmeri and co-workers^{15,16} and to the result of the conventional mean-field theory. We then build on our exact results to propose a modified mean-field model that accounts for the relevant physical ingredients and, particularly, for the presence of ion pairs.

The paper is organized as follows: in Sec. II, we present the one-dimensional model and its solution within a path-integral formalism. The reader interested only in the physical outcomes may skip directly to Sec. III, where we discuss the model’s prediction for the ion concentration within the channel, compare it to the mean-field solution, and interpret it in terms of tightly bound Bjerrum pairs. In Sec. IV, we establish a modified mean-field theory, based on the notion of *phantom pairs*, which reproduces our exact solution. The mean-field theory allows us to determine the number of unpaired ions and produces experimentally relevant predictions for a nanochannel’s ionic conductance. Section V establishes our conclusions.

## II. 1D COULOMB GAS MODEL

### A. Confined interaction

We consider a cylindrical channel of radius *R* and length *L*, connected to macroscopic reservoirs [Fig. 1(a)]. We first assume, for simplicity, that the channel is filled with water that has an isotropic dielectric permittivity *ϵ*_{w} = 80 and that it is embedded in an insulating medium with a much lower permittivity *ϵ*_{m} (for a lipid membrane,^{6} *ϵ*_{m} ∼ 2). The effective Coulomb interaction *V*(*x*) between two monovalent ions separated by a distance *x* on the channel axis can be computed exactly by solving Poisson’s equation.^{7,11,12} A simple approximate expression can be obtained for *x* ∼ *R* (Ref. 2),

where *α* is a numerical coefficient that depends on the ratio *ϵ*_{w}/*ϵ*_{m} (*α* = 6.3 for *ϵ*_{w}/*ϵ*_{m} = 40). The reinforcement of electrostatic interactions compared to the usual *e*^{2}/4*πϵ*_{0}*ϵ*_{w}*r* Coulomb interaction that ions experience in bulk water can be interpreted in terms of image charges within the channel walls [Fig. 1(b)]. Two confined ions interact not only with each other but also with their respective image charges.

We introduce the parameters *ξ* ≡ *αR* and *x*_{T} ≡ 2*πϵ*_{0}*ϵ*_{w}*R*^{2}*k*_{B}*T*/*e*^{2}: both have the dimension of a length. With these notations,

The effects of ion valence and the anisotropic dielectric response of confined water can be taken into account by adjusting *ξ* and *x*_{T}.^{12} The dielectric response of water is expected to become anisotropic in channels with radius smaller than $\u223c5nm$:^{11,24} the component perpendicular to the channel axis is reduced, while the parallel component is largely unaffected. This results in ionic interactions that are even stronger than predicted by Eq. (2), described, in the framework of Eq. (3), by a larger *ξ* and smaller *x*_{T}.^{12}

Formally, the expression in Eq. (2) is valid for any channel radius. Yet, it is only physically relevant if at *x* ∼ *R*, the interaction is significant compared to *k*_{B}*T*, which restricts, in practice, the applicability of Eq. (2) to *R* ≲ 2 nm. In such extreme 1D confinement, we may neglect the ions’ degrees of freedom perpendicular to the channel axis and assume that they are constrained to move in one dimension. The partition function of such a 1D electrolyte may be computed exactly, as detailed in Sec. II B.

### B. Path integral formalism

Here, we detail the analytical solution for the partition function of a 1D Coulomb gas-like system that was first introduced in Ref. 12. We set *k*_{B}*T* = 1 until the end of Sec. II. We start from a lattice model in order to rigorously establish a path integral description in the continuum limit.

Our computation is inspired by the original solution of the 1D Coulomb gas model by Lenard and Edwards^{25} and subsequent studies by Démery, Dean, and co-workers,^{18,20,26,27} and Shklovskii and co-workers.^{21,22} We consider a one-dimensional lattice with sites 1, …, *M* as a model for the nanochannel of radius *R* and length *L*. Each lattice site *i* carries a spin *S*_{i}, which takes the values {0, 1, −1}, corresponding to no ion, a positive ion, or a negative ion occupying the site. We model the surface charge distribution as an extra fixed charge *q*_{i} added at each lattice site. The spins interact with the Hamiltonian

The system is in contact with a particle reservoir at concentration *c*_{out}. Here, the parameters *ξ* and *x*_{T} are dimensionless and expressed in the number of lattice sites.

The grand partition function is given by

with *z* = *c*_{out}*πR*^{2}*L*/*M* the fugacity. The matrix *C* can be analytically inverted,

Hence, we can carry out a Hubbard–Stratonovich transformation; that is, rewrite the partition function as a Gaussian integral, introducing the integration variable *φ*,

with $det(C)=e1/\xi 2sinh(1/\xi )\u22c5\xi (1\u2212e\u22122/\xi )M$. After performing the sum over the spins, which is now decoupled, we obtain

We now take a continuum limit of the lattice model. We call *a* the physical lattice spacing and let $\xi \u0303=a\xi $, $x\u0303T=axT$, and $z\u0303=Mz/L$. We then let *a* → 0 and *M* → ∞ while keeping the physical length of the system *L* = *aM* constant. We then drop the tilde sign to lighten the notation and obtain

with

*q*(*x*) is the one-dimensional density corresponding to the surface charge, and *z* ≡ *πR*^{2}*c*_{out}. At this point, *ξ* and *x*_{T} have the dimension of length. The path integral measure is defined as

We now define the propagator *P*(*φ*, *x*|*φ*_{0}, 0), or simply *P*(*φ*, *x*), as

Considering an infinitesimal displacement Δ*x*,

Expanding the propagator as *P*(*φ* − Δ*φ*, *x* − Δ*x*) = *P*(*φ*, *x*) − Δ*x*∂*P*/∂*x* − Δ*φ*∂*P*/∂*φ* + (1/2) (Δ*φ*^{2})∂^{2}*P*/∂*φ*^{2}, and carrying out the Gaussian integrals, we obtain

*P*(*φ*, *x*),thus, solves the partial differential equation

with initial condition *P*(*φ*, 0) = *δ*(*φ* − *φ*_{0}), which is the equivalent of a Schrödinger equation for the path integral representation (9). The partition function can, thus, be computed as

where *P*(*φ*, *L*|*f*_{0}) is the solution of (15) with initial condition $P(\phi ,0)=f0(\phi )\u2261e\u2212xT\phi 2/4\xi $.

### C. Transfer operator

We introduce the Fourier transform of *P* with respect to *φ*,

Then, $P\u0303(k,x)$satisfies

From now on, we restrict ourselves to an uncharged channel (*q* = 0). We then define the operator $T$ such that

which plays the role of a functional transfer matrix. Recalling Eq. (16), the partition function then reads

with $f0(k)=e\u2212\xi k2/xT$ and ⟨*f*(*k*)|*g*(*k*)⟩ ≡ ∫d*kf**(*k*)*g*(*k*).

Now, in the limit *L* → ∞, we may consider the largest eigenvalue *λ* of the operator $T$ and the associated eigenfunction *χ*,

Then, up to an exponentially small correction,

### D. Ion concentration

Our aim is to compute the salt concentration *c*_{in} in the nanoscale channel, given a salt concentration *c*_{out} in the reservoir. At the level of the lattice model, the probability to find, say, a positive ion at position *k* can be computed by replacing a factor (1 + 2*z* cos *φ*_{k}) by $zei\phi k$ in Eq. (8). In the continuum limit, we obtain the positive (negative) ion linear density at position *x* by inserting the operator *ze*^{iφ} (*ze*^{−iφ}) at position *x*,

Upon Fourier transformation, the insertion of *e*^{iφ} amounts to a shift by unity. Introducing the operator,

the concentrations are given by

since *z* = *c*_{out}*πR*^{2}. In the thermodynamic limit, and using Eq. (22) for the partition function, we obtain

## III. PHYSICS OF ION FILLING

### A. Debye–Hückel solution

We now go back to the ion filling problem [Fig. 1(a)] and present first a one-dimensional mean-field solution. Typically, the mean-field solution of an electrolyte problem is obtained by solving the Poisson–Boltzmann equation.^{28,29} For the conventional Poisson–Boltzmann equation to apply, we would need to consider the full three-dimensional geometry of our problem, and the effective interaction of Eq. (3) would be introduced implicitly through the boundary conditions at the channel walls.^{14} In order to obtain a mean-field solution directly in the 1D geometry, we need to introduce a modified Poisson’s equation for the electrostatic potential Φ whose Green’s function coincides with Eq. (3),

with *ϕ* ≡ *e*Φ/*k*_{B}*T* the dimensionless potential. Imposing that the ions follow a Boltzmann distribution (*c*_{±} = *c*_{in}*e*^{∓ϕ}, where *c*_{in} is understood as the average concentration inside the channel), we obtain the analog of the Poisson–Boltzmann equation in our 1D geometry,

In order to proceed analytically, we make a Debye–Hückel-type approximation and linearize Eq. (28) with respect to *ϕ*. Then, the potential around an ion placed in the channel at *x* = 0 is given by

with

The chemical potential inside the channel is the sum of an ideal gas entropic part and an excess part due to interactions,

with

Λ being the de Broglie thermal wavelength of the ions. *μ*_{ex} can be obtained via a Debye charging process,^{30}

We determine *c*_{in} by imposing equality of the chemical potentials between the channel and the reservoir,

which yields

Evaluating analytically the integral in Eq. (33), we obtain an implicit equation for *c*_{in}. With the notation $c\u0302in\u2261\pi R2cin$,

In Figs. 2(b) and 2(c), we plot the ratio *c*_{in}/*c*_{out} as a function of *c*_{out}, as obtained by numerically solving Eq. (36). We fix *ξ* = 7 nm (which corresponds to a channel with *R* ≈ 1 nm and strong dielectric contrast) and vary *x*_{T} to set the ionic interaction strength. The interaction strength may be quantified through the self-energy barrier, *E*_{s} = *k*_{B}*T* × *ξ*/(2*x*_{T}). The limiting behavior of *c*_{in}/*c*_{out} may be understood directly from Eq. (36). When *c*_{in} is small, Eq. (36) reduces to the Arrhenius scaling $cin=coute\u2212Es/kB\u2009T$: this result typically holds for biological ion channels that may contain either 0 or 1 ion at any given time, and the effect of inter-ionic interactions is negligible. When *c*_{in} is large, we recover *c*_{in} = *c*_{out}. Indeed, the excess term in the chemical potential vanishes at high concentrations, which is then dominated by the entropic term. The fact that *μ*_{ex} → 0 as *c*_{in} → ∞ is non-trivial: it can be seen, physically, as resulting from the Coulomb potential of each ion being perfectly screened by the other ions. At small values of *E*_{s}, Eq. (36) has a single solution for all values of *c*_{out}, which interpolates smoothly between the two limiting regimes. However, for *E*_{s} ≳ 5*k*_{B}*T*, it has three solutions in a certain range of *c*_{out}, pointing to a pseudo-first-order phase transition between a low-concentration and a high-concentration phase, similar to the one predicted by Dresner^{14} and Buyukdagli *et al.*^{15} The transition occurs at $c\u0302in\u223cxT/\xi 2$; as per Eq. (30), this corresponds to the concentration where the effect of the screening cloud on an ion’s Coulomb potential becomes significant.

### B. Full Poisson–Boltzmann solution

The physical content of the mean-field solution presented above is similar to the one of Dresner’s solution, based on a linearized Poisson–Boltzmann equation.^{14} The difference in geometry and the fact that he foregoes the use of the Debye charging process do not seem to play a significant qualitative role. The solution of Buyukdagli *et al.*^{15} takes ionic correlations into account to some extent, yet it still involves a Debye–Hückel-type linear equation for the mean-field interaction potential between the ions.

One may ask whether the same phenomenology persists if one does not linearize the Poisson–Boltzmann equation. The full Poisson–Boltzmann equation cannot be solved analytically but supports the following integral form:

where we have used the fact that *ϕ* should vanish at *x* → ∞. For *x* → 0, the solution of Eq. (37) should reduce to the unscreened potential in Eq. (3) up to an additive constant so that

Once again, one may express the excess chemical potential of the confined ions through a Debye charging process,

This result is the analog of Eq. (33), with *ϕ*_{λ}(0) now being the solution of an implicit non-linear equation, so that *μ*_{ex} must be determined numerically. As before, the concentration inside the channel is then given by

The prediction of the full Poisson–Boltzmann equation is shown in Figs. 2(b) and 2(c): we find *c*_{in} to be a smooth function of *c*_{out} for all values of parameters, in contrast to the linearized solution. We may not, however, unambiguously conclude that the filling transition is an artifact of linearization since the non-linear solution still involves a mean-field approximation and is not guaranteed to yield the correct result.

Interestingly, the “physically motivated” mean-field solution in Eq. (28) differs from the mean-field limit of our exact solution. It is obtained by taking the saddle-point approximation in the path-integral expression of the partition function [Eq. (9)]. The Euler–Lagrange equation for the minimizer *φ*(*x*) of the action *S*[*φ*] in Eq. (10) is, upon identifying *ϕ* = −*iφ*,

This is Eq. (28) with *c*_{in} replaced with *c*_{out} and corresponds to a first order treatment of interactions. Indeed, if the ions are non-interacting, *c*_{in} = *c*_{out}. By solving the mean-field equation, we determine how the ions’ chemical potential is affected by Debye screening, which then results in a value of *c*_{in} that is different from *c*_{out}. Within a straightforward interaction expansion procedure, one should determine the effect of screening assuming the zeroth order value for the ion concentration inside the channel, which is *c*_{out}: this corresponds to Eq. (41). Equation (28) contains an additional self-consistency condition as it assumes the actual value of *c*_{in} for the ion concentration, which is not known until Eq. (28) is solved. One may draw a loose condensed matter physics analogy, where Eq. (41) resembles the Born approximation for impurity scattering, while Eq. (28) is analogous to the self-consistent Born approximation.^{31}

### C. Exact solution

We now turn to the exact solution obtained in Sec. II to unambiguously solve the ion filling problem. We determine *c*_{in} according to Eq. (26) as

where *χ*(*k*) is the highest eigenfunction of the transfer operator in Eq. (19), determined, in practice, by numerical integration. The exact results for *c*_{in}, with the same parameter values as for the mean-field solution, are shown in Figs. 2(b) and 2(c). When interactions are weak [small values of *E*_{s}, Fig. 2(b)], the exact and mean-field solutions are in good agreement. Notably, all solutions smoothly interpolate between the bulk scaling *c*_{in} = *c*_{out} at a high concentration and the Arrhenius scaling $cin=coute\u2212Es/kBT$ at a low concentration. Conversely, in the strongly interacting case [large *E*_{s}, Fig. 2(c)], the exact result yields a much larger ion concentration than the mean-field solutions for intermediate values of *c*_{out}. In this intermediate regime, *c*_{in} remains a smooth function of *c*_{out} and obeys the scaling $cin\u221dcout2$.

Such a scaling is the signature of the formation of tightly bound Bjerrum pairs of positive and negative ions—strongly correlated configurations that are not taken into account by mean-field solutions. Indeed, let us assume that the channel contains an ideal gas of ion pairs at concentration *c*_{in}. We further assume that in a pair, the distance between the two ions is uniformly distributed in the interval [−*x*_{T}/2, *x*_{T}/2] and the binding energy of a pair is *k*_{B}*Tξ*/*x*_{T} = 2*E*_{s}. Then, the grand partition function reads

where we recall that *z* = *πR*^{2}*c*_{out} and *β* ≡ 1/(*k*_{B}*T*). Using that

we obtain

We recover, indeed, the quadratic scaling.

We may check that the prefactor in Eq. (46) is the correct one by evaluating analytically the expression in Eq. (26) in the low concentration limit *z*_{T} ≡ *zx*_{T} ≪ 1. An analytical expansion of the function *χ*(*k*) in powers of *z*_{T} was derived in Ref. 12. Substituting it into Eq. (26), we obtain

The first term in the expansion corresponds to $cin=coute\u2212\beta Es$. At the lowest salt concentrations, forming Bjerrum pairs is too entropically unfavorable and the concentration inside the channel is controlled by the self-energy barrier. However, as the salt concentration increases, there is no abrupt transition to a highly screened concentrated phase inside the channel; instead, the channel is progressively filled by Bjerrum pairs. This corresponds to the quadratic term in the expansion, with the prefactor agreeing, indeed, with Eq. (46). (This justifies *a posteriori* our choice of [−*x*_{T}/2, *x*_{T}/2] as the interval in which a paired-up ion is allowed to move.) The expansion in Eq. (47) reproduces quite well the low-concentration behavior of the exact solution as shown in Figs. 2(b) and 2(c). However, it fails at high concentrations, where it does not recover *c*_{in} = *c*_{out}.

Our exact analysis of the ion statistics in a nanoscale channel has revealed that Bjerrum pairs are a crucial ingredient of the filling process. We now develop a modified mean-field theory that accounts for the presence of Bjerrum pairs and compares it to the exact solution.

## IV. PAIR-ENHANCED MEAN-FIELD THEORY

### A. Debye–Hückel–Bjerrum theory

The traditional mean-field treatment of electrolytes is incapable of taking Bjerrum pairs into account, as it naturally neglects any strong ion–ion correlations—pairing being a fundamentally discrete phenomenon. An idea proposed by Bjerrum to amend the Debye–Hückel theory was to introduce ion pairs as a separate species encapsulating all “strong” ion–ion correlations.^{32,33} More precisely, any two oppositely charged ions that are closer than some minimum distance can be considered a single neutral entity—a Bjerrum pair. The remaining “free” ions should then only experience weak interactions with each other and can be treated at the mean-field level. Importantly, this last remark justifies the Debye–Hückel linearization, as all non-linear effects are assumed to be hidden in the definition of ion pairs.

As before, we consider that pairs behave like particles of an ideal gas and that their maximum extension is given by *x*_{T}. Defining $cinp$ the concentration pairs inside the channel, the chemical potential of pairs is given by

where the geometrical factor inside the logarithm accounts for the internal degrees of freedom of a pair. The chemical potential only has an entropic term because the binding energy of the pair exactly compensates for the self-energy of the two separate ions. The chemical equilibrium between free ions and pairs inside the channel can be written as

where $\mu in+$ and $\mu in\u2212$ are the chemical potentials of cations and anions, respectively. We then obtain, using the Debye–Hückel solution for *μ*_{in} [Eqs. (31)–(33)],

which is the result obtained in Sec. III C. The average concentration of free ions $cinf$ is not modified compared to the Debye–Hückel solution and is, therefore, the solution of the self-consistent Eq. (36). One can then compute the total concentration inside the channel as $cin=cinf+cinp$ or, explicitly,

In other words, the only impact of pairs in Bjerrum’s computation is to add a quadratic term $2\pi R2xTcout2$ to the Debye–Hückel result, matching with the expansion (47) of the exact solution up to order 2 in the bulk concentration. We compare the two predictions in Fig. 3(b). The Debye–Hückel–Bjerrum solution is found to match the exact one quite well at low and intermediate concentrations. This result is, however, unphysical for *c*_{out} ≳ 1/*πR*^{2}*x*_{T}: *c*_{in} is found to grow much faster than the bulk concentration. One solution would be to consider higher-order terms in the mean-field treatment through the inclusion of triplets, quadruplets, etc., of ions and all possible interactions between these entities. Truncating the sum at any finite order, however, would not yield a valid solution in the entire range of concentrations, nor is it guaranteed to converge to the exact solution. This approach is also unsatisfactory as it would not yield a closed-form expression for *c*_{in} and would not allow for a qualitative understanding of the underlying physics.

Instead, we develop a different method that, through physics-driven arguments, prevents the divergence of *c*_{in} at high bulk concentrations and quantitatively reproduces the exact solution.

### B. Phantom pairs

Equation (51) overestimates the number of Bjerrum pairs in the channel because it fails to account for the presence of Bjerrum pairs in the reservoir. The electrolyte in the reservoir is treated as an ideal gas: the ions are non-interacting and they cannot form actual, tightly bound pairs. Nevertheless, we have defined any two oppositely charged ions that find themselves in a cylinder of radius *R* and length *x*_{T} to be a separate chemical species. Such configurations may arise in the reservoir simply out of statistical chance: we dub them *phantom pairs*. For our mean-field theory to be consistent, these phantom pairs need to be taken into account.

We incorporate phantom pairs in our model using a thought experiment: let us assume that the nanochannel does not directly connect to the reservoir but is, instead, in contact with another 1D channel of similar geometry, where interactions are turned off. This intermediate 1D reservoir is, in turn, connected to the real, 3D reservoir. Let $coutp$ (resp. $coutf$) be the concentration of the aforementioned phantom pairs (resp. free ions) in the 1D reservoir. The chemical equilibrium between phantom pairs and free ions imposes

In addition, one has $coutf+coutp=cout$: since the interactions are negligible, the total concentration of ions in the intermediate reservoir should be equal to the salt concentration in the real reservoir. Imposing this condition yields

We use this result to control the proliferation of pairs in the channel: we now equilibrate the free ions inside the nanochannel with only the free ions in the 1D reservoir,

which corresponds to Eq. (35) with *c*_{out} replaced by $coutf$. Equation (54) is again a self-consistent equation; this time on the concentration of free ions $cinf$ that must be solved numerically. Finally, equilibrating pairs with free ions inside the channel (or, equivalently, pairs inside with pairs outside), we obtain

where the second term corresponds again to Bjerrum pairs. Equations (53)–(55) constitute the main result of our modified mean-field theory. Note that *μ*_{ex} may be determined at the Debye–Hückel level [Eq. (33)] or by solving the full Poisson–Boltzmann equation [Eq. (39)]. In what follows, we will only discuss the latter, as it offers greater accuracy; however, the Debye–Hückel prediction provides reasonable results even in the case of strong interactions and yields a convenient analytical expression for *μ*_{ex} as a function of $cinf$.

The prediction of our phantom pair Poisson–Boltzmann model is compared to the exact solution (26) in Fig. 3(b). The two solutions are found to be in near perfect agreement for all values of the parameters, even in the strong coupling limit *E*_{s} ≫ *k*_{B}*T*. This result shows the importance of taking into account phantom pairs in the reservoir to solve the ion filling problem.

### C. Conductance

A strength of our modified mean-field model is that it offers insights into the physical properties of the confined system beyond the value of the ionic concentration. In particular, the decomposition of the electrolyte into free ions and bound pairs allows us to estimate the channel’s conductance. Tightly bound Bjerrum pairs are electrically neutral so that they do not contribute to the ionic current to the first order in an applied electric field: it would then be straightforward to assume that the channel’s conductance is proportional to the concentration of free ions. However, the reasoning needs to be more subtle, since the channel, in the same way as the reservoir, may contain non-interacting phantom pairs. Indeed, we have decomposed the confined electrolyte into tightly bound pairs that have no ionic atmosphere and free ions that are dressed by a Debye screening cloud. As the concentration increases, the interaction between dressed ions becomes weak and two of them may find themselves within a distance *x*_{T} without actually binding. Such a phantom pair is expected to still contribute to the conductance. The concentration of phantom pairs in the channel is obtained by imposing their chemical equilibrium with the free ions treated as an ideal gas. Thus, we estimate the channel’s conductance as

where *D* is the diffusion coefficient of ions and the second term corresponds to the contribution of phantom pairs. In Fig. 4(a), we compare this result to the Ohm’s law prediction where pairs are neglected and one assumes *c*_{in} = *c*_{out}. Ohm’s law is found to greatly overestimate the conductance at low concentrations. In the dilute limit, we, instead, recover the Arrhenius scaling, where one assumes $cin=coute\u2212Es/kBT$.

Finally, we stress that Eq. (56) only accounts for the electrophoresis of free ions and is, therefore, only valid in the limit of weak external electric fields. Stronger voltage drops will result in the breaking of ion pairs, causing a conductivity increase in a process known as the second Wien effect. This phenomenon is described in Refs. 12 and 13 and has been used to create solid-state voltage-gated nanochannels.^{34} Furthermore, stronger electric fields may result in an unconventional transport mechanism involving the breakdown of the water-ion chain inside the channel, in apparent violation of the Nernst–Einstein relation.^{35}

### D. Effect of a surface charge

Until now, we have restricted ourselves to channels with uncharged walls. However, in most experimentally relevant situations, the channel walls bear a surface charge density Σ, which strongly impacts nanofluidic transport. While introducing a surface charge is tedious within the exact framework, we may readily assess the effect of a surface charge in the interaction confinement regime using our pair-enhanced mean-field theory.

In the limit where the channel’s radius is smaller than the Debye length, we assume that the presence of the surface charge amounts to a homogeneous Donnan potential drop *V*_{D} inside the channel, which we do not need to determine explicitly. Then, the chemical potential of ions inside the channel reads

Note that the concentration in free anions $cin\u2212$ and cations $cin+$ are now distinct so that *μ*_{ex} is defined as a function of the average free ion concentration $cinf=(cin++cin\u2212)/2$. In a channel of finite length, the ion concentrations do not necessarily obey local electroneutrality.^{36} It may, indeed, be more favorable for counterions to form a cloud at the channel mouths instead of entering the channel. At low salt concentrations, this can occur even in very long isolated channels, but the effect is mitigated in a membrane configuration with many parallel channels.^{37} Such non-neutrality effects are, however, beyond the scope of this work, and we assume in the following that the channel is sufficiently long for local electroneutrality to hold,

Imposing chemical equilibrium with the reservoir, we obtain a modified version of the Donnan result [Eq. (1)]:

with $coutf$ given by Eq. (53).

One can again obtain the channel’s conductance through Eq. (56), which we compare to the Donnan/Ohm’s law result in Fig. 4(b). Importantly, the Donnan result predicts that conductance becomes independent of concentration for *c*_{out} ∼ 2Σ/*R* [see Eq. (1)]. In practice, this result is commonly used to estimate experimentally the surface charge as Σ ∼ *Rc**/2, where *c** is the reservoir concentration for which the conductance levels off. In contrast, in the interaction confinement regime, we predict that the transition occurs, instead, at $cinf\u223c2\Sigma /R$—corresponding to a higher reservoir concentration, due to the self-energy barrier. In this case, Donnan’s prediction overestimates the surface charge by typically one order of magnitude, as shown in Fig. 4(b).

Finally, let us stress that we considered here a charge homogeneously distributed along the channel’s surface. This assumption is relevant in the case of conducting wall materials, such as systems where the charge is imposed via a gating electrode connected to the channel walls. This situation may, however, be different in experimentally available devices, where the surface charge generally consists in localized charged groups and defects on the channel walls. In this case, the physics become more involved as ions may form bound pairs with the fixed surface charges. Some of these physics have been revealed by the exact computations of Shklovskii and co-workers;^{8,21} a technically simpler approach to these physics using our pair-enhanced mean-field theory would be possible but extends beyond the scope of the present work.

## V. DISCUSSION AND PERSPECTIVES

We have determined the salt concentration inside a nanometric channel connected to reservoirs filled with electrolytes. In the case of a fully 1D geometry, corresponding to a nanotube of radius *R* ∼ 1 nm, we developed an exact field-theoretical solution that allowed us to compute channel concentration *c*_{in} as a function of the reservoir concentration *c*_{out}. This solution clears up the ambiguities of pre-existing mean-field theories and contradicts the naive expectation that *c*_{in} = *c*_{out}. In particular, the concentration inside the nanochannel is found to be always lower than in the bulk, as the confinement of electrostatic interactions creates an energy barrier for ions to enter the channel.

Yet, we found that *c*_{in} is, in fact, higher than the prediction of the mean-field Debye–Hückel theory, as ion pairing counterbalances, to some extent, the energy cost of interaction confinement. Such strong ion–ion correlations cannot be directly accounted for in a mean-field theory, and the filling transition that emerges in the Debye–Hückel theory appears to be an artifact of linearization. To overcome this issue, one can add Bjerrum pairs as a separate chemical species within the Debye–Hückel model. Carefully accounting for the statistical formation of unbound *phantom pairs*, we obtain a modified mean-field theory that reproduces the result of the exact computation with a near-perfect accuracy and that can be extended to account for a non-zero surface charge on the channel wall.

Despite the concurring results, the two original formalisms developed in this work serve distinct purposes. The field-theoretical solution plays the role of a touchstone model owing to its exact treatment of all many-body interactions. Modeling electrolytes is a notoriously hard problem in statistical physics, and simplified models often lack a reference solution for benchmarking their approximations. This difficulty is lifted in the 1D geometry: thanks to the existence of the exact solution, we have been able to build a quantitatively precise mean-field model, adding step-by-step the qualitative ingredients necessary to reproduce the exact result.

Moreover, the field theory formalism gives access to the entire statistics of the system, including finite-size effects that elude any mean-field treatment. The latter are expected to be relevant in many experimental situations, as a substantial amount of current work focuses on short pores, where the length of the channel is comparable to its radius. For instance, one can expect shorter channels to deviate from electroneutrality^{36}—something entirely impossible in the limit of infinitely long channels.

On the other hand, our modified mean-field formalism has the advantage of mathematical simplicity, allowing for convenient physical interpretations. The simple distinction between free ions and Bjerrum pairs can be used to straightforwardly estimate the channel’s conductance. The influence of ion–ion correlations on the conductivity is of particular importance as conductance measurements underpin many nanofluidic experiments. In contrast, the exact solution does not provide any such insight into transport properties, as it is limited to thermal equilibrium.

Furthermore, the mean-field model may easily be adapted to other geometries, whereas an exact treatment is only possible in the strictly 1D case. Extensions of our results to 2D nanochannels would be of significant interest. In particular, 2D nanochannels can be made out of various materials with different electronic properties, which directly impact the confined ionic interactions.^{5} Therefore, 2D nanochannels could serve as a platform for exploring the impact of wall metallicity on the ion filling problem.

Both our exact and mean-field solutions can be expected to fail at very high concentrations. Indeed, our work relies on a simplified picture of electrolytes, where all steric effects are discarded. We considered point-like ions with no short-distance repulsion; therefore, no effects such as saturation or layering can be accounted for. Similarly, we neglected any interaction with the solvent—for example, we did not consider the decrement in relative permittivity at high salt concentrations.^{38} However, since all electrostatic interactions are screened in the limit of high concentrations, such considerations should not impact the conclusions of the present work: particularly, we would still expect that *c*_{in} = *c*_{out} at high concentrations.

Finally, let us briefly recall our results for the ion filling problem. In channels larger than a few nanometers, the conventional mean-field picture is valid so that in the absence of any surface charge, the salt concentration inside the channel equals that of the reservoirs: *c*_{in} = *c*_{out}. For nanometer-scale confinement and low concentrations, interaction confinement amounts to a finite energy barrier for ions to enter the channel: $cin=coute\u2212Es/kBT$. As concentration increases, more ions are able to overcome the barrier by forming Bjerrum pairs, neutralizing the electrostatic cost of confinement, at the price of entropy: $cin\u221dcout2$. Only at high concentrations can one recover the intuitive estimate *c*_{in} = *c*_{out}, as intense screening cancels out all electrostatic interactions. Overall, the interaction confinement has a significant impact on the properties of nanofluidic systems and the assumption that *c*_{in} = *c*_{out} should be questioned any time the system’s size reaches the nanometer scale.

## ACKNOWLEDGMENTS

N.K. acknowledges the support from a Humboldt fellowship. L.B. acknowledges the funding from the EU H2020 Framework Programme/ERC Advanced Grant Agreement No. 785911-Shadoks. The Flatiron Institute is a division of the Simons Foundation.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Paul Robin**: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Adrien Delahais**: Investigation (supporting). **Lydéric Bocquet**: Conceptualization (equal); Writing – review & editing (equal). **Nikita Kavokine**: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.