Multiscale response of ionic systems to a spatially varying electric field

In this paper the response of ionic systems subjected to a spatially varying electric field is studied. Following the Nernst-Planck equation, two forces driving the mass flux are present, namely, the concentration gradient and the electric potential gradient. The mass flux due to the concentration gradient is modelled through Fick's law, and a new constitutive relation for the mass flux due to the potential gradient is proposed. In the regime of low screening the response function due to the potential gradient is closely related to the ionic conductivity. In the large screening regime, on the other hand, the response function is governed by the charge-charge structure. Molecular dynamics simulations are conducted and the two wave vector dependent response functions are evaluated for models of a molten salt and an ionic liquid. In the low screening regime the response functions show same wave vector dependency, indicating that it is the same underlying physical processes that govern the response. In the screening regime the wave vector dependency is very different and, thus, the overall response is determined by different processes. This is in agreement with the observed failure of the Nernst-Einstein relation.


I. INTRODUCTION
Ionic liquids, molten salts and ionic solutions show a response to the application of an external electric field. In the case of small field amplitudes the response is typically modelled by linear constitutive relations, for example, the response manifested by a charge current is related to the local field by Ohm's law and the charge density to the external field through the charge-charge correlation function [1,2]. The response is characterized by different response functions (or transport coefficients) like the electric conductivity, electric permittivity and the charge-charge correlation. One can often find relations between the different response functions [2]; at least in some limits. A more famous one is the Nernst-Einstein equation that relates the self diffusion coefficient to the electric conductivity [3], i.e., single particle flux to the charge current. This is a quite surprising relation as the particle flux is a single particle phenomenon whereas the charge current is a collective phenomenon. The Nernst-Einstein equation is then only valid when the ion cross-correlations can be neglected [1,4].
One example where this assumption is not valid is where the flux of ion-pairs contributes to the mass flux, but not to the charge current as the charges cancels [2]. The deviation can be determined from simulations or experiments and is often quantified by a deviation parameter [2], which, interestingly, Harris et al. [5,6] have expressed in terms of the velocity cross-correlation functions. Importantly, the failure of the Nernst-Einstein equation means that the particle flux due to the electric field cannot be modelled through Ohm's law directly.
Rather than approaching this problem through the deviation parameter it is appealing to take one step back and propose a linear constitutive relation that involves a new response function relating the mass flux to the external field directly. This is done in this paper.
The system's response is dependent on the wavelength of the external field, and this can be modelled through wavevector dependent response functions [2,7]. Investigating the wavevector dependence is relevant as the response can vary as function of length scale [8].
Also, this provide a way to probe a characteristic correlation lengths for a given system [8,9]; if the characteristic length scales are different for the different response functions this indicates that different physical underlying mechanisms are responsible for the system response. This is also addressed here.
The paper is organized as follows: In the next section the theory for the response of an ionic system subjected to a static sinusoidal external field is presented. In Sect. III molecular dynamics simulation results are presented and discussed, and, finally, in the last section conclusions from the work are drawn.

II. THEORY
We consider an ionic system composed of one cation and one anion specie. The ions are rigid meaning that any higher order induced effects and electron transfer mechanisms are ignored. The charges are ±q, respectively. Let i indicate either a cation or an anion, i.e., i = + or −, then the number density n i follows the balance equation [10] where n i c i is the diffusive flux and n i u the advective flux. The production term σ i accounts for additional forces that generate a local change in n i ; this includes application of an external electric field. The terms on the right-hand side of Eq. (1) can be expressed as the divergence of fluxes such that if one writes the production term as σ i = −∇ · j e i and n i c i = j d i we have for zero advection The system is kept away from equilibrium by application of a static spatially varying external electric field. The field points and varies in the direction parallel to the system z-direction, i.e, the non-zero z-component of the external field reads where k n = 2πn/L z is the wavevector, n = 1, 2, . . ., and L z is the system length in the zdirection. m is either 0 or 1. The experimental realization of this field is not straightforward.
Here it is considered as we are interested in the wavevector dependent response and as such this resembles the sinusoidal transverse and longitudinal force field methods (STF and SLF), see for example Refs. 11-13. The corresponding electric potential is using φ ext (0) = 0. Note, m = 0 corresponds to a wavevector independent field amplitude and m = 1 to wavevector independent potential amplitude.
It is in place to discuss the Maxwell equations. First, the induced/screening field is E = E(z) and according to Gauss' law dE/dz = ρ q /ǫ 0 , where ρ q is the charge density given by the induced ionic density, ǫ 0 is the electric permittivity of free space. From the Maxwell-Faraday equation ∇ × E = −Ḃ = 0, that is, the field due to the screeing does not result in any change in the magnetic field B. Then Gauss' law for the magnetic field is fulfilled, ∇ · B = 0. Furthermore, since ∇ × B = 0 andĖ = 0 there are no net charge current (Ampere's circuital law). The system is therefore in a steady state.
To proceed one needs to relate the fluxes with the corresponding forces [10]. For sufficiently small force amplitude this is done through the generalized linear response theory.
Consider the mass flux in the z-direction j i to depend on N forces X n , n = 1, 2, . . . N then we have in the homogeneous situation where χ ′ n is the response function relating the flux j i to the force X n . Since the system is in a steady state we can safely ignore time memory effects and, furthermore, assuming isotropy the response functions can then be written as The final expression is due to the steady state conditon . This generalized response formalism can be applied to the present situation. The flux is proposed to be given by the two terms where D i is the diffusion response function (or diffusion coefficient) and χ i is the response function that relates the mass flux to the external field. The first relation is simply a generalized version of Fick's law, but the second relation is not a generalization of Ohm's law as χ i relates the mass flux directly to external electric potential. Also, the force in Ohm's law is given by the local electric potential, i.e., the sum of the screening potential and the external potential. The χ i -response function can be interpreted as the system response to an external field excluding the effects from diffusion. Note that Eq. (7) is a generalized form of the Nernst-Planck equation [14].
In the steady state j d i + j e i = 0, and one has In Fourier space by the convolution theorem for wavevector k = (0, 0, k n ) this reads or Equation (10)  respectively. From this result one can also find the Fourier coefficients for the charge density, ρ q . First, it is observed that due to symmetry the number density follows a sine series, i.e., The Fourier components of the charge density is then For small field strengths and negligible screening only the fundamental mode k n = 2πn/L is excited and we have that In the following, focus is on the case where Eq. (13) is true and where the two ionic species, + and −, have same transport properties χ i = χ, and D i = D. Then Eq. (12) reduces to From linear response theory [2] the Fourier components for the charge density is related to the charge-charge correlation function (or charge-charge structure) S ZZ by where n = n + + n − . We then have an expression for χ in terms of the diffusion coefficient and the charge-charge structure The charge-charge structure is a collective property, and from Eq. (16) one can see that χ relates this collective property to the single particle property governed by the diffusion coefficient.
It is worth noting that in the Debye-Hückel regime, k B T ≫ qφ, the charge-charge correlations are negligible, i.e., S ZZ (k) = 1. This corresponds to the limit of zero screening and a relative permittivity of unity. Equation (16) then reads which is equivalent to the Nernst-Einstein equation [3] and χ can in this limit be interpreted as the ionic electric conductivity. The charge density Fourier components are in this limit , they only dependent on amplitude of the external field. For systems This means that the wavevector dependent response in the presence of an external electric field is dominated by the screening effects.

A. Simulation details
The response is evaluated for two simple models: (i) one model for molten salt proposed by Hansen and McDonald [1] and (ii) one modified model for ionic liquids used by Chapela et al. [15]. For the molten salt the ions are simple spherical particles with same mass and point charges ±q. The van der Waals interaction is the inverse power law function V (r) = ǫ(σ/r) 9 , where r is the distance between two ions, ǫ and σ define the energy and length scale, respectively. The Coulomb interactions are calculated through the shifted force method [16,17], F(r) = q i q j (1/r 2 − 1/r 2 c )r/r, for r ≤ r c . Here r is the vector of separation with magnitude r, and r c is the cut-off radius set to r c = 3σ; this cut-off distance is also used for the van der Waals interactions. The positions of the particles are integrated forward in time with the leap-frog algorithm [18] and the temperature is controlled using a Nosé-Hoover thermostat [19,20]. In all simulations the total ion number density is n = 0.368σ −3 ; the number of ions are 1000, giving 500 ion-pairs. Two different temperatures are simulated, T = 0.0177ǫ/k B and 1.0177ǫ/k B , the former being a realistic temperature for the model. To simulate the Debye-Hückel regime k B T ≫ qφ the ion-ion Coulomb interactions are removed whilst keeping the temperature fixed at T = 1.0177ǫ/k B ; this system is symbolized using T ∞ . Alternatively, one can perform simulations at very high temperatures, but this will result in numerical instabilities. In the following all quantities are given in units of σ, ǫ, q, and mass m, and as it is common practise these the units are not written explicitly.
For the simple molten salt system the shifted force method can be tested against the direct Ewald summation method [21]. From equilibrium simulations it was found from the structure that the Ewald method converges satisfactory using 124 replica systems and that it agrees with the data from the shifted force method, see also Ref. 17. For the nonequilibrium situation at T = 1.0177 the Ewald and shifted force methods yield same results for all wavevectors tested 0 < k < 2.2.
The modified ionic liquid model is composed of cations with a spherical point charge particle (head group) and two spherical non-charged tail particles. The particles in the cation are linearly connected using a simple spring force F = −k(r − 1)r/r, where k = 100 is the spring constant. Anions are simple spherical point charge particles [15]. Rather than a hard-sphere type potential in the original model, the van der Waals interactions are here given through the Weeks-Chandler-Andersen potential [22] V (r) = 4((1/r) 12 where the cut-off is set at r c = 2 1/6 .The Coulomb interaction is given by the Yukawa potential V (r) = q 2 e −r/λ D /r, with λ D = 1/2 corresponding to a relative small Debye screening length and the reduced charge is q = 4. The cut-off distance for the Yukawa potential is set to r c = 2.5. The state point is (n, T ) = (1, 1) and the simulation method is the same as for the molten salt simulations. This choice of parameters gives, qualitatively, the fluid structure observed in different ionic liquid [23,24]. The number of particles are 864, that is, 216 ion pairs.
Simulations of the non-equilibrium system is also performed. Here an additional force from the external field, Eq. (3), is added to the total force experienced by the ions F ext i = q i E ext k, where k is the unit vector parallel to the z-axis.

B. Results: Molten salt
The wavevector dependent diffusivity can be obtained as follows. The Gaussian approximation [2,25] relates the diffusion coefficient to the incoherent intermediate scattering function (or the self-part of the density-density correlations), so in the diffusive regime, i.e., for large t, this is here generalized to The Fourier-Laplace transformation is which gives an expression for the wavevector dependent diffusivity in the limit of zero frequency Microscopically the intermediate scattering function is defined from the ensemble average where N is the number of ions and is thus a single particle property. In Fig. 1  Also, shown as punctured lines f (k, t) = e − 1 6 ∆r 2 k 2 t , where ∆r 2 is the particle mean square displacement. It is seen that the Gaussian approximation holds surprisingly well for this model validating Eq. (19). The data are Fourier-Laplace transformed and the Gaussian diffusion kernel is found from Eq. (21); the results are plotted in Fig. 1 (b). The function is fitted to data where the zero wavevector diffusion coefficient, D 0 , is found from the mean square displacement ∆r 2 = 2D 0 t for t → ∞. It is observed that the normalized kernel is  Next the charge-charge structure is evaluated. This is defined as [2] S ZZ (k) = 1 N ρ q (k, 0)ρ q (−k, 0) , where ρ q (k, 0) = i q i e −ik·r i . Note, this is a collective property. For non-zero wavevectors S zz (k) can also be calculated from the radial distribution functions, see e.g. Ref. 1, where ∆g(r) is the difference between the cation-cation and cation-anion radial distribution functions, ∆g(r) = g ++ (r) − g +− (r). The charge-charge structure is plotted in Fig. 2  the response function is monotonically decaying with respect to wavevector. This behavior is typically observed for the diffusion and viscosity kernels [8]. For non-zero screening the response features a maximum depending on temperature; the characteristic wave length l = 2π/k max where k max is the wavevector corresponding to maximum in χ, is approximately l = 2.5 for T = 1.0177 and l = 1.6 for T = 0.0177. This means that application of an external field will result in a relatively small flux, j e i , on large length scales and a maximum for wavelength of roughly 2 atomic diameters. For T = 0.0177, we have that lim k→0 χ(k) = 0 which means that the charge density is zero at these length scales; this is in agreement with perfect screening. Another important point is that lim k→∞ S zz = 1, and if lim k→∞ D = 0 as indicated in Fig. 1 we have that lim k→∞ χ = 0 according to Eq. (16). show good collapse, that is, there exists a master curve response function. This identical wavevector dependence indicates that the response functions are governed by the same underlying process. Specifically, it is here conjectured that the χ-response is given by the diffusion processes in the system, i.e., cross correlation effects can be ignored. From Fig. 3 (a) one can immediately see that this collapse is not found for the T = 1.0177 and T = 0.0177 cases, hence, different processes are involved.
The theory is compared with the non-equilibrium simulations. Figure 4 (a) shows the charge density profile, ρ q , for two wavevectors k = 2π/L and k = 8π/L at T = 0.0177. The system length is L = 13.955 and m = 1, hence, the potential field amplitude is constant. It is observed that the charge density amplitude is larger for smaller wavelengths as expected.
For k > 12π/L a simple spectral analysis shows that higher order modes are excited compromissing Eq. (13) and only results for k < 12π/L is shown. Figure 4 (b) compares the amplitude for all three temperatures with the predictions from the theory, Eq. (14). The agreement is excellent. Of course, this comparison is equivalent to test the linear response, Eq. (15). The case of m = 0 is also shown, however, the agreement is less satisfactory for low wavevectors, which is due to the diverging amplitude in the limit of zero wavevector causing a non-linear response and failure of the constitutive relation, Eq. (7).
In Fig. 5 (a) the diffusion kernels are is shown for the ion liquid model. These are evaluated as explained in Sect. III B. One sees that within statistical uncertainty the diffusion kernel is wavevector independent, at least up to k = 5. Beyond this wavevector value the statistical error increases dramatically and the results are non-conclusive. The charge-charge structure, Fig. 5 (b), is calculated from the direct definition Eq. (24). It features relatively strong structure, that is, the system is in the screening regime. We can therefore expect the response function χ to resemble low temperature molten salt response function.
For the ionic liquid Eqs. To investigate if the two kernels can be mapped onto the same master curve, the results from the cation kernel is normalized with respect the maximum. The normalized result is shown in Fig. 5 (c) as squares connected with a punctured line. To a reasonable agreement the two kernels do follow a master curve which indicates that the underlying mechanisms responsible for the response are the same. This contrasts the wavevector independent diffusion kernel, that is, the system response seen in the mass flux due to the density gradient.
Therefore, the physical mechanisms for the two fluxes j d i and j e i are fundamentally different; at least in the screening regime.

IV. CONCLUSION
In this paper the mass flux of an ionic system due to a spatially varying electric field is studied. Following the Nernst-Planck equation two forces are present in this system: (i) the concentration gradient and (ii) the gradient of the electric potential. The two response functions (or kernels) that account for the system response to these forces are the diffusion- and χ-response functions; the χ-response function relates the mass flux with the external electric field excluding the contribution from the concentration gradient (here modelled through the self-diffusion). Note, this differs from the charge-charge response function, S zz , which relates the charge density to the electric field including all underlying processes, and the ionic conductivity that relates the charge current to the local field.In the limit of zero screening the χ-response function is directly related to the conductivity, on the other hand, in the large screening regime the response function is related to the charge-charge structure.
The spatial correlations in the system are manifested in the wavevector dependence of the kernels. The molecular dynamics simulation data show the diffusion and χ-kernels feature very different wavevector dependence in the screening regime. Interestingly, in the screening regime both the molten salt and ionic liquid feature a wavevector independent diffusion kernel and the response to the external field is dominated by the charge-charge structure. This latter quantity is a collective property of the system. In the non-screening regime, on the other hand, the response to the external field is closely related to the ionic conductivity and in this regime the Nernst-Einstein relation holds to a good approximation, i.e., cross-correlation effects are negligible.