Charge polarization, local electroneutrality breakdown and eddy formation due to electroosmosis in varying-section channels

We characterize the dynamics of an electrolyte embedded in a varying-section channel under the action of a constant external electrostatic field. By means of molecular dynamics simulations we determine the stationary density, charge and velocity profiles of the electrolyte. Our results show that when the Debye length is comparable to the width of the channel bottlenecks a concentration polarization along with two eddies sets inside the channel. Interestingly, upon increasing the external field, local electroneutrality breaks down and charge polarization sets leading to the onset of net dipolar field. This novel scenario, that cannot be captured by the standard approaches based on local electroneutrality, opens the route for the realization of novel micro and nano-fluidic devices.

The transport of ions, molecules and polymers across constrictions such as pores, membranes or varyingsection micro-nano-channels is crucial for several biological as well as synthetic systems. For example, in biological cells ion channels control the uptake of ions from the environment [1] whereas, in resistive pulse sensing techniques, the interactions of colloidal particles [2] or macromolecules [3][4][5][6][7] with the nano-or micro-pore are measured from the variation of the electric conductance of the pore induced by the presence of the particle. Moreover, electro-osmotic flux can play a crucial role in molecule capturing in nanopores [8][9][10]. Recent studies have shown that ionic transport and electro-osmotic flow in micro-and nano-fluidic circuitry can be controlled by tuning the geometry of the micro-nano-channel [11][12][13]. Indeed, conical [14][15][16] or heterogeneously charged [17] pores have been used to realize nano-fluidic diodes and to rectify electro-osmotic flows [18,19]. Moreover, periodic varying-section channels have been used to realize nanofluidic transistors [20]. Similarly, recent contributions have shown that the shape of the confining vessel and the boundary conditions therein imposed can be exploited to control the flow. Indeed, electro-osmotic transport can be strongly enhanced by grafting charged brushes on the channel walls [21] and by hydrophobic surfaces [22].
In this contribution, we study, via molecular dynamics simulations the electro-osmotic flow of an electrolyte embedded in a varying-section channel. The advantage of our approach, as compared to others based on the solution of some continuum models [23][24][25][26][27], is that the ionic densities are left free to relax according to the interactions among ions and between ions and the channel walls. Therefore our approach allows us to critically discuss, for example, the local electroneutrality assumption and its possible breakdown. Our results show that Lx . Periodic boundary conditions are applied in the x and y direction. B) Snapshot of molecular dynamics set-up. The solid walls (gray) are constituted by Lennard-Jones atoms. The atoms of the layer exposed to the liquid are charged (blue). Water molecules are not reported while blue and red spheres represent K + and Cl − ions. An external electrical field parallel to the x-axes in applied.
when the systems is driven by a constant electrostatic field acting along the longitudinal axis of the channel, inhomogeneous ionic and charge densities are induced due to the variations of local channel section. This phenomenon is similar to concentration polarization (CP) reported for electrolytes transported across ionic-selective membranes [23,25,28] (see also Ref. [29,30] for recent reviews) i.e., in open circuit conditions and under severe modulations of channel section [25,27,28]. Interestingly, our results show that CP can be obtained also for close circuit conditions and for smooth variations of channel section. In particular, we observe CP when the Debye length is comparable with the channel bottlenecks, i.e., in the entropic electrokinetic regime [31]. Concerning local electroneutrality, our numerical results confirm that even in the entropic electrokinetic regime, the local elec- troneutrality assumed by previous works [26,27,30,32] is fulfilled for mild values of the external field. However, upon increasing the external field our numerical simulations show that local electroneutrality breaks down and charge polarization (QP) sets, leading to the onset of a net dipolar contribution to the electrostatic field. Interestingly, a similar phenomenon has been recently observed for pressure-driven flows across conical pores [33].
Results We study the electro-osmotic flow (see Suppl. Mat. for the details) of a KCl water solution across a channel with half section with average section h 0 = 38.5Å and modulation h 1 = 32.5Å and periodicity L x = 836.6Å, see Fig. 1. Channel section is constant along the y direction, with thickness L y = 60.37Å, and periodic boundary conditions are applied along x and y. Both channel walls are covered with a constant charge density σ = 0.15 C m 2 . The system is globally electroneutral since we compensated the wall charge with additional Cl − ions. After equilibration, a homogeneous and constant external electric field E = (E x , 0, 0) is applied to the whole system.
We begin our analysis by focusing on the case of larger ionic concentration ρ 0 3M . The density profile of both Cl − and K + is expected to decay over a length scale comparable to the Debye length, λ = ε β(ze) 2 ρ0 , where e is the elementary charge, z is the valence of the ions (z = 1, in our case), β = (k B T ) −1 with k B the Boltzmann constant and T the absolute temperature and ε = 80 · ε 0 is the dielectric constant with ε 0 the vacuum dielectric constant. For 3M ionic solution we estimate a Debye length of λ 1.8Å and therefore λ h min being h min = h 0 − h 1 = 6Å the half-section calculated at channel bottlenecks. In such a case, since there is only a small overlap between the Debye layers of the two facing walls, we do not expect the onset of any entropic electrokinetic effects [31]. As expected, the accumulation of Cl − ions at the positively charged solid wall induces and electro-osmotic flow, opposed to the direction of the external electrostatic field, that is almost symmetric with respect to the z−axis (see panels a)-b) of Fig. 2).
We have then increased the Debye length by reducing the ionic concentration ρ 0 . For ρ 0 ∼ 1M we have λ 3.1Å. Hence, the Debye layers of the facing walls overlaps in the narrower sections (λ/h min ∼ 1, h min = 6Å) highlighting that the system is within the entropic electrokinetic regime [31]. Panel b) of Fig. 2 shows that, for ∆V = 2.4 [V], the ionic densities are quite affected by the flow. In particular, Cl − concentration largely increases in the channel bottleneck, see e.g. the region x/L x ∈ (0.45, 0.5) in Fig. 2b. This increase in Cl − concentration is associated to K + depletion, that will be discussed more in details in next paragraphs. This feature is associated to the onset of eddies in the electro-osmotic velocity profile, as shown in panels c)-d) of Fig. 2. These eddies form for sufficiently large driving forces. Indeed, -- while for ∆V = 0.6 [V] we do not observe major discrepancies with the previous case, for ∆V = 2.4 [V] two eddies form inside the channel. Interestingly, such eddies break down the left-right symmetry of the channel, for instance, Fig. 2d, shows that the eddies are shifted in the direction of the volumetric fluid flow, i.e. negative x in our reference frame. This occurrence is in contrast to the prediction obtained in linear regime [31] for which the eddies center is in the channel center, x/L x = 0. Then we further increase the Debye length by setting ρ 0 ∼ 0.3M , for which λ 5.6Å . In such a regime panels e)-f) of Fig.2 show that the zone in which K + are depleted is enhanced as compared to the previous cases. Moreover, comparing panels e)-f) to c)-d) in Fig.2 we notice that, for ρ 0 ∼ 0.3M , the onset of the eddies occurs for smaller values of the external force as compared to ρ 0 ∼ 1M .
In order to quantitatively capture the accumulation of ions density, we have analyzed the dependence of the ionic densities averaged over the channel section as a function of the longitudinal position. Fig.3 shows that the dependence of normalized densities profilesρ ± (x) de-fined asρ on the longitudinal position strongly depends on the value of the Debye length as compared to the channel section at the bottleneck. Indeed, when λ/h min 1, i.e. for ρ 0 = 3M , panel a) of Fig.3 shows that the density profiles are almost constant along the channel for all the values of the external force we have tested. The only relevant difference occur close to the narrowest sections x/L x ±0.5. In contrast, when λ/h min ∼ 1, i.e. for ρ 0 = 1M and ρ 0 = 0.3M panels b)-c) of Fig.3 show the onset of a region where the coions (K + ) are depleted and the left-right symmetry is broken. Once we have analyzed the density of Cl − and K + separately we move to the local total charge and local ionic densities defined as: We stress that we choose a different normalization for q(x) andρ(x) so that in a plane channel |q(x)| always match the total surface charge of channel walls 2|σ| whereasρ(x) is the average density that, in the Debye-Hückel regime, determines the local value of the Debye length 1 . Interestingly, Fig. 4.a shows that in the limit of small Debye lengths, λ/h min < 1, local electroneutrality is recovered, i.e. q(x) 2σ, and the ionic density profile is almost constant along the channel. Only at very high voltages a small deviation is present at x/L x 0.4. In contrast, for larger values of the Debye length, λ/h min 1, neither q(x) norρ(x) are constant. Such inhomogeneities trigger the onset of modulations in the magnitude of the local electrostatic potential. In particular, an accurate inspection of panels b)-c) of Fig.4 reveals that weaker values of the external field triggers solely an inhomogeneity in the ionic density but do not affect the local electroneutrality. Interestingly, upon increasing the strength of external electric field local electroneutrality breaks down and on the top of the well-known concentration polarization (CP) [26,27,30,32,34], a charge polarization (QP) appears. We remark that charge polarization sets for smaller values of the external force for smaller values of the Debye length. This can be due to finite liquid slippage at solid wall that cannot be disregarded in the regime under study. Indeed, slippage 1 We stress that according to Eq.(4) we have that  is commonly described in term of the Navier boundary condition that, for a plane channel, reads u t | wall = L s ∂u t /∂n with u t the component of the velocity field tangent to the wall, n the normal to the wall and L s the slip length [35]. Atomistic simulation showed that for smooth hydrophilic and slightly hydrophobic (contact angle < 120 • ) surfaces L s , hardly exceeds a nanometer [36][37][38]. In addition, the presence of a strong surface charge, further reduces the slip length for hydrophobic surfaces [39]. When comparing simulations performed with different ionic concentrations it should be taken into account that the slip length L s in the three setups may be slightly different. Indeed, the relevant parameter ruling the effect of slippage on the electroosmotic flow is the ratio L s /λ D [40]. This feature is emphasized by Fig.S1 in Suppl. Mat. that shows that the mismatch between the prediction of the analytical model (see Ref. [31]) and the numerical results increases upon decreasing λ D .
In conclusion, we have reported on numerical simulations concerning a KCl solution embedded in a varyingsection channel under the action of a constant electrostatic field. Our simulations show that, when the Debye length is comparable to the width of the channel bottlenecks, the system is in the entropic electrokinetic regime that is characterized by the onset of eddies [31]. In this perspective we observe, in agreement with what has been reported in the literature, the onset of a concentration polarization and local recirculation of the fluid velocity that comes along with the onset of a standing shock in the ionic concentration. Surprisingly, for stronger external fields the local electroneutrality breaks down and an additional charge polarization (QP) sets in. Such a novel phenomena has been observed thanks to our microscopic approach based on Molecular Dynamics simulations in which the ionic densities are not constrained. In this perspective, our results show that for mild external fields local electroneutrality is recovered. This can justify a posteriori the assumption of local electroneutrality in these regimes. However, for larger external fields, local electroneutrality does not hold and a net electric dipole sets in inside the channel.
Acknowledgments PM acknowledges Dr. Mathijs Janssen for useful discussions.
This research used the computational resource from CINECA (NATWE project), and the Swiss National Super-computing Centre (CSCS), project ID sm11.

Appendix A: Ion and electroosmotic currents
Here, we discuss the dependence of ionic and electroosmotic currents on the applied voltage. Figures 5.a and  5.b show the electro-osmotic flow Q and the ionic current I as a function of the applied electrostatic potential drop ∆V = L x E x along with the predictions of Poisson-Nernst-Plank (PNP) theory presented in [31]. Both panels of Fig 5 show a significant mismatch between the numerical results (symbols) and the PNP model (dashed lines). Interestingly the mismatch increases upon increasing the salt concentration, i.e. decreasing the Debye length. As discussed in the main text, the effect of the slip on the walls of the channel increases as the ratio L s /λ D where L s is the so-called slip length and it captures the magnitude of the slip. In this perspective, Fig.5.a makes us speculate that the mismatch between the numerical results and the PNP model is due to the partial slip of the fluid on the channel walls. Indeed, for high concentration (low λ D ), the MD electro-osmotic flux (blue triangles) is significantly larger than the PNP prediction (dashed blue line). The partial slip can also explain the mismatch shown in Fig.5.b. While in a flat channel an enhancement of the volumetric fluid flow Q leads to an increased electric current (see Ref. [41]) the situation is different in the case of a corrugated channel. Indeed, in the latter case the K + ions have to overcome a potential barrier at the channel bottleneck whereas Cl − can move along the channel walls at almost constant potential. Therefore the enhanced Q will reduce the magnitude of the effective barrier that K + have to overcome. Since the flow is quite sensitive to the high of the barrier (for very large barrier the time scales exponentially with the high of the barrier) the enhancement of Q will be more beneficial to the flow of K + as compared to the gain for Cl − . Since the electric current is defined as I = J K − J Cl the enhancement of Q leads to a hampered electric current, i.e. the electric currents from MD (symbols) are smaller than the PNP predictions (dashed lines), in particular for larger ionic concentrations.

Appendix B: Methods
System set-up and equilibration. The system is constituted by a KCl water solution confined by two solid walls. The solid walls are formed by atoms arranged in a face centered cubic (fcc) structure with lattice side a = 4.92Å. The mass of each atom is m w = 70 a.m.u. and they interact via Lennard-Jones potential with potential well depth w = 0.152 Kcal/mol) and van der Waals radius σ w = 3.15Å. The curved walls are obtained by solid slabs cut parallel to 111 planes of the fcc solid. The two slabs are then bend to follow the prescribed channel shape with a proper change of the z-coordinate. Each atom exposed to the liquid has a charge q w .
Water molecules and ions at a molar concentration ρ 0 are added using VMD [42]. In particular, indicated with N Cl − and N K + the number of Chloride and Potassium ions, we have N Cl − = 1749 and N K + = 789 for ρ 0 = 0.3M , N Cl − = 3327 and N K + = 2367 for ρ 0 = 1M , and N Cl − = 8062 and N K + = 7102 for ρ 0 = 3M . In each case, N Cl − − N K + = 960. Since 960e is the total charge on the walls, the total charge of the system is zero.
All MD simulations are performed using NAMD software [43] by implementing periodic boundary conditions being L x = 836.6Å, L y = 60.37Å and L z = 200Å the box sizes. L x and L y are dictated by the size of solid walls while L z is arbitrary set to 200Å to avoid nonbonded interactions among the different images along the z-axes. The total number of atoms is ∼400 000. A snapshot of the system is reported in figure 1.b. Particle mesh Ewalds (PME) summation method was employed for the electrostatics [44]. TIP3P model [45] was employed for water while CHARMM36 force field [46,47] with NBFIX corrections is used for the ions [48]. More accurate water models are available, see e.g. the TIP4P 2005 [49] often used for nanofluidics [50,51]. However, these models are much more computational demanding than TIP3P, hence, in this work, beside the well known limitations of TIP3P (e.g. the viscosity is lower than the experimental value) we preferred TIP3P since it allowed us to run longer simulations and to partially reduce the statistical error.
Walls atoms are constrained at their initial position by harmonic springs (spring constant k = 10 Kcal/(mol A 2 )). A first 0.8 ns NVT simulations (T = 310K, time step 2 fs) was run and pressure is estimated from the average force acting on the walls. The distance between the upper and lower wall is then gradually reduced and further 0.8 ns simulations are run until the pressure reached a value close to 1 atm.
Production runs. As common in the MD studies on ionic and electroosmotic transport, a homogeneous and constant external electric field E = (E x , 0, 0) is applied to the whole system [52][53][54]. Snapshots are acquired every ∆t = 2ps.
The water velocity field in the x,z plane is calculated as follows. The velocity v i of the i-th water molecule is estimated as v i = [x i (t + ∆t) − x i (t)]/∆t where x i is the position of the oxygen atom of the i-th water molecule. The mean velocity of each particle v i is associated to the pointx i = [x i (t+∆t)+x i (t)]/2. The value of the velocity field u in a given point x is hence calculated as the mean of the velocities of the particles in the neighborhood of x, i.e. the molecule for whichx i belongs to a region [x − ∆x/2, x + ∆x/2] ∩ [z − ∆z/2, z + ∆z/2]. The production runs span in a range of 40−90 ns with longer simulations used for smaller forcing and lower ionic concentration.