ELECTROSTATIC CORRELATIONS AND NON-ELECTROSTATIC EFFECTS ON ION DYNAMICS IN ALTERNATING CURRENT VOLTAGES

RESUMO – We present a modified Poisson-Nernst-Planck (mPNP) model to include both non-electrostatic effects and electrostatic correlations on ion dynamics near charged parallel blocking electrodes. The first contribution appears directly on ion chemical potential, whereas the latter is addressed by a modified Poisson equation. The model consists in a partial-differential-algebraic equation system, solved by the DASSLC code, implemented in EMSO simulator. Monte Carlo simulation results were properly described without the aid of fitting parameters. The model predicts overscreening for 2:1 electrolytes, which was not accessed by classical PNP-like approaches. Alternate voltage simulations show a meaningful effect of electrostatic correlations on mobile charge dynamics for 2:1 electrolytes, where ion-ion correlations are stronger near the electrodes. This is the first report accounting for electrostatic correlations on ion dynamics via PNPlike models under alternating current (ac) voltages to the best of our knowledge.


INTRODUCTION
Ion dynamics under the influence of external electric fields is a fundamental feature common to all electrokinetic modeling strategy.Using mean-field approaches, this problem could be addressed through Poisson-Nernst-Planck (PNP) type models, in which the focus is on mobile charge dynamics subjected to external electric fields, without considering bulk fluid flow.Most of the limitations of classical theories are related to this part of the problem, since the ions are usually treated as point-like charges, with no dispersion interactions, no hydration (solvent) effects, and no ion-ion correlations.Whereas electrokinetic flow is usually laminar (low Reynolds number), in which numerical procedure for solving Navier-Stokes equation is standard, the difficulties in treat the complete electrokinetic phenomenon are frequently related to the complex coupling with PNP equations.
The limitations of classical PNP model for ion transport modeling are well known.In addition to those inherent to mean-field approaches (the solvent is a continuum without any local correlation), the ions are point-like charges, without fluctuations due to ion-ion size and electrostatic correlations, and also dispersion (van der Waals) interactions between each ion and the electrode are not considered (Alijó et al., 2014a).Although modified PNP (mPNP) derived by the inclusion of size correlation effects have experienced some success, such models fail to describe short-range Coulomb correlations, the so called electrostatic correlations.In many relevant situations, such as concentrated and/or multivalent electrolytes, ionic liquids, molten salts, or plasma systems, ion-ion correlations may appear, leading to overscreening of a charged surface, which is the formation of subsequent layers that overcompensates the charge of the previous layer, and so on (Bazant et al., 2011;Storey e Bazant, 2012).
The main purpose of this work is to analyze the effects of electrostatic correlations on ion dynamics in ac voltage blocking electrodes.Thus, we have proposed a mPNP model that accounts for both ion size effects and electrostatic correlations in a thermodynamically consistently way, enabling to account for ion size asymmetries as well as other non-electrostatic contributions.We show that electrostatic correlation improves the description of available MC simulation data of Boda et al. (2002), especially for multivalent ions in which ion-ion correlations are stronger, capturing ion density oscillations in a simple continuum framework without fitting parameters.Moreover, ac voltage simulations show that electrostatic correlations must be accounted for a complete description of mobile charge dynamics, even at moderate to low voltage gaps.This is the first report accounting for electrostatic correlations on ion dynamics via PNP-like models under ac voltages to the best of our knowledge.

METHODOLOGY
Ion dynamics is considered here in a simple picture sketched in Fig. 1: parallel-plate electrodes separated by L and immersed in electrolyte solutions suddenly subjected to a step direct current (dc) voltage (Fig. 1a), or to an ac voltage (Fig. 1b).The fundamental difference between these two cases is in system input disturbance.In the former we apply a dc step voltage of 0 s VV  in both electrodes at the initial time, giving an overall potential drop across the cell of 0 2V .In the latter, the harmonic is applied to the electrodes, so the overall potential drop across the cell is Vt  , corresponding to a 0 4V peak-to-peak voltage.The solvent is represented only by the dielectric permittivity  in both cases.We consider "blocking electrodes", i.e., no Faradaic reactions occur (no ion flux through the electrodes).
Here we skip the details of model derivation and show only the final formulation of mPNP model equations.For a complete description, we refer the reader to our previous work on steric effects in dc voltages (Alijó et al., 2014a), and to the manuscript that we recently submitted to Langmuir for publication (Alijó et al., 2014b).
The PNP model is a PDAE system involving variables and parameters in several orders of length and time scales, which may complicate the numerical solution.These aspects encourage problem re-parameterization, in order to give non-dimensional equations with gain in terms of smoothly and accuracy of the numerical solution.Therefore, we define the following dimensionless variables for the spatial coordinate is a characteristic length used to rescale some variables, where I is the ionic strength, related with bulk ion concentrations   As a first approximation, and also to be consistent with the Monte Carlo simulation results used to adjust the model, the ions are supposed to be equally sized ( i   , for all ions).Because of that, ion diffusion coefficients i D are assumed to be the same and equal to 0 D for all ions i. Applying the rescaling factors listed above to mPNP model equations, we obtain the following PDAE system (Alijó et al., 2014b): , for 1,..., ,0 1 ln ln , for 1,..., The mPNP model proposed here needs   boundary conditions (BCs) and c n initial conditions (ICs) to be solved.Since we are dealing with blocking electrodes, we have no ion flux through them: 01 0, for 1,..., Other two BCs come from (no correlations at the electrode interface), by applying the re-parameterizations proposed here and differentiating  with respect to the spatial independent variable  , giving: 01 0, which is equivalent to written in terms of y and  , i.e., the third derivative (with respect to  ) of the electric potential y is null at the electrode surfaces ( 0   and 1
The last two BCs are related to the electric potential at electrode surfaces.We follow a recent work (Olesen et al., 2010) and emulate a compact layer imposing a thin dielectric coating between the electrode and the electrolyte, which leads to mixed boundary condition at the electrodes: where  is a parameter that characterizes the ratio of the compact-layer capacitance to that of the diffuse layer in low voltage limit (Olesen et al., 2010).In the limit 0   , this compact layer vanishes, and a diffuse layer occupies the full domain between the electrodes.
The external applied voltage     0 sin s yy    is generally an oscillating time function with frequency  and amplitude 0 y .This is the case of ac voltages applied to the electrodes (Fig. 1b), which is the main focus of this work.In contrast, for dc voltages (Fig. 1a), the external electric potential applied is constant, i.e.,   0 s yy   .We assume that all ions have bulk concentrations at 0   .Therefore, initially all ions are uniformly distributed, as if no electric field is applied to the electrodes, i.e.,   0 1 i     , for 1,..., c in  .
The mPNP model described in this section is numerically integrated in time after spatial approximation via polynomial interpolation in finite elements, using the Jacobi polynomial, ( , ) () The main advantage of this approach, rather than a classic finite-difference scheme, is the accuracy with much less discretization points, without facing oscillatory instabilities of classical orthogonal colocation methods.Problem domain has been divided into a non-uniform grid that is denser in the vicinity of the electrodes and sparser towards the symmetry axis.Mesh convergence depends on the applied voltage, ion valences, and salinities.For the conditions studied here, mesh convergence is ensured with at most 26 finite elements with 1 n  internal point and Jacobi weights parameters 1  ab .The time derivatives are then integrated using DASSLC code (Secchi, 2010), implemented in the dynamic simulator EMSO (Environment of Modeling, Simulation and Optimization) (Soares e Secchi, 2003).In this work, we account for the size effects using the lattice-based Bikerman's model (Bikerman, 1942).The correspondent excess chemical potentials ( , ex Bik i


) is given elsewhere (Storey e Bazant, 2012).This is directly included on the chemical potential of each ion i in solution as a nonelectrostatic term, that is , which represent the non-electrostatic contribution taken into account here.

Comparison with Monte Carlo simulation results
Before presenting results for oscillatory electric fields, we apply a dc voltage to a single plate electrode and observe each ion concentration profile after all transients have decayed (stationary condition).The main purpose here is to adjust a single mPNP-model parameter to give a better description of MC simulation data for several electrolyte solutions, i.e., 1:1 and 2:1 salts.Table 1 list the values of some constants used in the simulations presented in this section.In Fig. 2 we compare ion distributions ( ,0 ) calculated by our mPNP model with those predicted by Monte Carlo simulations of Boda et al. (2002) (spheres) for 1:1 (Fig. 2a) and 2:1 salts (Fig. 2b) with (blue continuous line) or without (black dashed line) electrostatic correlations.
As observed by Storey e Bazant (2012), we find that mPNP offers a good description of ion distributions for 2:1 electrolytes, comparing with Monte Carlo simulation results, only when the electrostatic correlations are accounted, since classical continuum models cannot capture the essential physics of overscreening, and are limited to the description of ion density profiles which decay monotonically.Better agreement is achieved with an electrostatic correlation length equal to the ionic radius ( 0.5 ), providing a qualitative description of the double-layer with no fitting parameter.Moreover, as expected for simple 1:1 salts, we observe a modest effect of ion-ion correlations in Fig. 2a.However, the electrostatic correlations improve the description of MC data.
Table 1 -Specified constants used in mPNP analysis.

Alternating current (ac) voltage effects
In order to properly observe electrostatic correlation, we present dynamic profiles for 2:1 electrolytes, subjected to ac voltage fields.We show surface electric potential (y) and cation (  Therefore, we show that electrostatic correlations must be considered in order to properly understand and describe ion dynamics under ac voltages, especially for multivalent electrolytes, in which ion-ion correlations may be strong.The Poisson-Nernst-Planck framework presented here, as well as other mean-field models with different methodologies to describe non-electrostatic effects, have a wide potential of application, with little additional effort in terms of computational cost for the numerical solution.Finally, we should emphasize that the methodology presented here improves the description of our previous model (Alijó et al., 2014a) by including electrostatic correlations under ac voltages; also allowing the exploration of new modifications, representing a possible way for future advances on the description of mobile charge dynamics via Poisson-Nernst-Planck approach.

CONCLUSIONS
Mobile charge dynamics under ac voltage blocking electrodes is treated in a Poisson-Nernst-Planck framework.We have proposed a new modified PNP model to account for ion size correlation (as well as other non-electrostatic) effects directly on chemical potential, in a thermodynamically consistent way.Electrostatic correlations are also introduced by means of a modified fourth-order Poisson equation.The model describes ion concentration profiles in accordance with Monte Carlo simulation results for 1:1 and 2:1 electrolytes, capturing the essential physics of overscreening close to charged surfaces.Here we propose a universal electrostatic correlation length 0.5 c   .Results for ac voltage mobile charge dynamics show that ion-ion correlations play an important role in order to properly model transient behavior of charging and discharging processes.Finally, we shall emphasize that the novel methodology proposed here turns possible to include new contributions on PNP-like approaches, paving the way to better understanding mobile charge dynamics analysis in electrokinetic and colloid transport phenomena.

Figure 1 -
Figure 1 -Sketches of the ion dynamics in 1D model problem.Electrolyte solutions confined between parallel-plate blocking electrodes separated by L under the influence of: (a) dc voltage; and (b) ac voltage harmonic potential.
resulting from nondimensionalization.The chemical potential i  is given by:

Figure 2 -
Figure 2 -Comparison of the mPNP model with the electrostatic correlations (blue solid line), and neglecting electrostatic correlations (black dashed line), to Monte Carlo simulations of Boda et al. (2002) (spheres) for (a) 1:1 electrolyte; and (b) 2:1 electrolyte.We shall emphasize that the mPNP model proposed here, with a characteristic correlation length of 0.5 c Fig. 3.The continuous lines represent the external ac voltage   s y  , whereas the surfaces are the system response after numerical solution of mPNP model.Concentration profile in Fig.3dshows an oscillatory behavior near the electrodes, which is a strong evidence of overscreening.

Figure 4 -
Figure 4 -mPNP model numerical solutions for a 2:1 electrolyte for: (a) and (b) electric potentials; (c) and (d) cation concentrations.Surface profiles presented in (a) and (c) are for solutions neglecting electrostatic correlations ( 0 c  ); whereas surface profiles presented in (b) and (d) are for solutions considering electrostatic correlations ( 0.5 c