Confinement-controlled rectification in a geometric nanofluidic diode

Recent experiments with electrolytes driven through conical nanopores give evidence of strong rectified current response. In such devices, the asymmetry in the confinement is responsible of the non-Ohmic response, suggesting that the interplay of entropic and enthalpic forces plays a major role. Here we propose a theoretical model to shed light on the physical mechanism underlying ionic current rectification (ICR). By use of an effective description of the ionic dynamics we explore the system's response in different electrostatic regimes. We show that the rectification efficiency, as well as the channel selectivity, is driven by the surface-to-bulk conductivity ratio Dukhin length rather than the electrical double layer overlap.


I. INTRODUCTION
Dating back to the famous thought experiment of Maxwell's demon (1867), the dream of designing force-free transport devices has permeated different branches of physics, including nanofluidics. In the context of nanofluidics, one can imagine the ionic diode 1 , a nanofluidic device exhibiting ionic currents of unequal magnitude under voltages of equal magnitude and opposite polarity, as a realization of such a demon. The first realization of such a nanometric ionic diode was reported by Siwy and Fuliński in a geometrically asymmetric nanochannel obtained by asymmetric chemical etching of a polymer foil 2 . Their conical channel demonstrated a strongly non linear ionic current under ac voltage, resulting in a net average current under zero average forcing. Ionic current rectification (ICR) in conical nanochannels has since been extensively studied experimentally [3][4][5][6][7] , thanks to the considerable progress made over the last twenty years in nano-fabrication technologies 8 . ICR has also been observed in symmetric channels subject to a concentration gradient 9 and in the presence of a surface charge discontinuity 10 . Empirically, the two features necessary to observe current rectification have been identified as the presence of surface charge and broken symmetry in the direction of transport, irrespectively of the nature of the broken symmetry. Alongside practical applications in macromolecular sensing and manipulation 11,12 , energy harvesting [13][14][15] and water desalination 16,17 , the phenomenon raises fundamental questions on the nature of ionic transport at the nanoscale. At this lengthscale, surfaces and entropic confinement strongly influence mass transport leading to the emergence of non-linear and exotic responses 18 , of which ICR is a prominent example. A rationalization of the latter would then be a testbed for understanding more complex behaviour occurring at the nanoscale 18 such as that of biological functionalized protein channels [19][20][21] .
In nano-sized fluidic diodes, electrostatic interactions between charged species play a key role. In the presence of a surface charge density σ in contact with an electrolyte solution, an electrical double layer (EDL) builds up inside the channel with a characteristic decay length given by the Debye length, over which the imbalance of charge due to the channel walls is screened, Here, k B is the Boltzmann constant, T is the temperature, 0 and w are respectively vacuum and relative water permittivities, c s is the bulk electrolyte concentration, z is the electrolyte valency and e is the elementary charge. At room temperature w ≈ 80 and the Debye length can span from tens of nanometers down to a feẘ Angstroms depending on the salt concentration. Within the EDL an excess counterion concentration screens the surface charge giving rise to an electrically charged region. In the so-called entropic electrokinetic regime 22 the Debye length is comparable to the tip of the nanopipette, i.e., the smallest aperture. This is typically the case in most synthetic realizations of nanochannels 6,23,24 as well as in biological ion channels. Notably, measures of ICR in micrometer-sized systems have been reported more recently in the literature 25,26 . It is convenient to introduce a second electrostatic length known as the Dukhin length: which quantifies the relative importance of surface compared to bulk transport. Contrary to λ D , the Dukhin length is a phenomenological length: it does not directly correspond to a physically observable length in the system. Therefore it can be much larger or smaller than the system's size 18 . Eq. (2) indicates that the Dukhin length can be understood as the ratio between two different lengthscales, namely the Debye length (1) and the Gouy-Chapman length l GC = 2 0 w k B T /z 2 e|σ|, defined as the typical length at which the surface electrostatic potential energy equals the thermal energy. Nonetheless, the electrostatic phase space, associated with the surface charge σ and the bulk concentration c s degrees of freedom, is determined only by two independent lengths. As will become clear in the following, the choice of λ D and l Du as model parameters is convenient for the problem at hand. In analogy to colloidal science, we can also introduce a dimensionless Dukhin number 27 : whereh is the average half height of the channel. Du 1 identifies the regime globally dominated by surface transport. The theoretical literature on ICR has been confined mostly to numerical simulations of the ion dynamics using the classical Poisson-Nernst-Planck (PNP) equations for dilute electrolyte solutions [28][29][30][31] . Such a framework has quantitatively captured the phenomenon, demonstrating that a mean field continuum description is still valid for ionic dynamics down to a few nanometers. An early qualitative interpretation of ICR is traceable back to a paper of Dietrich Woermann 32 who rationalized the phenomenon in terms of ionic transference asymmetry between the ends of the channel. At the same time, the study of particle transport over entropic barriers has attracted the attention in non-equilibrium statistical physics 22,[33][34][35] . The first attempt to characterize transport in confined systems dates back to the early work of Jacobs 36 and Zwanzig 37 who proposed the so-called Fick-Jacobs approach (FJ) to account for the transport of Brownian particles geometrically confined in a quasi-one-dimensional system. Under the assumption of a separation of scales between the longitudinal and the transversal coordinates, the latter is integrated and the description is reduced to an effective 1D equation now containing an entropic term. The validity of the approach has been tested both in the case of free diffusion 38,39 and in the presence of an external force 40 and demonstrated to be quantitatively accurate for channel geometries with smoothly varying cross-section under moderate external field, typically requirements that are satisfied in nanofluidic setups 6,8,41 . Overall the FJ approach represents a well-established systematic framework to describe transport in the presence of entropic barriers, and it has been recently extended to the regime of competition between energetic and entropic interactions in electrolyte dynamics 42 . Our goal in the present work is to gain insights on the fundamental mechanism controlling current rectification in a geometric diode, i.e., a conical channel with uniform charge density in contact with two reservoirs held at the same electrolyte concentration. Such a configuration corresponds to an extensively studied nanopipette experimental setup. Moreover, it represents the conceptually FIG. 1. Schematic view of the channel in contact with two reservoirs at fixed salt concentration. The channel width is assumed to be constant along the z direction pointing out of the page. The channel walls carry a uniform negative charge density and a electrically charged double layer forms over a characteristic length λD.
intriguing case in which symmetry breaking originates only from the geometric confinement; such a system is thus able to harness entropy to rectify ionic current.
To address the problem we adapt the FJ approach to a 2D conical slab geometry. Contrary to previous works considering channels much larger than the Debye length 7,43 , the present formalism allows us to investigate the regime of finite λ D where partial Debye overlap occurs inside the channel, and to fully capture the interplay between energetic and entropic contributions. Furthermore, we are able to derive analytical predictions for the limiting conductance in the regime of strong EDL overlap which, to the best of our knowledge, have not yet been derived for the geometric diode. Finally, our results assess the key role played by the Dukhin number in the microscopic mechanism of rectification providing further insight on the nature of ICR.

II. IONIC DYNAMICS
As shown schematically in Fig. (1), we consider an open asymmetric channel with a slab geometry characterized by longitudinal size L, width L z and an x-dependent height whereh is the half-aperture of the channel and k = |d x h| = (h L − h R )/L is the difference between the left h L and the right h R channel half-heights in units of channel length. In the following sections the channel slope is varied by keeping fixed its half-heighth in order to compare systems with the same aspect ratio. The channel is filled with a symmetric monovalent electrolyte composed of species having equal diffusion coefficient D, in contact with two reservoirs at fixed temperature T and ionic strength c s . Each wall bears a uniform negative surface charge of density σ < 0. We assume L z h so that we can neglect the z dependence of any variables of the model and the resulting system is effectively 2D. In order to characterize the ionic dynamics we derive effective one-dimensional transport equations for the ionic concentration profiles c ± . The approach relies on the constraint of a small aspect ratio =h/L 1 , i.e. a slowlyvarying channel geometry. In this case, the transversal relaxation dynamics with characteristic relaxation time τ y ∼h 2 /D is decoupled from the longitudinal relaxation dynamics with τ x ∼ L 2 /D and the ions are assumed to instantaneously adjust to the Boltzmann distribution at each cross-section. Such separation of scale is known in the literature as local thermodynamic equilibrium 44 (LTE). Under these assumptions the steady-state Nernst-Planck equation for the positive and negative ionic species reads: where j ± is the constant mass flux density along x, β = 1/k B T and Φ(x, y) is the total electrostatic potential inside the channel. We have neglected in (5) the advective flux which proved to be minor compared to the electrophoretic contribution for moderate surface charge densities and moderate external fields 29 . Eq. (5) must be supplemented by the Poisson equation relating the electrostatic potential to the spatial charge distribution q = e(c + − c − ) inside the channel: In the next section we reduce (5) to an effective 1D equation by introducing the FJ ansatz for the ionic concentration profiles as explained in II A. For consistency, the same approximation is applied to the Poisson equation together with the assumption of small transversal variation of Φ (see section II B), which allows to formally integrate (6).

A. The Fick-Jacobs approach
Since the ionic transversal and longitudinal dynamics are assumed to be decoupled, it is convenient to introduce the marginal concentration as the cross-sectional integral of the volumetric concentration: Moreover, following the approach of Zwanzig 37 we define x−dependent free energies A ± (x) via: From the hypothesis of LTE we may factor the volumetric concentrations c ± (x, y) into the product of equilibrium normalized conditional densities ξ ± (y; x) and the marginal concentrations, dy e ∓βeΦ(x,y) ·c ± (x).
(9) Eq. (9) represents the key ansatz of the FJ approach. Martens et al 33 proved that (9) can be recovered as the zero-order term of a perturbative expansion in series for the geometrical parameter k around the zero-transversal-flux solution. Notably for the case of a conical channel, where |d x h(x)| = const, taking into account the extra x−dependence of the diffusivity D(x) amounts to a rescaling of the diffusion coefficient thus making the theory here developed valid up to k ≤ 1 39,45 .
In the present work we examine the zero-order FJ approximation and we leave to future work the discussion of higher order corrections.
By integrating Eq. (5) in the y coordinate and using (9) as a closure for c ± (x, y), an effective one-dimensional equation is obtained, (10) where J ± = dy j ± is the longitudinal mass flux per unit width for each species. In Eq. (10) the concentrations c ± (x) are the marginal ones; in the following, we refer to the marginal concentrations unless the both x− and y−dependences are explicitly noted. Now we introduce dimensionless variables. As reported in table I we rescale the x coordinate by the total length of the channel L and the coordinate y as well as the Debye length λ D and the channel profile h(x) by the half-height h. In this way the channel profile reads where we also introduced a rescaled channel slope κ = k/ . Since we keep fixed the half heighth allowing for variation in the degree of corrugation we note that κ < 2 for geometrical consistency. The electrostatic potential is rescaled by the thermal one k B T /e and the volumetric concentrations c ± by the concentration in the bulk c s . Consequently, the charge density q is rescaled by ec s , the mass flux J per unit width by Dc sh /L and the conductance per unit width, G = ∂I/∂∆V with I the total ionic current and ∆V the applied potential drop, by the bulk conductance Dc sh e 2 /k B T L.
In dimensionless form Eq. (10) now reads: where we have introduced the (dimensionless) electrochemical potential µ ± = log c ± + βA ± . In Eq. (12) the electrophoretic contribution now appears in terms of the previously introduced effective free energies A ± (x) . For a neutral species the effective free energy reduces to the standard Boltzmann entropy βA(x) = − log 2h(x). In this case d dx βA(x) is referred to as an entropic force, originating from the variation in phase-space volume available for free diffusion along the channel. For a charged species A embeds both enthalpic and entropic contributions. Eq. (12) must be integrated with the appropriate boundary conditions, i.e. by imposing continuity in the electrochemical potential at the ends of the channel 46 . The discontinuity in the surface charge distribution at the channel ends, and the consequent readjustment of ions within the diffusive layer results in an apparent local discontinuity in the concentration and electrostatic potential profiles. In the present framework, such discontinuities are treated as point-like discontinuities, which stands for the fact that these entrance effects are O( ), hence they fall within the level of our approximation.

B. Local Debye-Hückel approximation
In dimensionless units, the Poisson equation (6) reads It is convenient to decompose the electrostatic potential as where φ = 1 is the average potential across y, ψ = φ − φ in the excess potential at each section and φ ext = −∆V x − 1 2 is the potential drop applied externally, resulting in a constant electric field directed in the x direction. By using FJ approximation into Eq. (13) together with (14) and by linearizing in ψ under the assumption of small potential variation in the transversal direction, we reduce (13) to We refer to the linearization used to derive Eq. (15) as a local Debye-Huckel (DH) approximation: the potential is linearized with respect to the local cross-sectional average preserving therefore global nonlinearity. We stress that the assumption of small ψ is more general than standard DH, which requires small ζ potential everywhere (typically 47 ζ ≤ 25mV ). In fact it allows to explore the ideal gas 48 regime for arbitrarily high Dukhin number, where typically the global Debye-Hückel assumption would fail. The small aspect ratio constraint allows for a lubricationlike approximation of (15) which reduces to a linear equation for ψ: Consistently, the scaling argument applies as well to the electrostatic wall boundary condition, which after neglecting terms of O(k 2 ), and introducing rescaled variables, reduces to We shall note here that the LTE hypothesis previously introduced implies local electroneutrality, in which the integrated charge density balances the surface charge density at each cross-section, In fact, by integrating Eq. (15) in the y−coordinate, using (17) and neglecting O( 2 ) terms, Eq. (18) is obtained.
To first order approximation the FJ ansatz, the lubrication approximation and local electroneutrality are different naming for the same unique assumption, i.e. separation of transversal and longitudinal scales, applied to different physical properties 49 . The reduced Poisson equation can be formally integrated leading to In Eq. (19) the potential is naturally expressed in terms of a local Dukhin number and a local Debye length respectively defined as: Du(x) = Du c vol (x) (21) in terms of the total average volumetric concentration We recognize the first term on the rhs of Eq. (19) to be the Debye-Huckel potential carrying an extra x−dependence due to the varying channel geometry. The second term on the rhs ensures local electroneutrality. Eqs. (12) and (16) need to be solved numerically. It is convenient to rewrite Eq. (12) in terms of ψ : so that the coupling between the concentration profiles and the electrostatic potential is made now explicit. We use finite-element simulations (COMSOL) to solve the system of Eqs. (22a)(22b)(22c) in order to look at the electric current I = J + − J − generated by the applied potential drop ∆V . (See appendix for details on the numerical simulations).
The expression for the electric current obtained by formally integrating Eq. (22a) and (22b) reads where the denominator is responsible for the non-linear (rectified) response of the channel, as it expresses the coupling between the dissipative dynamics (thermodynamic forcing) and the geometric asymmetry. For a flat channel, Eq. (23) reduces to the standard ohmic response 50 (per unit width) which in dimensional unit reads Eq. (23) is valid for slowly varying channels under the assumption of small potential variation in the transversal direction. Hence it represents a well-grounded expression for the ionic current allowing to span across different regimes in the electrostatic phase space in both λ D and Du.
Previously proposed analytical approaches 32,43,51 assume that λ D is the relevant controlling parameter by treating separately the case of no overlap λ D 1 and strong overlap λ D 1. This is not necessary in the present framework, where λ D can vary continuously. Nevertheless, it is useful at this stage to introduce the regime of strong Debye overlap as it represents a well-known scenario which we will use as a benchmark to compare with numerical results. Let us consider the regime in which the channel height is much smaller than the Debye length. The EDL extends all throughout the interior of the confined electrolyte, rendering the channel perfectly charge-selective. Both the electrostatic potential, Φ(x), and the ionic concentration profiles, c ± (x), are assumed to be uniform in the transversal direction allowing for a substantial simplification of the mathematical problem at hand. We stress that the concentrations c ± (x) here are not the marginal concentrations but the total concentrations which in this limit are independent of y. Together with local electroneutrality which in this case reads continuity in the chemical potential provides an expression for the Donnan potential at either end of the channel 17,28,50 , Notably, already at equilibrium the varying geometry results in a non-uniform tilted potential across the channel. Analogously, the channel's junctions concentrations read Hence, a jump in concentration profiles builds up at each junction of the channel to compensate for the potential discontinuity (26). Such a local balance is known in the literature as local Donnan equilibrium. These expressions will allow for asymptotic analytical predictions for the conductances when ∆V → ±∞. The equation of motion (12) for λ D 1 reduces to which, rewritten in terms of the total mass flux J = J + + J − and electric current I, becomes In (29a-29b) local electroneutrality (25) has been used to further simplify the expressions.

A. Current response and limiting conductances
We focus first on the current response obtained by numerically solving the system (22a-22c) under an applied potential difference ∆V . The two reservoirs are kept at the same ionic strength so that the only thermodynamic force at play is a constant electric field along the longitudinal coordinate. A positive (negative) ∆V corresponds respectively to the anode being placed at the left (right) reservoir. A standard measure of ionic rectification is given by the current-voltage (I-V) curve which we report in Fig. (2) for the case of λ D = 1/2 and Du = 1/2 and for different values of the channel slope. This regime corresponds to the case of partial Debye overlap inside the channel. For instance, in the case of k = 3/2 the local ratio λ D h(x) spans from ∼ 0.3 at the base junction up to ∼ 2 at the tip junction. Therefore by moving from left to right ions experience the building up of a Donnan potential passing from a region (left) where the bulk dominates to a region (right) where the EDL dominates. The non-linear curves in Fig. (2) display the usual diode-like behaviour reported in the literature, with a preferential direction of ionic current. When the electric field is applied parallel to the x−direction with the counterions moving from base to tip, the current is suppressed with respect to the Ohmic response (grey curve) and the system is said to be in a low conductance state. On the contrary, when the electric field is applied antiparallel with respect to x with the counterions moving from tip to base, the current is magnified and the system is said to be in a high conductance state. The rectification magnitude is monotonous in the degree of asymmetry in the system. The greater the channel's slope the larger the rectification. This must come as no surprise since the channel slope is the only element introducing asymmetry in the system. For k → 0 the channel is flat and it behaves like a standard Ohmic resistor. The numerical I-V curves can be compared with analytical predictions of the limiting differential conductances For strong Debye overlap the equations of motion reduce to (29a-29b). By neglecting the diffusive contribution to Dimensionless current I as a function of ∆V for a channel with λD = 1/2, Du = 1/2 at different value of channel slope, respectively κ = 0, 1, 3/2, 1.8. We recognize two different conductance states. For positive voltage drop (positive electric field) the system is in a low conductance state, the current being smaller than the Ohmic one (grey line). On the contrary for negative voltage drop (negative electric field) the current is magnified and the system is said to be in a high conductance state.
the mass flow with respect to the electrophoretic contribution in (29a) and by integrating in x we obtain Combining Eqs. (29a) and (29b) we solve for d x c in terms of the ratio I J 2hd x c + (2Du) 2 2hc d x log 2h = 2Du 2hc which is bound asymptotically, ∆V → ±∞, if the prefactor on the rhs vanishes, i.e. 2Du 2hc I J → 1. Accordingly the limiting conductance, G ±∞ , reduces to because the diffusive contribution to the ionic flux for very large fields is negligible. Eq. (33) implies that the marginal concentration inside the channel approaches a uniform value in the limit ∆V → ±∞. When λ D 1 we have analytical expressions for the marginal concentration at the channel's ends where, due to the channel geometry, the left end is characterized by the higher marginal concentration while the right end fixes the lower value. Hence, from Eqs. (27b-27b) (see Discussion section for further details) In Fig. (3) we show the I-V curves for λ D = 2 and Du = 1, i.e. in the regime of strong overlap. For k = 3/2 we report the analytical predictions for the asymptotic curves I ±∞ = ±G ±∞ ∆V with the limiting conductances obtained from (34), showing that these analytical expressions accurately capture the numerical results. Further discussion on the saturation mechanism for the conductance are reported in the Discussion section. From the comparison between Fig. (2) and Fig. (3) we observe that the quantitative structure of the I-V curves does not change respectively for partial Debye overlap with λ D = 1/2 and strong overlap with λ D = 2. It follows that the Debye length seems not to play a primary role in governing rectification. Notably, this is at odd with previous understanding of ICR which relies on λ D as the main controlling parameter.
In the next session this observation is further explored and clarified by looking closely to the dependence of ICR on the electrostatic lengthscales.

B. Current rectification ratio
In order to gain further insights on the rectified behaviour of the present system we introduce the rectification ratio η defined as the ratio between the absolute value of the current for opposite polarity of the external field. In the  Interestingly, η shows a strongly non-monotonic dependence on Du with a maximum of rectification approximately at Du ≈ 1/2. For Du 1 or Du 1 the rectification ratio goes to one and the standard ohmic behaviour is recovered. For values of Du close to unity the rectification ratio reaches a maximum which depends on the strength of the applied field upon reaching a saturation value as shown in Fig. (4). The saturation value of ∆V is itself modulated by Du. Fig. (5) shows that Du is a critical parameter controlling rectification, in contrast with λ D that seems not to be an adequate parameter to describe ICR. This is further illustrated by looking at Fig. (6), where η is plotted as a function of λ D for three different values of Du. We report a dashed line when we enter the regime in which linearization in ψ is no further justified. This happens in the limit of small λ D when the potential at the centerline vanishes and ψ ∼ ζ. In the regime of partial and strong overlap no significant dependence on λ D is shown. Albeit not quantitative, our results suggest that ICR decreases while approaching the limit of vanishing λ D . In this limit it is known that ICR approaches a non-zero asymptotic value 52 .

IV. DISCUSSION: THE ROLE OF DUKHIN NUMBER
IThe results of the previous section show that ICR is not primarily governed by the Debye length but rather by the Dunkhin length. This suggests that the Dukhin number directly controls the high (low) conductance state, for negative (positive) potential drop. This can be understood in terms of ionic concentration enrichment and depletion for opposite polarity of the external field, as discussed in previous works 23,32,51 . The panel in Fig. (7) shows the volumetric cross-sectionally averaged concentration c vol along the channel axis for two different regimes of λ D . In both figures we observe an overall increase (decrease) of ionic concentration for negative (positive) ∆V with respect to the equilibrium profile, represented by the grey line. Therefore the high conductance state for negative ∆V is due to an increase in ionic concentration inside the channel. The larger the external forcing, the stronger the accumulation of ions. On the contrary, when a positive voltage drop is applied the electrical conductance decreases due to the decrease of ionic concentration. In order to understand the phenomenon of salt accumulation and depletion we now turn our attention to the behaviour of the marginal concentration for large fields. In the previous section we already anticipated that in the limit of very large potential drop we expect the marginal concentration to saturate to a uniform value along the channel axis. hand, for positive potential drop the saturation value is bounded to the boundary condition at the right site (tip). Fig. (8) also shows an overshoot in the marginal concentration for large (but finite) negative ∆V . The overshoot is not present in the case of positive ∆V which stands as an additional sign of the asymmetry in the system. The microscopic mechanism causing it is still not clear and requires further investigations. Fig. (9) reports the marginal concentration profiles for increasing slope of the channel showing a significant dependence of the overshot on κ.
Asymptotically, Du 1, c L → c R , i.e. η → 1. In this regime transport is controlled by the diode surface, where entropic interactions are negligible with respect to electrostatic interactions and ions do not feel the symmetry breaking originated from the confinement. That is to say, enthalpy wins. The local marginal selectivity, γ ± (x) (directly proportional to the ionic marginal concentrations), constitutes a second, relevant quantity. For the counterions, the local selectivity at either end of the channel respectively reads making transparent the key role of the Dukhin number in controlling the local channel selectivity. Eq. (37) quantifies the relative importance of the counterion flux over the total transport. Due to the conical shape of the channel, γ R + is larger than γ L − , meaning that counterion transfer in presence of an external driving is larger at the tip than at the base. Such imbalance in selectivities results in a transient ion readjustment when an external driving is switched on. In the case of counterions moving from tip to base (negative ∆V ) this imbalance in selectivities results in a transient accumulation of ions inside the channel. On the contrary, when counterions move from base to tip (positive ∆V ) there will be a relative larger amount of ions leaving than entering the channel resulting in an overall decrease of salt concentration. In either case, the stationary state is reached when the nonequilibrium accumulation/depletion dynamics counterbalances the asymmetry of local selectivity induced by the geometry. Eq. (37) implies that the imbalance in selectivities is controlled by the asymmetry between Du/h L and Du/h R . Both Du 1 and Du 1 result in a uniform selectivity between the two ends of the channel, i.e. no rectification (see Fig. (10)). High Dukhin number, Du 1, means that the selectivity of counterions at either end tends to one (that is the selectivity of coions tends to zero): the coions are completely excluded from the system and the geometrical asymmetry is nullified by the perfect selectivity of the channel. No bulk transport is present so that the entirety of transport takes place in the EDL. On the other side, for Du 1 the selectivity at either ends tends to its bulk value 1/2. In this regime, irrespectively of the physical extension of the EDL the entirely of transport takes place in the unselective bulk and the omhic bulk response is restored. The asymmetry between Du L /h L and Du R /h R is maximized for Du ∼ 1 (in our case Du ∼ 1/2 because of the normalization used for the marginal concentrations). The qualitative interpretation of ICR caused by an asymmetry in the local selectivity at either end of the nanochannel is qualitatively consistent with the pioneer proposal of Woarmann 32 . However, our analysis provides a fresh interpretation of an old puzzle. We have shown that Du is the principal electrostatic parameter that locally controls the channel selectivity, with a secondary effect due to λ D , while Woermann pointed at λ D as the main length to be compared with channel confinement. Although it may fly against intuition, it is not the physical size of the EDL that determines the system capability to rectify ionic current.  10. The counterion selectivity at the tip γ R + (blu) and at the base γ L + (red) of a channel with k = 3/2 as a function of the reference Dukhin number. For Du ∼ 1/2 the difference between the two selectivities (yellow curve) is maximized, leading to a maximum of rectification.

V. CONCLUSION
In summary, we have presented here a theoretical analysis to address the phenomenon of ionic current rectification in nanometric channels. We have specifically focused on the case of a geometric ionic diode where the symmetry breaking is caused only by the conical geometry of the system. The theoretical framework mainly relies on two assumptions: a slowly varying channel geometry and a small electrostatic potential variation in the transversal direction. These ingredients allow us to derive formal expressions for the electrostatic potential, Eq. (19), and for the ionic current, Eq. (23), and to explore the response of the system for different values of λ D and Du. The main outcome of the work is the identification of the Dukhin length as the primary electrostatic length scale controlling rectification. It follows that rectification is expected to be measured in systems with size comparable to the Dukhin length, which remarkably can reach the micrometer scale 18 . This fact may explain recent experimental works 25,26 in which ICR is observed in mesoscopic pores. To conclude by misquoting Wolfrang Pauli 53 , it is a dynamical usage of surfaces that let the nanofluidic diode succeed where demons don't.

VI. APPENDIX
Here we report some details of the implementation in COMSOL for the numerical integration of the following equations: First of all let us recall here the appropriate boundary conditions for the system at hands. In both ends of the channel we have to impose continuity in the electrochemical potential for each species. Starting from the left side we write in adimensional variables: where by definition: By substituting eq.(40) into (39) we obtain: where the boundary condition for c ± (0) is expressed in terms of the function φ L (y) ≡ φ(0, y). The latter is obtained by solving the following transversal equation at x = 0: where φ L = 1 2h L h L 0 dyφ L (y) , h L = 1+ κ 2 and we made use of the fact that: Eq. (42) can be then numerically integrated using the standard electrostatic boundary conditions: Likewise we find the appropriate boundary value for c ± (1) using : Therefore the expressions (41) and (45) are now numbers which can be directly used as boundary conditions for the system in (38). It is also convenient in COMSOL to rescale the y variable in the following way: In this way we map the original domain to a square do-main substantially simplifying the COMSOL calculation. From the chain rule it follows: The only variables in the model that depend on y are ψ(x, y), φ L (y) and φ R (y). For each of them we apply (48) so that the electrostatic boundary condition for ψ (likewise for φ L and φ R ) become: ∂ y ψ (x, y = 0) = 0 (49) and the rescaled Poisson equation: