Local electroneutrality breakdown for electrolytes within varying-section nanopores

Abstract We determine the local charge dynamics of a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z-z$$\end{document}z-z electrolyte embedded in a varying-section channel. By means of an expansion based on the length scale separation between the axial and transverse direction of the channel, we derive closed formulas for the local excess charge for both, dielectric and conducting walls, in 2D (planar geometry) as well as in 3D (cylindrical geometry). Our results show that, even at equilibrium, the local charge electroneutrality is broken whenever the section of the channel is not homogeneous for both dielectric and conducting walls as well as for 2D and 3D channels. Interestingly, even within our expansion, the local excess charge in the fluid can be comparable to the net charge on the walls. We critically discuss the onset of such local electroneutrality breakdown in particular with respect to the correction that it induces on the effective free energy profile experienced by tracer ions. Graphical abstract

From a theoretical perspective, the usual approach is based on the Poisson-Boltzmann theory, within which the ions are regarded as point-like particles whose distribution fulfills the Boltzmann weight that is eventually determined by the solution of the Poisson equation.Such an approach can recover the dynamics of electrolytes confined between charged plates provided that the ions are sufficiently diluted and the surface charges are not too large.Within such conditions, the Poisson-Boltzmann theory recovers the global electroneutrality, namely the total charge in the liquid phase matches the one on the plates.However, recent experimental [15] and theoretical [16][17][18][19] results show that when a e-mail: p.malgaretti@fz-juelich.de (corresponding author) the plates are approaching at distances 10nm the reduced space and the steric and electrostatic interactions between the dissolved ions may lead to a breakdown of the global electroneutrality.This means that the net charge in the liquid between the plates does not balance the charge on the plates.However, since the electrolyte is typically in contact with a reservoir, the eventual global electroneutrality is attained when accounting for the charge distribution outside the channel [16,20].While global electroneutrality breakdown is associated with a significant increase in the total free energy of the system [21], local rearrangements of the charge may occur such that the electroneutrality is fulfilled globally but not locally, as it has been recently predicted for macroion solutions [22].
In particular, when the section of the confining vessel is not constant, novel dynamical regimes appear.Indeed, asymmetric pores have been used to rectify ionic currents [23][24][25], as well as to realize highly sensitive dopamine-responsive iontronic devices [26].Moreover, recirculation and local electroneutrality breakdown have been reported for electrolytes confined between corrugated walls [27,28] and the variation in channel section can tune their permeability [29,30] and even enhance the effective transport coefficients [31].
Accordingly, the question arises about the interplay between the geometry of the pore and the local electroneutrality breakdown.In order to be able to efficiently explore the parameter space we aim at an ana-lytical approach that can provide closed formulas capturing the dependence of electroneutrality breakdown on the parameters characterizing the system.Accordingly, we exploit the (linearized) Poisson-Boltzmann theory.Interestingly, our results show that even in the simplest scenario, i.e., at equilibrium and within the Debye-Hückel approximation, the interplay between the (varying) local section of the channel and the electrostatic forces leads to a breakdown of local electroneutrality and to the onset of an inhomogeneous excess charge for both dielectric and conducting channel walls in planar (2D) and cylindrical (3D) geometries.Interestingly, upon tuning the parameters, the local excess charge can be as large as the local charge on the walls.Once integrated along the transverse direction the local excess charge can be expressed via a multipole expansion whose leading term is a quadrupole.Finally, such an insight can be used to predict the corrections to the effective free energy profile experienced by a tracer ion induced by the local excess charge.

2D model
In the following, we analyze the case of a channel whose half-section h(x) (see Fig. 1) varies solely along the x direction.For later use, we introduce the entropic barrier [32][33][34] ΔS = ln h max h min . ( Here, h max ≡ h 0 (1 + h 1 ) and h min ≡ h 0 (1 − h 1 ) are, respectively, the maximum and minimum channel sections.The electrostatic potential φ inside the channel is determined by solving the (linearized) Poisson-Boltzmann equation, where k −1 = λ is the inverse Debye length.
In order to get an analytical insight into the solution of Eq. (3) we exploit the lubrication approximation, based on a separation between longitudinal and transverse length scales.In the system under study, the longitudinal length scale is captured by L which is the period of the channel section, whereas the transverse length scales are the channel average section h 0 and the Debye length λ.In typical situations we have that L λ and hence the Laplace operator can be simplified by assuming that the changes of φ along the longitudinal direction are much smaller than those along the transverse direction.Hence, we identify λ/L as the Fig. 1 Schematic view of the system: channel walls are in grey, the local surface charge is represented by blue crosses whereas the width of the Debye double layer is in red small parameter.Accordingly, we rewrite Eq. ( 3) as In the last expression, we have highlighted that the first term on the left-hand side is of O λ L 2 as compared to the second one.Accordingly, in the lubrication approximation, we treat the dependence on the longitudinal coordinate parametrically.The solution of Eq. ( 3) is governed by the boundary conditions on the channel walls that in the following are assumed to be the same on both walls.For conducting channel walls the potential equals the zeta potential ζ: Alternatively, for insulating channel walls, we have that where n is the unit vector normal at the channel walls and pointing inside the channel.
In particular, calling α = arctan(∂ x h(x)) the local slope of the channel walls and focusing on the upper channel wall (the same calculations can be re-derived for the lower one) we have: Accordingly, we get In order to determine the boundary condition at different orders in the lubrication expansion we recall that and hence we have that at leading order in lubrication the boundary condition reads whereas at the second order (i.e. next leading order) we have At equilibrium, we use the following ansatz for the linearized ionic number densities:

2D first order lubrication approximation
At leading order in the lubrication expansion, namely at order O(h 0 /L) 0 , the ionic number densities read ρ+,0 (x, y) = ρ +,0 (x) (1 − βezφ 0 (x, y)) , (13a) ρ−,0 (x, y) = ρ −,0 (x) (1 + βezφ 0 (x, y)) (13b) and the Poisson equation reads where we have introduced the notation The solution of Eq. ( 14) is then given by where is the, zeroth-order, inverse Debye length and A 0 is determined by the boundary conditions imposed by the channel walls.For conducting channel walls the potential equals the zeta potential ζ at the walls, φ 0 (x, h(x)) = ζ.Hence, we have Alternatively, for insulating channel walls and for smoothly varying-channels, ∂ x h(x) 1, the boundary condition can be expressed as [29] where we have substituted α ∂ x h(x) in Eq. ( 19).Accordingly, using Eq. ( 19) the electrostatic potential reads (20) At equilibrium, in order to solve for ρ +,0 (x), ρ −,0 (x), we impose that the electrochemical potential is constant [35]: Disregarding terms of order O(φ 0 ) 2 , by plugging Eq. ( 13) into Eq.( 21) leads to ρ ±,0 (x) = e β μ± (22) and the constant values of ρ ±,0 are set by the equilibrium chemical potentials μ± .For z − z electroneutral systems we have μ + = μ − and therefore This implies In particular, Eq. ( 24) leads to i.e., the Debye length is constant along the channel.Finally, the electrostatic potential for conducting channel walls reads whereas for insulating walls it is Substituting Eq. ( 27) or Eq. ( 28) into Eq.( 14) and using Eq. ( 25) the net local charge reads, respectively, which recovers the local electroneutrality of the system. 1 Equations ( 26)- (30) show that at leading order in lubrication, there is no correction to the Debye length, electrostatic potential or local charge and ionic density induced by the geometrical confinement.Such a tight relationship between the surface charge and the charge in the liquid phase is carved into the Debye-Hückel equation.In fact, multiplying by − and integrating the Debye-Hückel equation at first order in lubrication along the transverse coordinate leads to After integrating it follows that where we have identified that the right-hand side of Eq. ( 31) is indeed the total charge in the fluid and the left-hand side corresponds to minus the surface charge.Therefore, in order to check if local electroneutrality can be broken also at equilibrium, we need to account for higher order corrections in the lubrication approximation.

2D higher order corrections
It is clear that, at order O( h0 L ), Eq. ( 4) jointly with the boundary conditions, Eq. ( 19) leads to φ 1 (x, y) = 0. Therefore, the next leading order is O( h0 L ) 2 .Accordingly, the ansatz of the ionic number densities reads (33a) 1 We note that, at zeroth order in lubrication, the surface charge at each conducting wall can be obtained by At equilibrium the chemical potential is homogeneous.Therefore, since the zeroth order contribution to the chemical potential is already fulfilling the equilibrium conditions, see Eq. ( 21), we have that The expression for the second order contribution to the local chemical potential reads Using Eqs. ( 22),( 23),(33a),(33b) and keeping only terms linear in the electrostatic potential, Eq. ( 34) leads to from which we obtain Accordingly, the Poisson equation reads We remark that, due to Eq. (37), there is no second order correction to the Debye length, which indeed is clear from Eq. (38).The solution of Eq. (38), using Eq. ( 16), reads where A 0 (x) is the zeroth-order integration constant and it is determined by the following boundary conditions: for conducting walls and for dielectric walls.Finally, φ 2 (x, y) is determined by imposing the boundary conditions at the channel walls.
For conducting channel walls, at order O h0 L 2 , the boundary condition reads This leads to (44) In contrast, for dielectric channel walls the boundary conditions, at order O( h0 L ) 2 is given by Eq. ( 19) which leads to (similar results can be obtained for y = −h(x)) from which we obtain and accordingly we get

2D local electroneutrality breakdown
a. Dielectric channel walls Using Eq. ( 47) we can calculate the corrections to the local charge for dielectric channel walls 2 : For a dielectric channel, at second order in lubrication O(h 0 /L) 2 , the charge per unit area on the channel walls reads Accordingly, we can define the excess charge Using Eqs. ( 30),(48),(49), at second order in lubrication, the last expression reduces to This can be rewritten as Figure 2 shows the dependence of Δq σ on the longitudinal position.Interestingly, there is a depletion (with respect to the case of planar channel walls) of counterions close to the bottleneck, at x/L = 0, that leads to a net local positive charge.In order to fulfill global electroneutrality, this is compensated for by a net excess of counterions in the rest of the channel.As expected due to symmetry, the distribution of the excess charge is symmetric about the center of the channel and, in a far-field expansion, it leads to the onset of a net quadrupolar contribution.We remark that when the local excess charge Δq σ (x) is integrated over the channel period it leads to a vanishing contribution, i.e., global electroneutrality is retrieved.We quantify the magnitude of the local electroneutrality breakdown by computing the dimension-Fig.3 2D dielectric channel (σ > 0).Top left: absolute value of the quadrupole moment, as defined in Eq. ( 53), Q σ as a function of k0h0 for different values of ΔS as reported in the legend and with h0/L = 0.01.The grey dashed line is proportional to ∝ (k0h0) −1 and the dashed-dotted one to ∝ (k0h0) −2 .Top center: absolute value of the quadrupole moment, Q σ , as a function of ΔS for different magnitudes of kh0 as reported in the legend and with h0/L = 0.1.The thin grey dotted line is a guide for the eye and it is proportional to ΔS 2 .Top right: absolute value of the quadrupole momentum, Q σ , as a function of h0/L for different magnitudes of ΔS as reported in the legend and with kh0 = 1.The thin grey dotted line is a guide for the eye and it is proportional to (h0/L) 2 .Bottom: ratio ΔF σ 2 /Q σ for the same values of the parameters of the corresponding top panel less quadrupolar moment The top panel of Fig. 3 shows the dependence of Q σ on k 0 h 0 .Interestingly, Q σ shows a twofold scaling: for smaller values of k 0 h 0 , Q σ decays as k 0 h −2 0 , whereas for larger values of k 0 h 0 it decays as k 0 h −1 0 .In particular, the crossover between the two regimes is around k 0 h 0 1, i.e., when the Debye length is comparable to the average channel section.The central panel of Fig. 3 shows that Q σ grows almost linearly with the entropic barrier ΔS.Finally, the bottom panel of Fig. 3 shows that Q σ depends quadratically on h 0 /L.This is expected since these results were derived at the second order in lubrication.b.Conducting channel walls Using Eq. ( 44) we can calculate the corrections to the local charge for conducting channel walls 3 : (54) 3 We recall that For a conducting channel, we first define the effective surface charge as which up to second order contributions in lubrication reads as This can be rewritten as We recall that the last expression accounts for the surface charge density along the surface of the channel.However, when comparing the surface charge to that in the liquid phase we have to account for the fact that the latter is per unit length dx along the longitudinal axis of the channel.Accordingly, when computing the local electroneutrality we have to multiply Eq. ( 57) by the local area 1 + 1 2 (∂ x h(x)) 2 .We then get + sinh(k0h(x)) . (58) We define the excess charge as In this case, we chose a different normalization as compared to the dielectric case because at zeroth order in lubrication, the local surface charge is not homogeneous and this would lead to unphysical contribution when assessing the global electroneutrality.Finally, using Eqs.( 29), ( 54), ( 58) and (59) reduces to As for the dielectric case, we note that Eq. ( 60) shows the onset of local electroneutrality breakdown but also the fulfillment of global electroneutrality once Δq ζ is integrated over the channel period.By comparing Fig. 2 with Fig. 4 we note that while for the dielectric channel there is a clear and sharp peak at the channel bottleneck, for the conducting channel the peaks are located where the slope of the channel walls is maximum, i.e., x/L = ±0.25.At the same time, on the top of the shift of the maxima, we also observe that the magnitude of the peak is reduced in the case of conducting as compared to dielectric channel walls.This may be due to the fact that, even at zeroth order in lubrication, the local surface charge density on the channel walls is not constant: Eq. ( 29) indeed shows that, for k 0 h 0 1, the local charge density is minimum at the channel bottleneck.Accordingly, a smaller local surface charge leads to a minor departure from local electroneutrality.
For what concerns the magnitude of the local excess charge captured by the quadrupolar moment, the top panel of Fig. 5 shows that Q ζ decays as 1/k 0 h 0 for larger values of k 0 h 0 (as it is for dielectric channel walls) whereas for k 0 h 0 1, at variance with dielectric walls, Q ζ attains a plateau.Finally, the central and bottom panels of Fig. 5 show a behaviour similar to that observed for dielectric channel walls.

2D free energy barrier
While the leading term for the local electroneutrality breakdown is quadratic in h 0 /L this is not the case for the equilibrium free energy profile.For a tracer ion with elementary charge e, the local free energy can be obtained from the logarithm of the partition function and within the Debye-Hückel approximation it becomes 1 − βeφ(x, y)dy By expanding the logarithm, Eq. ( 62) can be decomposed into the following contributions: with where βF gas (x) is the free energy profile of an uncharged point particle, whereas βF 0 (x) and βF 2 (x) are, respectively, the leading order and higher order correction for In particular, we are interested in quantifying the contribution of the quadrupole moment to such a correction to the free energy barrier.The bottom rows of Figs. 3 and 5 report the ratio between the secondorder correction to the free energy difference and the quadrupole moment.As shown in the figures, the ratio of ΔF 2 and Q is generally sensitive to both k 0 h 0 and ΔS hence highlighting the relevance of higher-order multipoles in the free energy difference.Finally, as expected ΔF 2 and Q have the same scaling with h 0 /L and hence their ratio is insensitive to it.

3D first order lubrication approximation
Exploiting the experience gathered for the 2D case we can straightforward write down the solution of the Debye-Hückel equation in the case of axially symmetric channels which, at first order in lubrication, reads The solutions to it are where I n are modified Bessel functions of the first kind of order n.As for the 2D case, at linear order in the lubrication approximation the solution of the Debye-Hückel equation, Eq. ( 70), fulfills local electroneutrality: This indeed is minus the charge on the wall.

3D higher order corrections
At second order in lubrication the Debye-Hückel equation becomes whose solution reads with and B 2 is determined by the boundary conditions.

3D local electroneutrality breakdown
a. Dielectric channel walls For dielectric channel walls the boundary condition is the same as the one derived for the 2D case, Eq. ( 19), and reads Hence, Accordingly, the second-order correction to the integrated local charge can be obtained by substituting Eq. (78) into Eq.(70) and integrating along the radial direction (see "Appendix A"): Finally, recalling that, at second order, the surface charge per unit longitudinal length is the local charge excess, at second order in lubrication, reads (see "Appendix B") Figure 6 shows the dependence of Δq on the position.Interestingly, as for the 2D case, Δq displays a maximum at the channel bottleneck, x/L = 0.5.However, a part of the different shape of the profile, the main striking difference between Figs. 2 and 6 is the difference in magnitude of the effect.In fact, while for the 2D case the excess charge at the peak is comparable to the net charge density on the wall, for the 3D case (with similar geometry, ΔS, and Debye length, k 0 h 0 ) the effect is weaker.Hence we do expect that local electroneutrality breakdown to be more prominent for slablike channels then for cylindrical pores.More in detail, Fig. 7 shows again a similar trend as compared to the 2D case, Fig. 3, the only major difference being the nonmonotonous dependence on ΔS shown for large values of k 0 h 0 .b. Conducting channel walls For conducting channel walls the boundary condition is the same as the one derived for the 2D case, namely φ ζ 2 (x, h(x)) = 0, and hence Eq. (74) leads to  The second order correction to the wall charge can be obtained as in the 2D case, Eq. ( 56), where we change y for r and we multiply by the local area 2πh(x)(1 + 1/2(∂ x h(x)) 2 ): Similarly to the 2D case, the net charge in the liquid phase reads Combining Eqs. ( 71), ( 83) and (84) leads to the local excess charge density as shown in Fig. 8,9

3D free energy barrier
Similarly to the 2D case, at equilibrium, the local free energy of a tracer ion with elementary charge e is given by βF Within the Debye-Hückel approximation and expanding the logarithm, Eq. ( 62) can be decomposed into the following contributions: βF DH (x) βF gas (x) + βF 0 (x) + βF 2 (x), (87) with where βF gas (x) is the free energy profile of an uncharged point particle whereas βF 0 (x) and βF 2 (x) are, respectively, the leading order and higher order correction for charged particles.In order to assess the impact of the local excess charge on the dynamics of a tracer ion, the bottom rows of Figs. 3 and 5 report the ratio between the second-order correction to the free energy difference and the quadrupole moment.As already mentioned for the quadrupole, the overall behaviour of ΔF 2 resembles the one observed in the respective 2D cases.

Conclusions
In this contribution, we focus on the case of an electrolyte embedded between corrugated channel walls.In order to gain analytical insight we restrict our analysis to channels whose section is varying smoothly enough so that we can exploit the lubrication approximation to solve for the Poisson equation.Under such approximation, we have derived closed formulas for the corrections induced by the varying section of the channel to the local charge distribution.In particular, at equilibrium, the Debye length keeps homogeneous even when second-order corrections in the lubrication expansion are accounted for.At variance, while at first order in lubrication the local electroneutrality of the system is recovered, this is not so for second order corrections.This implies that upon reducing the length scale separation between the longitudinal and transverse direction the local excess charge will grow and this will also induce additional corrections to the effective free energy profile experienced by a tracer ion.Such local charge reorganization within corrugated channels has been observed so far only in out of equilibrium situations [28] where the advection of the ions plays a major role.Our results show that such a phenomenon occurs also at equilibrium and hence solely due to the interplay between the geometry of the channel and the electrostatic forces.In order to assess the robustness of our results we have derived such corrections for both dielectric and conducting channel walls in both planar (2D) and cylindrical (3D) geometries.Indeed, our results show quite a remarkable similarity between the 2D and 3D cases for both dielectric and conducting channel walls.Interestingly, the dielectric case shows an enhanced sensitivity to the dimensionality, as compared to the conducting case.In particular, the local (integrated) excess charge attains its maximum at the channel bottleneck (x/L = 0), for 2D dielectric walls, and it can be as large as the bare charge (density) on the walls.This is not the case for cylindrical channels (3D) for which the local excess charge is 100 times smaller than the local charge.The difference in the location of the excess charge along the channel for dielectric and conducting walls indicates that the specific boundary conditions play a relevant role in the transport properties of confined electrolytes and ions.At variance, for conducting channels, the maxima of the local excess charge are located where the slope of the channel walls is maximum (x/L = ±0.25)for both 2D and 3D cases.All in all, the magnitude of the corrections that we report on are not very large and indeed this is expected since they are obtained via an expansion.However, this may not be the case when the longitudinal length is comparable to the Debye length.For typical values of the Debye length, λ 10 − 100nm, this implies to have channels with length L 100nm.This is the case for many biological ionic channels as well as synthetic pores and membranes.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Fig. 2
Fig. 2 2D dielectric channel (σ > 0).Local excess charge Δq(x) as a function of the normalized position x/L along the channel axis for different values of ΔS as detailed in the legend with h0/L = 0.1, and kh0 = 1

Fig. 4
Fig. 4 2D conducting channel (ζ > 0).Local excess charge Δq(x) as a function of the normalized position x/L along the channel axis for different values of ΔS as deailed in the legend with h0/L = 0.1 and kh0 = 1

Fig. 5
Fig.52D conducting channel (ζ > 0).Top left: absolute value of the quadrupole moment as defined in Eq. (53), Q σ as a function of k0h0 for different values of ΔS as reported in the legend and with h0/L = 0.01.The grey dashed line is proportional to ∝ (k0h0) −1 and the dashed dotted to ∝ (k0h0) −2 .Top center: absolute value of the quadrupole moment, Q σ , as a function of ΔS for different magnitudes of kh0 as reported in the legend and with h0/L = 0.1.The

Fig. 6
Fig. 6 3D dielectric channel (σ > 0).Local excess charge Δq(x) as a function of the normalized position x/L along the channel axis for different values of ΔS as deailed in the legend with h0/L = 0.1 and kh0 = 1

Fig. 8
Fig. 8 3D conducting channel (ζ > 0).Local excess charge Δq(x) as a function of the normalized position x/L along the channel axis for different values of ΔS as detailed in the legend with h0/L = 0.1 and kh0 = 1

Fig. 9
Fig. 9 3D conducting channel (ζ > 0).Top: local excess charge Δq σ as a function of k0h0 for different values of ΔS as reported in the legend and with h0/L = 0.01.The grey dashed line is proportional to ∝ (k0h0) −1 and the dashed dotted to ∝ (k0h0) −2 .Center: local excess charge Δq σ as a function of ΔS for different magnitudes of kh0 as reported 3D dielectric channel (σ > 0).Top left: local excess charge Δq σ as a function of k0h0 for different values of ΔS as reported in the legend and with h0/L = 0.01.The grey dashed line is proportional to ∝ (k0h0) −1 and the dashed dotted to ∝ (k0h0) −2 .Top center: local excess charge Δq σ as a function of ΔS for different magnitudes of kh0 as reported in the legend and with h0/L = 0.1.Top right: local excess charge Δq σ as a function of h0/L for different magnitudes of ΔS as reported in the legend and with kh0 = 1.The thin grey dotted line is a guide for the eye and it is proportional to (h0/L) 2 .Bottom: ratio ΔF σ 2 /Q σ for the same values of the parameters of the corresponding top panel