Transient fluctuation-induced forces in driven electrolytes after an electric field quench

Understanding how electrolyte solutions behave out of thermal equilibrium is a long-standing endeavor in many areas of chemistry and biology. Although mean-field theories are widely used to model the dynamics of electrolytes, it is also important to characterize the effects of fluctuations in these systems. In a previous work, we showed that the dynamics of the ions in a strong electrolyte that is driven by an external electric field can generate long-ranged correlations manifestly different from the equilibrium screened correlations; in the nonequilibrium steady state, these correlations give rise to a novel long-range fluctuation-induced force (FIF). Here, we extend these results by considering the dynamics of the strong electrolyte after it is quenched from thermal equilibrium upon the application of a constant electric field. We show that the asymptotic long-distance limit of both charge and density correlations is generally diffusive in time. These correlations give rise to long-ranged FIFs acting on the neutral confining plates with long-time regimes that are governed by power-law temporal decays toward the steady-state value of the force amplitude. These findings show that nonequilibrium fluctuations have nontrivial implications on the dynamics of objects immersed in a driven electrolyte, and they could be useful for exploring new ways of controlling long-distance forces in charged solutions.

Understanding how electrolyte solutions behave out of thermal equilibrium is a long-standing endeavor in many areas of chemistry and biology. Although mean-field theories are widely used to model the dynamics of electrolytes, it is also important to characterize the effects of fluctuations in these systems. In a previous work, we showed that the dynamics of the ions in a strong electrolyte that is driven by an external electric field can generate long-ranged correlations manifestly different from the equilibrium screened correlations; in the nonequilibrium steady state, these correlations give rise to a novel long-range fluctuation-induced force (FIF). Here, we extend these results by considering the dynamics of the strong electrolyte after it is quenched from thermal equilibrium upon the application of a constant electric field. We show that the asymptotic long-distance limit of both charge and density correlations is generally diffusive in time. These correlations give rise to long-ranged FIFs acting on the neutral confining plates with long-time regimes that are governed by power-law temporal decays toward the steady-state value of the force amplitude. These findings show that nonequilibrium fluctuations have nontrivial implications on the dynamics of objects immersed in a driven electrolyte, and they could be useful for exploring new ways of controlling long-distance forces in charged solutions.

I. INTRODUCTION
For more than a century, extensive research efforts have focused on understanding the properties of electrolytes and charged solutions [1][2][3]. Conventional mean-field theories have been remarkably effective in such endeavors and, in addition, the electrostatic correlation effects have been characterized to a great extent [4]. Similar studies of correlation effects and fluctuation phenomena in nonequilibrium electrolytes are, however, relatively scarce. The present manuscript focuses on one such aspect, namely fluctuation-induced forces (FIFs) [5] in electrolyte solutions driven out of equilibrium, by considering the transient behavior of the FIF after a quench by an electric field.
Mean-field descriptions of electrolytes are obtained by combining those equations that govern the electrostatic interactions among the charged particles (e.g., the Poisson equation) with the statistical weights that describe the distribution of the charges. In thermal equilibrium, where the statistical weights are given by Boltzmann factors, this combination leads to the so-called nonlinear Poisson-Boltzmann (PB) equation [1], while out of equilibrium one needs to consider more general distributions [6] which are governed by Fokker-Planck equations [7,8]. Even though the resulting nonlinear models in principle contain the full mean-field information about the system, it is often difficult to gain insights into the underlying physical mechanisms from them; as such, linear models have been of great importance in the conceptual developments toward understanding electrolytes and charged fluids. The celebrated Debye-Hückel (DH) theory [9], for instance, is obtained from linearizing the PB equation, and it shows how correlations in an electrolyte become short-ranged as a result of the screening effects of counterions (i.e., opposite charges). The simple picture provided by the DH theory serves as a general starting point to understand and construct models for a wide range of charged systems [1].
The mean-field theories, however, do not generally take into account the statistical correlations between the ions in the system. It is well-known that the correlations in electrolytes can give rise to a number of important phenomena, ranging from phase transitions in two-dimensional systems [10], to charge renormalization in colloids [11], and counterion condensation and correlation-induced interactions which are relevant biological processes [12][13][14][15][16]. In these cases, it can still be useful to combine the phenomenology of the correlation effects, which can lead to the breakdown of the mean-field assumptions, with the linearized DH equations to get mathematically tractable models [4]. Although some related phenomena have also been considered in dynamical settings [17][18][19], a general understanding of nonequilibrium electrostatic correlations is still far from complete and poses a number of outstanding open questions [7]. This is in part due to the inherent mathematical difficulty of characterizing the distribution functions of interacting ions with Brownian motions out of thermal equilibrium [6,20]. Langevin formulations of the electrolyte dynamics provide a straightforward way of taking into account the stochasticity in the motion of the ions and can be used to obtain the required correlation functions, e.g., to compute the power spectrum of nanopores or the conductivity of strong electrolytes [21][22][23]. Inspired by some of the recent investigations of electrolytes and ionic liquids [24][25][26][27][28][29][30][31], we use a similar mathematical framework to study how fluctuations in a strong electrolyte in the presence of an external electric field may modify the force that is exerted on boundaries that confine the electrolyte. This direction would also be relevant for designing more efficient and environmental-friendly batteries and electrochemical capacitors [32][33][34] as electrolytes and ionic liquids are essential to the underlying conduction processes and energy storage.
Fluctuation-induced forces (FIFs) comprise a remarkable aspect of fluctuation effects and arise when external objects disturb a correlated medium by imposing specific boundary conditions on its fluctuation modes [5]. If the correlations in the medium are long-ranged, they then give rise to FIFs that can persist between objects at large distances and exhibit a number of universal features [35][36][37][38]. Fluctuation forces have many applications in nanosciences [39] and colloidal systems [40] and they are also prevalent in nonequilibrium and driven settings [41][42][43][44][45][46][47][48][49][50]. For an electrolyte in thermal equilibrium, the spatial range of the FIF is limited by the Debye screening length, which is often of the order of a few nanometers [1,51,52]. Since nonlocal correlations generically emerge from conserved dynamics [53][54][55], it may not be surprising that nonequilibrium FIFs in driven electrolytes become long-ranged [56]. Since the underlying cause for this FIF, namely anisotropy in the conserved dynamics, is rather generic, similar nonequilibrium forces might also be relevant to modeling electrokinetic and biological systems where there are many instances of ionic currents in confined spaces [8,57].
In a previous work [56], we used the Dean-Kawasaki formalism [58][59][60] along with scaling arguments to investigate the steady-state fluctuations in a simple strong electrolyte that is driven by a constant external electric field. It was shown that the anisotropy caused by the electric field leads to the so-called generically scale invariant dynamics [61], and the correlation functions in the steady-state become long-ranged and take power-law forms, despite the screening effects. These correlations then lead to nonequilibrium FIFs on the confining boundaries, which in flat geometry and at the steady state, depends on the plate separation H as H −d , where d is the spatial dimensionality. The FIF in this case has unique non-monotonic changes as the strength of the applied electric field is increased. In the present manuscript, we look into the transient features of the FIF after the electric is switched on at t = 0. To elucidate the fundamental concepts, we keep our focus on the case of a strong binary electrolyte which is initially at thermal equilibrium and is confined between two neutral parallel plates. The electrolyte is acted upon at time t = 0 by a DC electric field parallel with the boundaries (see Fig. 1). We investigate the temporal variations of the FIF and show that the FIF exhibits a diffusive behavior in time where its difference with the steady-state value decays as t −d/2 in d spatial dimensions (with a possible crossover in the value of the exponent, see Section IV). These results cast some light on the emergent behavior in electrolytes out of thermal equilibrium and demonstrate that the dynamical and nonequilibrium effects in electrolytes can be manifestly different from their screened counterparts in equilibrium.
The manuscript is organized as follows: in Section II, we start from the Langevin equations that govern the dynamics of the ions and then derive the stochastic density equations based on the Dean-Kawasaki approach. This is then accompanied by a linearization and quasi-stationary approximation schemes. We also discuss the applicability of these approximations to settings where the electric field is a slowly varying function of time and, in addition, to solutions with different mobilities of charge species (which could be relevant for maintaining an unscreened electric field in the bulk). In Section III, we use the linearized density equations to derive the charge and density correlation functions; first, this is done for a bulk solution and then for an electrolyte that is confined between two flat boundaries. In Section IV, we compute the stress exerted by the electrolyte on the confining plates. To this end, we use Maxwell stress and make use of the obtained correlation functions to obtain analytical expressions for the normal stress on the boundaries. The steady-state and transient parts of the stress amplitude are investigated using asymptotic expansions and numerical evaluations. Section V contains the summary and final concluding remarks. There are also four appendices containing the scaling analysis of the nonlinearities (Appendix A), the calculation of the two-point correlation functions of the linearized dynamics at equal times (Appendix B) and at different times (Appendix C) without using the quasi-stationary approximation, and a simplification of the transient part of the FIF amplitude (Appendix D).

II. FORMALISM
To analyze the dynamics of the electrolyte in the presence of the external electric field, in Section II A we start from the single particle Langevin equations that describe the motion of the individual ions, and then we recast the set of Langevin equations in terms of the time evolution of their instantaneous density using the Dean-Kawasaki approach [58][59][60]. The density equations are then linearized in Section II B whose applicability on relaxing the assumptions of static electric field and equal ionic mobilities are discussed in Sections II C and II D.

A. Stochastic density equations
The electrolyte we study here is composed of a symmetric collection of charged particles with charges ±Q that have equal mobilities µ + = µ − = µ. These particles undergo Brownian motion due to the microscopic collisions with the implicit solvent molecules (we will neglect the hydrodynamic effects of the solvent, assuming a finite (screened) range of hydrodynamic interactions [62][63][64][65][66]). In d spatial dimensions, the overdamped trajectory r ± a (t) of a cation or anion, labeled by a (a = 1, 2, . . . , N ) is governed by the Langevin equatioṅ where E = Eê x is the external electric field along the x axis and η ± a represent independent Gaussian white noises with correlations η ± ai (t)η ± bj (t ) = δ ab δ ij δ(t − t ) and zero means. The fluctuation-dissipation relation connects the noise strength D to the mobility µ through the Einstein relation µ = βD, where β = 1/(k B T ) represents the inverse temperature of the system. The electrostatic potential field φ(r, t) is created by the ions in the system, and in the Gaussian units it satisfies the Poisson equation where S d = 2π d/2 /Γ(d/2) is the surface area of d-dimensional unit sphere, in is the permittivity of the electrolyte, and the instantaneous density of charge species is defined as Using the stochastic equations of Dean and Kawasaki, one then obtains the exact stochastic dynamics of C ± as continuity equations, i.e. ∂ t C ± + ∇ · J ± = 0, where the stochastic currents are given by Here, ζ ± (r, t) are uncorrelated Gaussian noise fields characterized by zero averages and ζ ± i (r, t)ζ ± j (r , t ) = δ ij δ d (r− r )δ(t − t ).
We next shift our focus to the dynamics in terms of the number density C = C + + C − and the charge density ρ = C + − C − , for which the density equations can be recast as where the density and charge currents are defined as J c = J + + J − and J ρ = J + − J − and they explicitly read Note that the noise fields √ 2DCζ c,ρ are obtained from the addition and subtraction of the Gaussian noise fields √ 2DC ± ζ ± ; consequently, these uncorrelated fields also have zero averages with correlations given by Substituting the density and charge currents J c and J ρ into the corresponding continuity equations, we find for the dynamics where now the electric potential satisfies −∇ 2 φ = (S d Q/ in )ρ. Note that the electric field introduces a source term ∝ ∂ x ρ to the dynamics of C and vice versa; this is the origin of the long-range behavior we will obtain later, since it leads to the modification of the effective density diffusion coefficient parallel to the electric field.

B. Linearized density equations
As Equations (7) and (8) contain nonlinear terms and multiplicative stochasticity, they are in general difficult to analyze. It is thus customary to focus on the fluctuations of the densities fields about a uniform background baseline C 0 (i.e, C ± = C 0 + δC ± with |δC ± | C 0 ), and to only keep terms in leading order of the fluctuations [21-23, 65, 67]. A scaling analysis of the nonlinearities also shows that they become negligible in the macroscopic limit (see Appendix A). We define the density fluctuation field c(r, t) and the charge fluctuation field (in units of Q) ρ(r, t) as Expanding Eqs. (7) and (8) and keeping the linear terms, one arrives at the following linear stochastic equations for the dynamics of c and ρ: In these equations, we have defined the Debye length κ −1 through and the linearized noise correlations now read η ρ (r, The linearized density and charge equations (10) and (11) will be used in Section III to compute the correlation functions. This linear set of equations allows for the full characterization of the correlation functions in the Fourier space (see Appendices B and C). However, since we are interested in the long-distance asymptotic limit, a convenient simplification can already be made by noting that for time and length scales beyond those set by the Debye screening processes, the charge fluctuation has a quasi-stationary solution given by which is obtained by neglecting the temporal and spatial derivatives of the charge density [68]. Since the typical Debye screening length is of the order of κ −1 ∼ 1 − 10 nm [1], such approximation is justified for studying the dynamics of the electrolyte beyond the screening scale. Accordingly, the corresponding FIFs that will be calculated below can in principle be realized in settings where the boundary separations are larger than the screening scale of the electrolyte which, for instance, may be the case for wet ion channels such as mechanosensitive channels [69], synthetic nanopores [70], and force measurement setups with large inter-plate separations [29,31]. Equation (13) shows that the charge density which persists beyond the relaxation time ∼ 1/(Dκ 2 ) is proportional to the gradient of the number density along the direction of the electric field. This is in principle due to the bias introduced to the dynamics of the charged particles by the external field, which competes with the electrostatic forces that tend to relax any excess charge in the bulk of the electrolyte. Upon substituting Eq. (13) back into the number density dynamics Eq. (10), we arrive at an anisotropic diffusion equation that reads where the dimensionless electric field E is defined as Equation (14) is the central result of this section and one of the main points of this work as it also underlies the long-range FIF that will be calculated later. Note that the Einstein relation µ = βD that holds between the mobility and the noise strength in the single particle Langevin dynamics (1), is no longer valid at the macroscopic level of Eq. (14). This introduces a mismatch between the noisy fluctuations and the dissipative forces in the conserved dynamics of the density fluctuations, as a result of which the dynamics of c represents a realization of generic scale invariance [53,54,61].
Finally, it is worth noting that E can be regarded as the ratio of the energy density of the electric field ( inE 2 2S d ) to the thermal energy density (C 0 k B T ) of the electrolyte, which is also the ratio between the typical values of the corresponding stress components. Alternatively, E could be interpreted as the average distance traversed along the electric field by a charged particle during the charge relaxation process (velocity × time ∼ (µQE) × (Dκ 2 ) −1 ), divided by the Debye length (κ −1 ) which determines the spatial extent of the counterion cloud around a cation or anion. As such, E qualitatively encodes the extent of the deformation of the counterion atmospheres from their equilibrium (symmetric) forms as a result of the applied external field [6].

C. Slowly varying electric fields
In the preceding analysis, we have assumed the electric field is DC, i.e., it is constant in time. However, maintaining a constant electric field in the bulk of an electrolyte is difficult in experimental setups, as ions will quickly accumulate on the electrodes with opposite charges and screen their electric field to a short-ranged residual field acting only on the boundary layers. This naturally raises the question of whether the analysis that led to derivation of Eq. (14) can be extended to cases where the driving field is time-dependent? To answer this, we note that the main simplifying approximation that allowed us to obtain the anisotropic diffusion process, Eq. (14), is the quasi-stationary relation between charge and density fluctuations Eq. (13). In physical terms, this approximation assumes the charge relaxation processes take place over time-scales that are much faster when compared to other dynamical scales. Such an approximation is therefore applicable to time-dependent electric fields as far as the variations in the field strength over the charge relaxation time is negligible; in that case, the ionic atmosphere surrounding each ion will have enough time to rearrange itself to the moving position of the central ion before the applied field changes considerably [20]. As the frequency set by the relaxation processes is often of order of few nano-seconds, the condition of slowly changing electric field is, in fact, often met in setups with oscillatory electric fields. In such cases, the anisotropic diffusion of the density fluctuation field c(r, t) is therefore given by with E(t) = βQE(t)/κ. This equations shows that with a slowly varying electric field, the diffusion coefficient for density fluctuations along the direction of the field becomes a time-dependent factor. It is finally worth noting that as far as the time-scale for the variations of this diffusion coefficient is slower with respect to the charge relaxation processes, it can, in principle, be faster than other (diffusive) time-scales in the system.

D. Unequal mobilities of cations and anions
Another simplifying assumption in deriving Eq. (14) was the equal mobilities and diffusivities for the cations and anions in the electrolyte solution. However, often the cations and anions in an electrolyte have different ionic radii and as such their mobilities and diffusivities are not the same [1]. A crucial observation is that in oscillatory electric fields, such differences in the ionic mobilities lead to a mismatch in the spatial range of the motion of the ions within the solution which, in turn, gives rise to an average electric field in the bulk of the solution [71]. In fact, this observation suggests a way of maintaining unscreened fields in the bulk of the electrolyte. Having this in mind, we now show that once an electric field is set in the bulk, a difference in the mobilities µ ± of the charge species would not change the macroscopic dynamics described by Eq. (14).
Let us consider the mobilities and diffusivities of the cations (+) and anions (−) to be given by and assume they are related through the Einstein relation µ ± = βD ± (note that δD and δµ may be positive or negative and are not necessarily small). Following the same steps outlined in Section II A, and using the definition of the number density C = C + + C − and the charge density ρ = C + − C − , we obtain the modified version of the full number density and charge currents as A similar linearization procedure as in Section II B yields the dynamics of the fluctuations in the number and charge densities according to It could be seen that the quasi-stationary approximation of Eq. (13) still holds for Eq. (21) since the added terms ∝ δD and ∝ δµ with higher number of derivatives are not relevant at macroscopic scales (this can be seen more rigorously through a scaling analysis similar to Appendix A). Therefore to obtain the dynamics of the density fluctuations c, we substitute Eq. (13) into the modified density dynamics given by Eq. (20); noting that δµ/δD = µ/D = β, we recover the same anisotropic diffusion equation (14) derived for the case with equal mobilities and diffusivities (with the only modification being that the diffusion coefficient should now be replaced by D = (D + + D − )/2). We conclude that as long as an electric field is set in the bulk of the electrolyte, the difference in the diffusivities of the cations and anions may only affect the microscopic dynamics of the particles, but it does not change the effective macroscopic equations governing the long-distance dynamics of the density and the charge fluctuations.

III. CORRELATION FUNCTIONS
We now turn to calculating the correlations of the density fluctuation field c(r, t) and the charge fluctuation field ρ(r, t) after a quench in the electric field. In Section III A, we first compute the correlation functions of the density fluctuations in the bulk, using the appropriate Fourier representation of Eq. (14), and then make use of Eq. (13) to obtain the bulk charge correlations. We then compute the correlations in the presence of neutral confining boundaries parallel with the electric field in Section III B. Our focus will be on the case that a constant (DC) electric field is switched on at t = 0, before which the electrolyte is assumed to have been at thermal equilibrium.

A. Bulk correlation functions
Let us first derive the bulk correlations in the absence of any confining boundaries. The translational symmetry makes it convenient to work with the spatial Fourier transform of Eq. (14) which reads We use the Fourier convention f (k) = d d r e −ik·r f (r) with r = (y, s) ∈ R d where s = (s 1 = x, s 2 , . . . , s d−1 ) ∈ R d−1 , and k = k s + k yêy with k s = (k s1 = k x , k s2 , . . . , k s d−1 ). The noise correlation in the Fourier space reads Using an integration factor, Eq. 22 can be solved as Note that for a time-dependent electric field, this expression should be modified to In other words, the expression for the density fluctuations when the electric field is time-dependent can be obtained from that of the static electric field by making the substitution . For instance, for a periodic drivin field E(t) = E 0 cos(Ωt), this substitution reads E 2 → (E 2 0 /2) (1 + sin(2Ωt)/(2Ωt)), which at long time scales simply reads E 2 → E 2 0 /2. Keeping in mind that this substitution is also be applicable to the results that will follow, from here on we will solely focus on the case of static electric fields.
Equipped with Eq. (23), it is now straightforward to obtain the density correlations as where the averaging is performed over the noise realizations, and we have assumed t ≥ t without loss of generality. The first line of Eq. (25) represents the decay of the initial conditions; the second line shows the diffusive propagation of the density fluctuations between two different times, as well as the establishment of the steady-state correlations; and the third line encodes both the steady-state and the transient correlations of the nonequilibrium fluctuations which vanish for E = 0. Hereafter, we restrict the calculation to electric field quenches applied to electrolyte solutions which are initially in thermal equilibrium at t = 0; performing an ensemble averaging over the thermal initial configurations (denoted by . . . IC ) implies A similar averaging on Eq. (25) for t = t gives the density correlation functions at equal times according to which also defines the distinct part of the bulk correlation function, c bulk (k, t) [72]. Note that the local part of the correlation function (∝ C 0 δ d (k + k )) represents the screened equilibrium correlations in the long distance limit (see Appendix B for the full expression without taking this limit). Transforming c (2) bulk (k, t) to the real space yields wherer is obtained from r by substituting x →x = x/ √ E 2 + 1. In Eq. (28), the first term on the r.h.s gives the nonequilibirum steady-state correlation function which vanishes without the external electric field, and in d spatial dimensions it decays as ∼ r −d with distance; moreover, this term is manifestly anisotropic with a dipolar character [61]. The second term on the r.h.s, on the other hand, represents the transient effects with their long-time decay governed by the power-law tail ∼ t −d/2 .
Next, we turn to the charge fluctuation correlation functions. Incorporating the quasi-stationary charge profile of Eq. (13), one obtain the bulk charge correlations upon taking derivatives of Eq. (27) as bulk (k,t) .
Note that we have not included the second term in the brackets in ρ (2) since it originates from the density selfcorrelations in Eq. (27); the exact expressions provided in Appendix B reveal that this term represents the asymmetry in the screened correlation functions caused by the external field and is another local contribution to the charge density (∝ ∂ 2 x δ d (r) in real-space). Such local terms will not contribute to the long distance behavior of the fluctuation-induced forces (Section IV), and therefore will be discarded from the subsequent calculations. The real-space form of the nonlocal charge correlations ρ (2) bulk can also be obtained from Eq. (28) by taking derivatives with respect to the x coordinate (i.e., along the electric field) in accordance with Eq. (13); it is therefore seen that such charge fluctuations are also long-range correlated in the electrolyte solution.
Equations (27) and (29) give both the transient as well as the steady-state correlation functions, within the length and time scales where the approximation Eq. (13) holds. In Appendices B and C, without making use of this approximation, we derive the density and charge two-point correlation at steady-state, for fluctuation field considered at equal times as well as at different times. One can go further and obtain the full correlations in the transient regimes by solving the stochastic dynamics of Eqs. (10) and (11), which in the matrix form read The formal solution of this matrix equation are given by which then allows for computing the full correlation functions. These expressions are rather cumbersome and will not be evaluated here. Nevertheless, it is clear that the eigenvalues of M control the approach of the correlation functions toward their steady-state form. These eigenvalues read The relative values of κ, E, and k x determine whether the eigenvalues have an imaginary part. Two different dynamical behaviors are inferred from these eigenvalues: for k x < κ/(2E), both eigenvalues remain real, and the full solution is the superposition of a (fast) decaying term with relaxation time 1/(Dκ 2 ), and a soft diffusive mode. For k x > κ/(2E), on the other hand, the eigenvalues also acquire an imaginary part; in this case, both contributions to the full solution are damped oscillatory with relaxation time 1/(Dκ 2 ). (A similar behavior also appears in two-point correlation functions evaluated at different times, while the equal-time correlations at steady state remain the same in both regimes, see Appendix C.) These results imply that Eq. (14) captures the diffusive dynamics of the electrolyte at length scales beyond 2E/κ, while for smaller length scales the dynamics is relaxational. This therefore introduces an additional scale for the dynamics of the solution should be taken into account together with the fact that the approximate charge profile given by Eq. (13)is already restricted to scales beyond κ −1 .

B. Correlation functions in flat confinement
We now turn to computing the density and charge correlation functions of the driven electrolyte in the presence of flat boundaries located at y = 0 and y = H (Fig. 1). The boundaries are assumed to be impenetrable (i.e., they are blocking electrodes), and therefore they impose no-flux Neumann boundary conditions on Eq. (14), namely ∂ y c| y=0 = ∂ y c| y=H = 0. One can construct the corresponding solutions that satisfy these boundary conditions by making use of the Neumann eigenmodes cos(p n y) with p n = nπ H through [73] c(r, t) = where s = (s 1 = x, s 2 , . . . , s d−1 ) ∈ R d−1 is the position along the boundary surfaces, k s = (k s1 = k x , k s2 , . . . , k s d−1 ) is the corresponding momentum vector (with k s = |k s |), and n indicates a summation where the n = 0 term takes an additional factor of 1/2. Performing the similar transformation on the noise term η c , the correlation functions between the density modes c n (k s , t) are then obtained, which after simplification read c n (k s , t)c n (k s , t ) = (2π) d−1 δ d−1 (k s + k s ) δ n,n (1 + δ n,0 ) This expression could be used in order to obtain the following partial Fourier transformation of the density correlation functions along the s coordinates: cos(p n y) cos(p n y ) Narrowing down to correlations of density fluctuations at equal times, similar to the case of bulk correlations we define c(y; k s ; t)c(y ; k s ; t) = (2π) d−1 δ d−2 (k s + k s ) 2C 0 δ(y − y ) + c (2) (y, y ; k s ; t) , and get cos(p n y) cos(p n y ) Here, the first line is an equilibrium-like diffusive correlation which vanishes at long times, and the second line represents the nonequilibrium transients that in the steady-state give rise to a long-ranged term. Upon taking derivatives of Eq. (35), we get for the charge correlation functions at equal times cos(p n y) cos(p n y ) Note that the first line hear is due to the short-ranged screened correlations (since ∞ n=0 (2/H) cos(p n y) cos(p n y ) = δ(y − y )), and the second line is the long-ranged nonequilibrium contribution. It is worth mentioning that the correlation functions for the two half-spaces out of the confined space (i.e., for y < 0 and y > H) are obtained from Eqs. (35) and (36) by making the substitution 1 H n g(p n ) → ∞ −∞ dp 2π g(p) where p n = nπ H and g(p n ) stand for the summand in Eqs. (35) and Eq. (36). These bulk correlations will be used to compute the stress exerted from the electrolyte outside the confinement, which are then subtracted from the internal stress to obtain the net total stress on the boundaries.

IV. STRESS TENSOR
In this section, we turn to calculating the stress exerted by the driven electrolyte on confining parallel boundaries. Since the system under consideration is out of thermal equilibrium due to the driving electric field, the stress or pressure cannot be calculated from thermodynamics. Based on the mechanical definition of stress, one can instead obtain stress formulas that also work in nonequilibrium conditions [74,75]; such procedure for electric forces yields the well-known (electrostatic) Maxwell stress [76,77] which also circumvents the ambiguity often faced in deriving the stress tensor from body forces (due to the freedom in inverting the divergence operator). In Section IV A, we derive a simplified formula for the noise-averaged Maxwell stress exerted by a general charge distribution confined between two parallel plates, assuming the charge distribution is invariant under translations along the boundaries. In Section IV B we implement the charge correlations given by Eq. (36) into this formulation to obtain the FIF, whose steady-state and transient parts are analyzed in Sections IV C and IV D.

A. General expression for Maxwell stress in plane parallel geometry
For an electrostatic potential field φ, the noise-averaged Maxwell stress tensor in d spatial dimensions reads where the electric potential satisfies the Poisson equation −∇ 2 φ = S d ρQ/ in . With free boundary conditions, the solution to this Poisson equation is given by φ free (r, t) = Q in(d−2) d d r ρ(r ,t) |r−r | d−2 (for d = 2, the potential is given by a logarithmic Coulomb form, but the electric field relation and the Maxwell stress formula remain unchanged.) To proceed with Eq. (37), the correlation function of the electric potential φ(r, t) in confinement are needed. We thus first construct the solutions of the Poisson equation taking into account the electrostatic boundary conditions. Assuming the boundaries located at y = 0 and H do not carry free charges, the electrostatic boundary conditions read where The potential φ that satisfies these boundary conditions can be obtained using the method of electrostatic image charges [76]. We group the (infinite number of) image charges into two sets: image charges in group I L are obtained from first imaging the source charge distribution in the left boundary (y = 0), and then in the right boundary, and so on. For this group, the location of the image charges are On the other hand, image charges in group I R are obtained by first imaging the source charge distribution in the right boundary (y = H), and then in the left boundary and so on; in this case, the image locations are given by (Note that . . . and . . . are the floor and the ceiling functions, respectively·) For both groups, the nth image has electric charge Q n = λ n Q src where we have defined the dielectric contrast λ between the electrolyte solution and the boundary material as ( in and out are the permittivities of the solvent and the boundaries, respectively) which determines the electric charge ratio between successive images. The electric potential φ that satisfies the boundary conditions (38) is then calculated as where the first term is the electric potential created directly by the source charge distribution, and we have made use of the translation invariance of the system along the direction parallel to the surfaces (i.e., the s coordinates) by taking the Fourier transform of φ free with momentum k s ∈ R d−1 along them. The noise-averaged normal component of the Maxwell stress corresponding to this electric potential is then calculated as To calculate the (excess) stress at the location of the boundaries, we implement Eq. (42) and perform the summation over the image charges, noting that One can readily show that the stress at the location of the two plates are equal, and after some algebraic manipulation it can be written as × e −ksy e 2ksH e 2ksH − λ 2 + e ksy λ e 2ksH − λ 2 e −ksy e 2ksH e 2ksH − λ 2 + e ksy λ e 2ksH − λ 2 .
Note that this expression holds for a generic charge distribution which is invariant with respect to translations along the s coordinates (i.e., parallel to the surfaces) and has the charge correlation function ρ (2) (y , y ; k s ; t).

B. Stress exerted by the confined driven electrolyte
We now substitute the charge correlation function Eq. (36) into the stress formula in Eq. (45) and perform the integrations over y and y to obtain the stress for the specific settings at hand (see Fig. 1). Note that the resulting expressions contain formally divergent terms that will be removed upon subtracting the bulk stress (which is necessary to obtain the net force acting on the boundaries). The net force per unit area of the right boundary is then obtained as with the dimensionless amplitudes A defined as Here, τ = Dt/H 2 and we have also defined where λR n represents the electrostatic response from the image charges (note that the amplitude vanishes for λ = 0 where there is no such response). Moreover, we used the definition where ν s = (ν s1 = ν x , ν s2 , . . . , ν s d−1 ) ∈ R d−1 is a dimensionless vector and ν s = |ν s |.
To analyze the time-dependent stress amplitude given by Eq. (47), it is convenient to separate its steady-state and transient parts as where A s denotes the long-time steady-state amplitude and is given by the general expression and A τ is the transient part of the amplitude which reads The transient part A τ determines the initial behavior of the force amplitude after the quench, but it vanishes at long times. As will be shown in the following section, for τ 1, the temporal decay of A τ takes a power-law form (with possible crossover regimes) which corresponds to the long-time tails of diffusion processes.

C. Steady-state stress amplitude
The steady-state properties of the stress and the amplitude were investigated in Ref. [56]. Here, we present a brief account of the results. Performing the summations involved in Eq. (51) and after some simplifications we arrive at the following expression for A in d spatial dimensions where ν s cos θ = ν x and Y ± (λ, ν s ) are defined in Eq. (48). For d = 3, the some integrations in Eq. (53) can be carried out, yielding where Li n (z) = ∞ k=1 z k k n is the polylogarithm function. This expression can be evaluated by numerical methods, and we also report simple asymptotic expressions for its limiting cases in Table. I (for more details on the derivation of the asymptotic formulas and also similar analysis for the d = 2 case, see the supplemental material of Ref. [56]).

D. Transient stress amplitude
We now focus on the transient part of the stress amplitude as per Eq. (52). First, we note that the initial rate of change of A τ is given by where we have used ν s cos θ = ν x . Note that this initial rate is independent of the electric field, and it is a function of the dielectric contrast λ only. For d = 3, carrying out the integrals in Eq. (56) yields which can be approximated for small dielectric contrasts to second order in λ as Since the total stress at t = 0 vanishes, this shows that for 0 λ 0.1 the stress amplitude initially decreases and becomes negative (i.e., the force between the plates is initially repulsive). Having calculated the initial rate of A τ and using the asymptotic expansions of A s given in Table I which determine the long-term behavior of the stress, one can see there are a number of different dynamical behaviors that the FIF exhibits • for λ −0.31 and λ 0.17, both the initial slope ∂ τ A τ =0 and the stready-state amplitude A s are positive.
• for −0.31 λ < 0, the initial slope of A τ is positive, and A s is negative for weak electric fields while it becomes positive for strong fields.
• for 0 < λ 0.1, the initial slope is negative, and A s is positive for weak fields and negative for strong fields.
• for 0.1 λ 0.17, the initial slope is positive, and A s is positive for small electric fields while it is negative for strong fields.
To investigate how A τ decays at long times, we defineν s = √ τ ν s upon which Eq. (52) reads (In the second line we have definedñ = √ τ n for the integration variable.) Since the exponential factors suppress the integrands for large values ofν s , the final outcome of the integration is effectively determined by the behavior of the integrand forν s ∼ O(1). In addition, for τ → ∞, to leading order the summation is given by the n = 0 term. We thus only keep this contribution and also expand the exponential factor e πνs/ √ τ to obtain From the second line, it becomes evident that the transient force amplitude due to the bulk electrolyte outside the boundaries decays as ∼ τ −d/2 ; this is, in fact, the usual power-law tail of the diffusion process in d spatial dimensions.
On the other hand, the temporal decay of the transient stress coming from the electrolyte confined between the plates exhibits two different regimes: for τ π 2 /(1 − λ) 2 the decay is governed by the power-law form ∼ τ −(d−1)/2 , while for τ π 2 /(1 − λ) 2 it follows the form ∼ τ −(d+1)/2 . Since these expressions are obtained for τ 1, the first of these two regimes is only accessible when 1 − λ 1. (However, note that A τ may also exhibit a sign change at times comparable to the crossover time which in principle can make it difficult to observe the crossover between the two regimes.) The above analysis shows that at the longest time scales, the decay of the FIF amplitude is governed by the diffusive tails of the bulk electrolyte outside the boundaries (and not those of the confined electrolyte between the plates). It is worth mentioning that in Eq. (60), the sign of this asymptotic long-time behavior is determined by λ. A comparison with the (sign of the) steady-state amplitude A s shows that for values of dielectric contrast and with weak applied electric fields E 1, the total amplitude A overshoots A s at a finite time before approaching it at long times. This phenomenon also happens for λ 0.17 with strong applied fields. In Figs. 2 and 3, the temporal variations of the full FIF amplitude A and its transient part A τ as obtained from the numerical evaluation of Eqs. (51) and (52)  for other values of λ, the decay is governed by τ −3/2 . The long-time limit of the full amplitude A corresponds to the FIF at the steady state and can change sign for −0.17 λ 0.31 [56]. Note that for λ = 0.12, the sign change is not observable for E = 3 since it requires very large values of the electric field strength.

V. CONCLUDING REMARKS
In this work, we studied the correlation functions and the fluctuation effects in a strong electrolyte in the transient regime after an electric field quench that drives the solution out of thermal equilibrium state. The density and charge fluctuations are generally long-range correlated both in this transient period as well as in the long-time nonequilibrium steady state, as a result of the generic scale invariance of the stochastic dynamics. Such fluctuations, in turn, give rise to novel long-ranged forces on the confining boundaries. We analyze these forces as a function of the time elapsed from the electric field quench, which together with the steady-state result of Ref. [56] provide a complete account of the forces at different time scales. We find that the FIF scales with the plate separation as H −d in d spatial dimensions, and in general it has a diffusive character in its approach toward the steady-state form (see Eqs. (46) and (47)). This diffusive approach gives rise to power-law temporal decays of the transient part of the force at long times. The strength and the direction (i.e., attraction or repulsion) of the steady-state force can be controlled by the strength of the applied electric field. On the other hand, the early time temporal variations of the force amplitude can be non-monotonic and in some cases, it can result in changes in the sign of the force. These rich features point towards unexplored methods of force manipulation in practical applications, for example to control neutral colloidal particles that are immersed in an electrolyte solution. Moreover, they resemble some of the experimentally observed temporal patterns of force variation in surface measurements [29,31]; this implies that fluctuation effects which are generally discarded in mean-field models can indeed be relevant to understanding the force generation mechanisms in charged solutions out of equilibrium.
Although to derive the correlation functions, which form the basis of the force calculation, we have assumed a static external electric field, in Section II C we argued these results are also applicable when the electric field varies slowly over time. Since the time-scale of the charge relaxation is often of order of a few nanoseconds, this condition is, in fact, met in many experiments where oscillatory electric fields are used. The long-time description employed here can then be used by making minimal modifications as described below Eq. (24). It has been shown that an oscillatory driving field accompanied by different mobility coefficients for cations and anions can give rise to steady electric fields within the electrolyte [71]; however, once an electric field driving the ions in opposite directions is set up in the electrolyte, the analysis in Section II D shows that the long-distance asymptotic limit of the fluctuations is governed equations similar to the case of equal mobilities of charge species.
The present work uncovers a novel dynamical mechanism for generating long-ranged forces in driven charged fluids which had remained obscure within mean-field approaches. Although we have focused on a very simple setup here (constant electric field, flat boundaries, symmetric binary electrolyte), these results can in principle be extended to more complex settings since the underlying notion of generic scale invariance in anisotropic conserved dynamics pertains to a wide range of driven charged systems. It would be particularly relevant for experimental setups to investigate how oscillatory electric fields that are perpendicular to the confining boundaries would change the FIF studied here (note that static fields will lead to charge accumulation on the electrodes which eventually screen the applied field withing a short distance).

Appendix B: Equal-time correlation functions
In this appendix, we provide a detailed computation of the density and charge two-point correlation functions at equal times, directly from the linearized stochastic dynamics of Eqs. (10) and (11), both in the bulk as well as in flat confinement.
These coupled equations can directly be solved for c(k, ω) and ρ(k, ω) in terms of the noise fields η c and η ρ , from which the correlation functions read To compute the equal-time correlation functions, we need to perform the frequency integrations. For the density correlations, these can be expressed in the following form c(k, t)c(k , t) bulk = (2π) d δ(k + k ) 4DC 0 k 2 dω 2π where α = D 2 (k 2 + κ 2 ) 2 + E 2 D 2 κ 2 k 2 x (recall that E = µQE/(Dκ)), and ω ± u and ω ± l are defined as which represent the frequency poles in the upper half and lower half of the complex frequency plane, respectively. The frequency integration can be carried out, with the final result reading: Here, the first term in the bracket is independent of the applied electric field and represents the local correlations in equilibrium, while the second part vanishes for E = 0 and is the nonequilibrium part of the density correlations. In the limit of k/κ 1, the above bulk density correlation is approximated by recovering Eq. (28) which was obtained using the quasi-stationary approximation of Eq. (13). For the charge correlations, a similar calculation yields Similar to the density correlation functions, the first term in the brackets represents the short-ranged equilibrium correlations (note that in real space this gives the sum of a delta function and a screened-Coulomb (Yukawa-type) term which decays exponentially with Debye screening length κ −1 ); the second term, on the other hand, is the nonequilibrium part due to the external field, and it vanishes for a non-driven electrolyte (i.e., when E = 0). For k/κ 1, expanding the charge correlation function yields .

(B9)
Note that this expression agrees with the long-distance limit of Eq. (29) at late times.
Finally, the y-dependent correlation functions are obtained from these expression via c(k s , y)c(k s , y ) = 4 cos(p n y) cos(p n y ) c n (k s )c n (k s ) , with a similar expression for the charge correlations.
For k x > κ/(2E), the charge correlations are given by the following damped oscillatory expression: