1. Introduction
Electrolyte solutions, consisting of a polarized solvent and ionic species, are extremely important in a wide range of fields, including chemical physics, biology and geochemistry (Guglielmi et al. Reference Guglielmi, Voermans, Gualandi, Van Engelen, Ferlini, Tomelleri and Vattemi2013; Molina-Osorio et al. Reference Molina-Osorio, Gamero-Quijano, Peljo and Scanlon2020). They are also of interest in engineering physics, such as aggregation and dispersion of charged liquid particles (Boutsikakis et al. Reference Boutsikakis, Fede, Pedrono and Simonin2020; Nieto et al. Reference Nieto, Agrawal, Bravo, Hofmeister-Mock, Pepi and Ghoshal2021) present in aerospace systems. Despite a long history of electrolyte solution studies, there are still important open questions associated with fluctuations and correlations of electrolyte bulk solutions (Donev et al. Reference Donev, Garcia, Péraud, Nonaka and Bell2019a). At the mesoscale (i.e. nanometre to micrometre length scales), thermal energies of electrolyte solutions are of the same magnitude as the characteristic energies of hydrodynamics and electrokinetics. Therefore, thermally induced fluctuations can play an important role in both equilibrium and non-equilibrium electrokinetic phenomena at the mesoscale (Ladiges et al. Reference Ladiges, Nonaka, Klymko, Moore, Bell, Carney, Garcia, Natesh and Donev2021). Quantifying the impact of fluctuations on mesoscale fluid systems is critical to understanding the large-scale dynamics of complex fluid systems designed from the nano/mesoscale in a bottom-up approach.
The essential feature of an electrolyte solution is that charged species interact with each other with long-range Coulomb forces, leading to a system whose properties are significantly different from those of electrically neutral solutions (Klymko et al. Reference Klymko, Nonaka, Bell, Carney and Garcia2020). The presence of a liquid–liquid or liquid–solid interface will change the dynamics of charged particles, and the behaviour of these charged particles near surfaces could differ significantly from the bulk properties as they form condensed ion layers (Larsen & Grier Reference Larsen and Grier1997) and ionic double-layer structures (Sidhu, Frischknecht & Atzberger Reference Sidhu, Frischknecht and Atzberger2018). Elimination of sound waves in low-Mach-number transport of charged species at the mesoscale with significant thermal fluctuations can yield a quasi-incompressible formulation of fluctuating hydrodynamics (Péraud et al. Reference Péraud, Nonaka, Chaudhri, Bell, Donev and Garcia2016). Discrete ion effects can induce depletion forces and significant nanoparticle–wall interactions (Asakura & Oosawa Reference Asakura and Oosawa1954, Reference Asakura and Oosawa1958). The differential capacitance of charged particles can vary from a minimum value for aqueous solution (Kilic, Bazant & Ajdari Reference Kilic, Bazant and Ajdari2007) to a maximum value for molten electrolytes (Lamperski & Kłos Reference Lamperski and Kłos2008). The changes from minimum to maximum values are driven by a change of reduced charge density from small, for aqueous solutions, to large, for molten electrolytes (Lamperski, Outhwaite & Bhuiyan Reference Lamperski, Outhwaite and Bhuiyan2009). In the electroneutral limit, the fluctuating Poisson–Nernst–Planck (PNP) equations are constrained to preserve charge neutrality by a variable-coefficient elliptic equation instead of the standard Poisson equation (Donev et al. Reference Donev, Nonaka, Kim, Garcia and Bell2019b). At the continuum scale level, fluctuating hydrodynamic and electrokinetic theories have been used to describe thermally induced fluctuations through random tensor/flux terms in the governing equations, formulated properly to satisfy the fluctuation–dissipation theorem (Landau & Lifshitz Reference Landau and Lifshitz1959; Ortiz & Sengers Reference Ortiz and Sengers2006). The electrostatic properties can be modelled by the Poisson–Boltzmann (PB) equation (Baker et al. Reference Baker, Sept, Joseph, Holst and McCammon2001). At the atomic scale, molecular dynamics (MD) simulations with explicit ions have also been used to study electrolyte solutions. MD models of electrolyte fluids can maintain a good accuracy compared to the classical density functional theory (cDFT) calculations, while PB models depart significantly from cDFT and MD at high charge densities (Lee et al. Reference Lee, Nilson, Templeton, Griffiths, Kung and Wong2012). However, the huge computational cost of MD simulations prevents their use at large length and time scales (Chen & Pappu Reference Chen and Pappu2007; Joung, Luchko & Case Reference Joung, Luchko and Case2013; Yoshida et al. Reference Yoshida, Mizuno, Kinjo, Washizu and Barrat2014). Mesoscale approaches smoothly bridge the gap between continuum and atomistic descriptions, and provide the possibility to integrate fluctuations consistently into macroscopic field variables. To this end, we developed a charged dissipative particle dynamics (cDPD) model (Deng et al. Reference Deng, Li, Borodin and Karniadakis2016) to simulate mesoscale electrokinetic phenomena with fluctuations.
Hydrodynamic fluctuations at mesoscopic scales have been studied extensively in the past decade (Bian et al. Reference Bian, Li, Deng and Karniadakis2015; Péraud et al. Reference Péraud, Nonaka, Chaudhri, Bell, Donev and Garcia2016; Yu et al. Reference Yu, Eckmann, Ayyaswamy and Radhakrishnan2016; Bian, Deng & Karniadakis Reference Bian, Deng and Karniadakis2018; Kim et al. Reference Kim, Nonaka, Bell, Garcia and Donev2018; Donev et al. Reference Donev, Garcia, Péraud, Nonaka and Bell2019a). However, the electrokinetic fluctuations in electrolyte solutions have been explored less, especially from a theoretical perspective. In the present work, we focus on theoretical analysis on electrokinetic fluctuations in bulk electrolyte solutions close to the thermodynamic equilibrium state, aiming to provide closed-form expressions for temporal correlation functions of electrokinetic fluctuations. We consider a continuum formulation and derive the coupled fluctuating hydrodynamic and electrokinetic equations. Based on the linear response assumption (Kubo Reference Kubo1982), we will linearize the fluctuating hydrodynamic and electrokinetic equations and derive analytical closed-form expressions for current correlation functions of electrolyte bulk solutions in Fourier space. Then we will perform both all-atom MD simulations with explicit ions and mesoscale cDPD simulations with semi-implicit ions to compute spatial and temporal correlations of hydrodynamic and electrostatic fluctuations directly from particle trajectories. These simulation results will be used to validate the closed-form expressions for correlation functions derived from the linearized fluctuating hydrodynamics theory by comparing simulation results against theoretical predictions.
The remainder of this paper is organized as follows. In § 2, we introduce the continuum formulation for fluctuating hydrodynamics and electrokinetics, and also the derivations of analytical solutions for current correlation functions using perturbation theory. In § 3, we describe the details for performing all-atom MD simulations and mesoscopic cDPD simulations, and we also present the simulation results comparing against the theoretical predictions. Finally, we conclude with a brief summary and discussion in § 4.
2. Continuum theory
2.1. Fluctuating hydrodynamics and electrokinetics
We consider a mesoscale system of electrolyte solution at thermal equilibrium in a periodic domain 
$\boldsymbol {\varOmega }$ of fixed volume. The system contains 
$S$ types of ionic species with concentration 
$c_\alpha (\boldsymbol {r},t)$ and charge valency 
$z_\alpha$ for the 
$\alpha$th ionic species. Let 
$e$ be the elementary charge. Then a global constraint is imposed by the charge neutrality condition
\begin{equation} \sum_{\alpha=1}^{S} \int_{\boldsymbol{\varOmega}} e z_\alpha\,c_\alpha(\boldsymbol{r},t)\,{{\rm d}}\boldsymbol{r} = 0. \end{equation}
The solvent molecules are represented implicitly through their electrostatic and thermodynamic properties, such as dielectric constant 
$\epsilon$, bulk viscosity 
$\zeta$, and shear viscosity 
$\eta$. From the continuum perspective, this system can be described by equations of classical fluctuating hydrodynamics with an additional electrostatic force:
$$\begin{gather} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{g} = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \boldsymbol{g}}{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{g}\boldsymbol{v}) ={-} \boldsymbol{\nabla} p + \eta\,\nabla^{2}\boldsymbol{v} + \left(\frac{\eta}{3} + \zeta\right) \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{v}) + \rho_e \boldsymbol{E} + \boldsymbol{\nabla} \boldsymbol{\cdot} \delta \boldsymbol{\boldsymbol{\varPi}}, \end{gather}$$
for velocity 
$\boldsymbol {v}(\boldsymbol {r}, t)$, pressure 
$p(\boldsymbol {r},t)$, mass density 
$\rho (\boldsymbol {r},t)$, momentum density 
$\boldsymbol {g}(\boldsymbol {r},t) = \rho (\boldsymbol {r},t)\, \boldsymbol {v}(\boldsymbol {r},t)$, and electric field 
$\boldsymbol {E}(\boldsymbol {r},t)$. The local charge density is given by 
$\rho _e(\boldsymbol {r},t) = \sum _{\alpha =1}^{S} e z_\alpha \, c_\alpha (\boldsymbol {r},t)$. The random stress tensor 
$\delta \boldsymbol {\varPi }$ is a matrix of Gaussian-distributed random variables with zero means and variances given by the fluctuation–dissipation theorem (FDT) (Ortiz & Sengers Reference Ortiz and Sengers2006)
\begin{equation} \langle \delta \varPi_{ij}(\boldsymbol{r},t)\,\delta \varPi_{kl}(\boldsymbol{r}',t')\rangle = 2 k_BT \boldsymbol{\mathsf{C}}_{ijkl}\, \delta(\boldsymbol{r}-\boldsymbol{r}')\,\delta(t-t'), \end{equation}
where 
$\boldsymbol{\mathsf{C}}_{ijkl}=\eta (\delta _{ik}\delta _{jl} + \delta _{il}\delta _{jk})+(\zeta -{2\eta }/{3}) \delta _{ij}\delta _{kl}$ is a rank-4 tensor. The fluctuating hydrodynamic equations are closed with the equation of state, e.g. 
$c_s^{2} = (\partial p/\partial \rho )_T$, with 
$c_s$ being the isothermal sound speed.
The Ginzburg–Landau free energy functional for an electrolyte solution is (Ginzburg & Landau Reference Ginzburg and Landau1950)
\begin{equation} G = \int_{\boldsymbol{\varOmega}} {{\rm d}}\boldsymbol{r} \left\{k_BT \sum_{\alpha=1}^{S} c_\alpha \ln c_\alpha + \sum_{\alpha=1}^{S} e z_\alpha c_\alpha \phi - \frac{\epsilon}{2}\,(\boldsymbol{\nabla}\phi)^{2}\right\}, \end{equation}
where 
$k_BT$ is the thermal energy, and 
$\phi$ is the electrostatic potential. Variation of the free energy with respect to the electrostatic potential yields the Poisson equation
\begin{equation} {-}\boldsymbol{\nabla}\boldsymbol{\cdot} (\epsilon(\boldsymbol{r})\,\boldsymbol{\nabla}\phi) = \sum_{\alpha=1}^{S} e z_\alpha\,c_\alpha(\boldsymbol{r},t), \end{equation}
by setting 
$\delta G/\delta \phi = 0$. Similarly, the electrochemical potential of the 
$\alpha$th ionic species can be derived by variation with respect to ionic concentration:
\begin{equation} \mu_\alpha = \frac{\delta G}{\delta c_\alpha} = k_BT \ln c_\alpha + z_\alpha e\phi. \end{equation}
The transport and dissipation of ionic species are driven by the fluid velocity, electrochemical potential and thermal fluctuations, which can be described in terms of the ionic concentration flux 
$\boldsymbol {J}(\boldsymbol {r},t)$ by
\begin{equation} \frac{\partial c_\alpha(\boldsymbol{r},t)}{\partial t} + \boldsymbol{v}(\boldsymbol{r},t) \boldsymbol{\cdot}\boldsymbol{\nabla} c_\alpha(\boldsymbol{r},t) ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{J}_{\alpha}(\boldsymbol{r},t) + \delta \boldsymbol{J}_{\alpha}(\boldsymbol{r},t) ), \end{equation}
where 
$\boldsymbol {J}_{\alpha }(\boldsymbol {r},t)$ is the dissipative flux, and 
$\delta \boldsymbol {J}_{\alpha }(\boldsymbol {r},t)$ is the random flux when the system is near thermodynamic equilibrium. The diffusion flux can be written as 
$\boldsymbol {J}_\alpha = - \sum _{\beta =1}^{S} M_{\alpha \beta }\,\boldsymbol {\nabla } \mu _\beta (\boldsymbol {r},t)$, in which 
$\boldsymbol {\nabla } \mu$ is the thermodynamic force for diffusion flux, and 
$M_{\alpha \beta }$ are the Onsager coefficients related to macroscopic ionic diffusion coefficients. In general, 
$M_{\alpha \beta } = M_{\beta \alpha } \neq 0$, as implied by reversal invariance (Onsager Reference Onsager1931a,Reference Onsagerb). The off-diagonal terms 
$M_{\alpha \neq \beta }$ describe mutual diffusion and are assumed to be relatively small compared to the diagonal self-diffusion terms. For example, according to experimental data (Chapman Reference Chapman1967), the self-diffusion coefficients of cation and anion in a 
$1\,{\rm M}$ NaCl solution are 
$1.16\times 10^{-9}$ and 
$1.99\times 10^{-9}\,{\rm m}^{2}\,{\rm s}^{-1}$, respectively, while the mutual diffusion coefficient is one order of magnitude smaller: 
$1.3\times 10^{-10}\,{\rm m}^{2}\,{\rm s}^{-1}$. Therefore, we assume that the mutual diffusion terms can be ignored in the present work for simplicity. In numerical experiments, this assumption is checked by computing the self-diffusion and mutual diffusion coefficients from MD data using the Green–Kubo relations (Zhou & Miller Reference Zhou and Miller1996; Wheeler & Newman Reference Wheeler and Newman2004), as presented in § 3.
By substituting the electrochemical potential expression into the diffusion flux, we obtain
\begin{equation} \boldsymbol{J}_\alpha ={-} M_{\alpha}\,\boldsymbol{\nabla} \mu_\alpha (\boldsymbol{r},t) ={-} \frac{M_{\alpha} k_BT}{c_\alpha}\,\boldsymbol{\nabla} c_\alpha - M_{\alpha}z_\alpha e\,\boldsymbol{\nabla} \phi. \end{equation}
In practice, it is more convenient to use the macroscopic diffusion coefficients 
$D_\alpha = {M_{\alpha } k_BT}/{c_\alpha }$ instead of the phenomenological coefficients 
$M_{\alpha }$. Therefore, the ionic concentration transport equation can be rewritten as
\begin{equation} \frac{\partial c_\alpha(\boldsymbol{r},t)}{\partial t} + \boldsymbol{v}(\boldsymbol{r},t) \boldsymbol{\cdot}\boldsymbol{\nabla} c_\alpha(\boldsymbol{r},t)= \boldsymbol{\nabla}\boldsymbol{\cdot} \left(D_\alpha\,\boldsymbol{\nabla} c_\alpha + \frac{D_\alpha e z_\alpha c_\alpha}{k_BT}\,\boldsymbol{\nabla} \phi +\delta \boldsymbol{J}_{\alpha}(\boldsymbol{r},t)\right). \end{equation}For a system near thermodynamic equilibrium, the random fluxes can be modelled as Gaussian random vectors with zero means and variance given by the generalized FDT as
\begin{align} \langle \delta \boldsymbol{J}_{\alpha} (\boldsymbol{r},t) \boldsymbol{\cdot} \delta\boldsymbol{J}_{\beta}(\boldsymbol{r}',t')\rangle &= 2k_BT M_{\alpha \beta}\,\delta (\boldsymbol{r} -\boldsymbol{r}')\,\delta(t-t') \nonumber\\ &= 2 D_{\alpha} c_{\alpha}\,\delta (\boldsymbol{r}-\boldsymbol{r}')\,\delta(t-t'). \end{align}The coupled equations (2.2)–(2.10) form the fluctuating hydrodynamic and electrokinetic equations.
2.2. Linearized theory of electrolyte solution
The stochastic partial differential equations (sPDEs) given by (2.2)–(2.10) could be solved through numerical discretization, with the FDT satisfied at the discrete level following the GENERIC framework (Grmela & Öttinger Reference Grmela and Öttinger1997; Öttinger & Grmela Reference Öttinger and Grmela1997). In general, it is very challenging to obtain an analytical solution of such sPDEs. However, if the local hydrodynamic and electrokinetic fluctuations are sufficiently small, then linear response theory can be applied to derive linearized equations to describe the relaxation process of electrolyte solution towards equilibrium. In the present work, we focus on the linearized equations for electrolyte solutions and their closed-form solutions.
The equilibrium state of a bulk electrolyte solution is characterized by its mean field properties, i.e. constant mass density 
$\rho _0$, constant pressure 
$p_0$, constant bulk ionic concentration 
$c_{\alpha 0}$, zero momentum field 
$\boldsymbol {g}_0 = 0$, and zero electrostatic potential field 
$\phi _0 = 0$. The local fluctuating hydrodynamic field can be expressed as the perturbation around the mean field state:
$$\begin{gather} \rho(\boldsymbol{r},t) = \rho_0 + \delta \rho(\boldsymbol{r},t), \end{gather}$$
$$\begin{gather}p(\boldsymbol{r},t) = p_0 + \delta p(\boldsymbol{r},t), \end{gather}$$
$$\begin{gather}\boldsymbol{g}(\boldsymbol{r},t) = \delta \boldsymbol{g} (\boldsymbol{r},t), \end{gather}$$
where the local perturbation of pressure can be related to density fluctuations via the equation of state 
$\delta p = c_s^{2}\,\delta \rho$ under isothermal conditions. Also, the local fluctuating electrostatic field can be decomposed as
$$\begin{gather} c_\alpha(\boldsymbol{r},t) = c_{\alpha 0} + \delta c_\alpha(\boldsymbol{r},t), \end{gather}$$
$$\begin{gather}\phi(\boldsymbol{r},t) = \delta \phi(\boldsymbol{r},t), \end{gather}$$
$$\begin{gather}\rho_e(\boldsymbol{r},t) = \delta \rho_e (\boldsymbol{r},t), \end{gather}$$
in which 
$\delta \rho _e(\boldsymbol {r},t) = \sum _{\alpha =1}^{S} e z_\alpha \,\delta c_\alpha (\boldsymbol {r},t)$ when the global electroneutrality condition 
$\sum _{\alpha =1}^{S} e z_\alpha c_{\alpha 0} = 0$ is imposed. The fluctuating electrostatic potentials and charge densities are related via the Poisson equation.
For small fluctuations of electrostatic field, the linearized fluctuating hydrodynamic and electrokinetic equations can be written as
$$\begin{gather} \frac{\partial (\delta \rho)}{\partial t} ={-}\boldsymbol{\nabla}\boldsymbol{\cdot} \delta \boldsymbol{g}, \end{gather}$$
$$\begin{gather}\frac{\partial (\delta \boldsymbol{g})}{\partial t} ={-}{c_s^{2}}\,\boldsymbol{\nabla} \delta \rho + \eta\,\nabla^{2} (\delta \boldsymbol{v}) + \left(\frac{\eta}{3}+\zeta\right) \boldsymbol{\nabla} (\boldsymbol{\nabla}\boldsymbol{\cdot}\delta \boldsymbol{v} ), \end{gather}$$
$$\begin{gather}\frac{\partial (\delta c_\alpha)}{\partial t} = \left(D_{\alpha}\,\nabla^{2} (\delta c_{\alpha}) + D_{\alpha}\, \frac{e z_{\alpha} c_{\alpha 0}}{k_BT}\,\nabla^{2} \delta \phi \right), \end{gather}$$
$$\begin{gather}-\boldsymbol{\nabla}\boldsymbol{\cdot} (\epsilon\,\boldsymbol{\nabla} \delta \phi) = \sum_{\alpha=1}^{S} e z_\alpha\,\delta c_\alpha(\boldsymbol{r},t), \end{gather}$$
where only first-order perturbation terms are kept in these linearized equations. It is important to note that the advection term 
$\boldsymbol {\nabla }\boldsymbol {\cdot }(\delta \boldsymbol {g}\,\delta \boldsymbol {v})$ and the electrostatic force 
$\delta \rho _e\, \boldsymbol {\nabla } \delta \phi$ in the momentum equation, and the convection term 
$\delta \boldsymbol {v} \boldsymbol {\cdot }\boldsymbol {\nabla } \delta c_\alpha$ in the transport equation, are high-order perturbation terms and are thus assumed to be negligible in the equations above. Therefore, the linearized hydrodynamic and electrokinetic equations become explicitly decoupled for fluid systems in thermodynamic equilibrium.
The linearized fluctuating hydrodynamic equations given by (2.13) can be transformed into 
$k$-space by a spatial Fourier transform and then solved in the 
$k$-space. Appendix A describes the detailed derivation of mass–momentum correlations by solving (2.13). The normalized temporal correlation functions of the mass–momentum fluctuations in the 
$k$-space are given by
$$\begin{gather} \frac{\langle \hat{\rho}(\boldsymbol{k},t)\,\hat{\rho}(\boldsymbol{k},0)\rangle}{\langle \hat{\rho}(\boldsymbol{k},0)\,\hat{\rho}(\boldsymbol{k},0)\rangle} = \exp(-\varGamma_T k^{2} t) \cos(c_s k t), \end{gather}$$
$$\begin{gather}\frac{\langle \hat{\boldsymbol{g}}_\parallel(\boldsymbol{k},t)\,\hat{\boldsymbol{g}}_\parallel(\boldsymbol{k},0)\rangle}{\langle \hat{\boldsymbol{g}}_\parallel (\boldsymbol{k},0)\,\hat{\boldsymbol{g}}_\parallel(\boldsymbol{k},0)\rangle} = \exp(-\varGamma_T k^{2} t) \cos(c_s k t), \end{gather}$$
$$\begin{gather}\frac{\langle\hat{\boldsymbol{g}}_\perp(\boldsymbol{k},t)\,\hat{\boldsymbol{g}}_\perp(\boldsymbol{k},0)\rangle}{\langle \hat{\boldsymbol{g}}_\perp (\boldsymbol{k},0)\,\hat{\boldsymbol{g}}_\perp(\boldsymbol{k},0)\rangle} = \exp(-\nu k^{2} t), \end{gather}$$
$$\begin{gather}\frac{\langle \hat{\rho}(\boldsymbol{k},t)\,{\rm i}\, \hat{\boldsymbol{g}}_\parallel(\boldsymbol{k},0)\rangle}{\langle \hat{\rho} (\boldsymbol{k},0)\,{\rm i}\, \hat{\boldsymbol{g}}_\parallel(\boldsymbol{k},0)\rangle} = \exp(-\varGamma_T k^{2} t) \sin(c_s k t), \end{gather}$$
where the symbol 
$\hat {~}$ indicates Fourier components, the wave vector 
$\boldsymbol {k} = (k, 0, 0)$ is defined along the (arbitrary) 
$x$-direction, 
$\nu = \eta /\rho$ is the kinematic viscosity, 
$\varGamma _T=2\nu /3+\zeta /2\rho$ is the sound absorption coefficient, and 
$c_s$ is the isothermal sound speed. Also, 
$\boldsymbol {g}_\parallel$ represents longitudinal momentum (sound mode) parallel to the wave vector 
$\boldsymbol {k}$, and 
$\boldsymbol {g}_\perp$ represents the transverse component (shear mode) perpendicular to the wave vector 
$\boldsymbol {k}$.
For simplicity in deriving closed-form solutions, we consider only two types of ionic species: the cation (denoted by 
$p$) and anion (denoted by 
$n$). By substitution of the Poisson equation into the linearized ionic transport equation, we obtain
$$\begin{gather} \frac{\partial(\delta c_p)}{\partial t} = D_p \left(\nabla^{2} (\delta c_p) + \kappa_p^{2}\,\delta c_p - \left| \frac{z_n}{z_p} \right| \kappa_p^{2}\,\delta c_n \right), \end{gather}$$
$$\begin{gather}\frac{\partial(\delta c_n)}{\partial t} = D_n \left(\nabla^{2} (\delta c_n) + \kappa_n^{2}\,\delta c_n - \left| \frac{z_p}{z_n} \right| \kappa_n^{2}\,\delta c_p \right), \end{gather}$$
where 
$\kappa _p^{2}$ and 
$\kappa _n^{2}$ are defined as
\begin{equation} \kappa_p^{2} \equiv \frac{D_p {z_p}^{2} c_{p 0}e^{2}}{k_BT \epsilon},\quad \kappa_n^{2} \equiv \frac{D_n {z_n}^{2} c_{n 0}e^{2}}{k_BT \epsilon}, \end{equation} with units of inverse-squared distance. For a periodic system, we can expand the solution in Fourier modes: 
$\delta c(\boldsymbol {r},t) = \sum _{\boldsymbol {k}}\delta \hat {c}(\boldsymbol {k},t) \,{\rm e}^{{\rm i}\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {r}}$. The spatial Fourier transformation of (2.15) gives two coupled ordinary differential equations
\begin{equation} \left.\begin{gathered} \frac{\partial \delta \hat{c}_p(\boldsymbol{k},t)}{\partial t} = D_p \left[(\kappa_p^{2}-k^{2})\,\delta \hat{c}_p(\boldsymbol{k},t) - \left| \frac{z_n}{z_p} \right| \kappa_p^{2}\,\delta \hat{c}_n(\boldsymbol{k},t)\right], \\ \frac{\partial \delta \hat{c}_n(\boldsymbol{k},t)}{\partial t} = D_n\left[(\kappa_n^{2}-k^{2})\,\delta \hat{c}_n(\boldsymbol{k},t) - \left|\frac{z_p}{z_n} \right| \kappa_n^{2}\,\delta \hat{c}_p(\boldsymbol{k},t)\right]. \end{gathered}\right\} \end{equation} The above equations can be written in matrix form as 
${\partial \boldsymbol {u}(\boldsymbol {k},t)}/{\partial t} + \boldsymbol{\mathsf{L}}\boldsymbol {u} = 0$, with the vector 
$\boldsymbol {u} = (\delta \hat {c}_p, \delta \hat {c}_n)^\textrm {T}$, and 
$\boldsymbol{\mathsf{L}}$ the 
$2\times 2$ matrix defined as
\begin{equation} \boldsymbol{\mathsf{L}} =
\left[\begin{array}{@{}cc@{}} (k^{2} - \kappa_p^{2})D_p &
\left|\dfrac{z_n}{z_p}\right|\kappa_p^{2} D_p \\
\left|\dfrac{z_p}{z_n}\right|\kappa_n^{2} D_n & (k^{2} -
\kappa_n^{2}) D_n \\ \end{array}\right].
\end{equation}
The matrix 
$\boldsymbol{\mathsf{L}}$ can be further split as 
$\boldsymbol{\mathsf{L}} = -\boldsymbol{\mathsf{L}}_0 + k^{2} \boldsymbol{\mathsf{L}}_1$, where
$$\begin{gather} \boldsymbol{\mathsf{L}}_0 = \left[ \begin{array}{@{}cc@{}} \kappa_p^{2} D_p & -\left|\dfrac{z_n}{z_p}\right|\kappa_p^{2} D_p \\ -\left|\dfrac{z_p}{z_n}\right|\kappa_n^{2} D_n & \kappa_n^{2} D_n \end{array}\right], \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{L}}_1 = \left[\begin{array}{@{}cc@{}} D_p & 0 \\ 0 & D_n \end{array}\right]. \end{gather}$$
These equations can be solved by a linear combination of the eigenvectors 
$\xi ^{(i)}(\boldsymbol {k})$, which satisfy the eigenvalue equation
\begin{equation} [-\boldsymbol{\mathsf{L}}_0 + k^{2} \boldsymbol{\mathsf{L}}_1]\,\xi^{(i)}(\boldsymbol{k}) = \lambda_i\,\xi^{(i)}(\boldsymbol{k}). \end{equation}
The conditions 
$\kappa _{p, n}^{2} \gg k^{2}$ hold in the continuum limit, when either the domain size 
$L$ or the charge concentration 
$z_\alpha ^{2} c_{\alpha 0}$ is sufficiently large. In this limit, the eigenvalue equation can be solved perturbatively by expanding 
$\xi ^{(i)}$ and 
$\lambda ^{(i)}$ in powers of 
$k^{2}$:
$$\begin{gather} \xi^{(i)} = \xi_0^{(i)} + k^{2} \xi_1^{(i)} + \cdots , \end{gather}$$
$$\begin{gather}\lambda^{(i)} = \lambda_0^{(i)} + k^{2} \lambda_1^{(i)} + \cdots. \end{gather}$$Substitution of (2.22) into the eigenvalue equation (2.21) gives the zero- and second-order perturbation theory equations:
$$\begin{gather} (\boldsymbol{\mathsf{L}}_0 - \lambda_0^{(i)}\boldsymbol{\mathsf{I}})\xi_0^{(i)} = 0, \end{gather}$$
$$\begin{gather}(\boldsymbol{\mathsf{L}}_0 - \lambda_0^{(i)}\boldsymbol{\mathsf{I}})\xi_1^{(i)} = (\boldsymbol{\mathsf{L}}_1 + \lambda_1^{(i)}\boldsymbol{\mathsf{I}}) \xi_0^{(i)}, \end{gather}$$
where 
$\boldsymbol{\mathsf{I}}$ denotes the identity matrix. The solution to order 
${O}(k^{2})$ is given by two real negative roots approximated as 
$\lambda _1 \approx -((k^{2}+\kappa _p^{2})D_p + (k^{2}+\kappa _n^{2})D_n)+k^{2}\, D_s(k)$ (fast decay) and 
$\lambda _2 \approx -k^{2}\,D_s(k)$ (slow decay), where the collective diffusion coefficient of cation and anion pairs is defined as
\begin{equation} D_s(k) = \frac{(k^{2} + \kappa_p^{2} + \kappa_n^{2})D_p D_n}{(k^{2}+\kappa_p^{2})D_p + (k^{2}+\kappa_n^{2})D_n}, \end{equation}
which has the same order as the diffusion coefficients 
$D_{\alpha }$. The solutions of (2.17) with initial conditions 
$\delta \hat {c}_p(\boldsymbol {k},0)$ and 
$\delta \hat {c}_n(\boldsymbol {k},0)$ can be written as
$$\begin{gather} \delta \hat{c}_p(\boldsymbol{k},t) = A_{1}\,\delta \hat{c}_p(\boldsymbol{k}, 0) + A_{2}\,\delta \hat{c}_n(\boldsymbol{k}, 0), \end{gather}$$
$$\begin{gather}\delta \hat{c}_n(\boldsymbol{k},t) = A_{3}\,\delta \hat{c}_n(\boldsymbol{k}, 0) + A_{4}\,\delta \hat{c}_p(\boldsymbol{k}, 0), \end{gather}$$with the coefficients
$$\begin{gather} A_{1} = \alpha_p \,{\rm e}^{\lambda_1 t} + \alpha_n \,{\rm e}^{\lambda_2 t}, \end{gather}$$
$$\begin{gather}A_{2} = \beta_p({\rm e}^{\lambda_2 t} - {\rm e}^{\lambda_1 t}), \end{gather}$$
$$\begin{gather}A_{3} = \alpha_n \,{\rm e}^{\lambda_1 t} + \alpha_p \,{\rm e}^{\lambda_2 t}, \end{gather}$$
$$\begin{gather}A_{4} = \beta_n({\rm e}^{\lambda_2 t} - {\rm e}^{\lambda_1 t}), \end{gather}$$where the dimensionless parameters are defined as
$$\begin{gather} \alpha_p = \frac{(k^{2} + \kappa_p^{2})D_p - k^{2}\,D_s(k)}{(k^{2} + \kappa_p^{2})D_p + (k^{2} +\kappa_n^{2})D_n - 2 k^{2}\,D_s(k)}, \end{gather}$$
$$\begin{gather}\alpha_n = \frac{(k^{2} + \kappa_n^{2})D_n - k^{2}\,D_s(k)}{(k^{2} + \kappa_p^{2})D_p + (k^{2} +\kappa_n^{2})D_n - 2 k^{2}\,D_s(k)}, \end{gather}$$
$$\begin{gather}\beta_p = \frac{\left|\dfrac{z_n}{z_p}\right|\kappa_p^{2} D_p}{(k^{2} + \kappa_p^{2})D_p + (k^{2} + \kappa_n^{2})D_n - 2 k^{2}\,D_s(k)}, \end{gather}$$
$$\begin{gather}\beta_n = \frac{\left|\dfrac{z_p}{z_n}\right|\kappa_n^{2} D_n}{(k^{2} + \kappa_p^{2})D_p + (k^{2} + \kappa_n^{2})D_n - 2 k^{2}\,D_s(k)}, \end{gather}$$
such that 
$\alpha _p + \alpha _n = 1$.
The fluctuations modelled by the above equations are transported and dissipated in time, and this can be shown via the time correlation of the fluctuating concentration field with different Fourier modes. The temporal correlations between species 
$\alpha$ and 
$\beta$ in the fluctuating concentration field are given by
\begin{equation} \varPsi_{\alpha \beta} = \frac{\langle \delta \hat{c}_\alpha(\boldsymbol{k},t)\, \delta \hat{c}_\beta^{*}(\boldsymbol{k},0)\rangle}{\langle \delta \hat{c}_\alpha(\boldsymbol{k},0)\,\delta \hat{c}_\beta^{*}(\boldsymbol{k},0)\rangle}, \end{equation}
where the symbol 
$^{*}$ denotes the complex conjugate. The temporal correlation function for specific ionic species can then be expressed as
$$\begin{gather} \varPsi_{p p} = \alpha_p \,{\rm e}^{\lambda_1 t} + \alpha_n \,{\rm e}{}^{\lambda_2 t} + \beta_p S_{n p} ({\rm e}^{\lambda_2 t} - {\rm e}^{\lambda_1 t}), \end{gather}$$
$$\begin{gather}\varPsi_{n n} = \alpha_n \,{\rm e}^{\lambda_1 t} + \alpha_p \,{\rm e}{}^{\lambda_2 t} + \beta_n S_{p n} ({\rm e}^{\lambda_2 t} - {\rm e}^{\lambda_1 t}), \end{gather}$$
$$\begin{gather}\varPsi_{p n} = \alpha_p \,{\rm e}{}^{\lambda_1 t} + \alpha_n \,{\rm e}{}^{\lambda_2 t} + \frac{\beta_p}{S_{p n}}\, ({\rm e}^{\lambda_2 t} - {\rm e}^{\lambda_1 t}), \end{gather}$$
$$\begin{gather}\varPsi_{n p} = \alpha_n \,{\rm e}{}^{\lambda_1 t} + \alpha_p \,{\rm e}{}^{\lambda_2 t} + \frac{\beta_n}{S_{n p}}\,({\rm e}^{\lambda_2 t} - {\rm e}^{\lambda_1 t}), \end{gather}$$
in which 
$S_{p n}$ and 
$S_{n p}$ are static structure factor-like terms for a wave vector 
$\boldsymbol {k}$ defined as
$$\begin{gather} S_{p n}(\boldsymbol{k}) = \frac{\langle \delta \hat{c}_p(\boldsymbol{k},0)\, \delta \hat{c}_n^{*}(\boldsymbol{k},0)\rangle}{\langle\delta\hat{c}_n(\boldsymbol{k},0)\, \delta \hat{c}_n^{*}(\boldsymbol{k},0)\rangle}, \end{gather}$$
$$\begin{gather}S_{n p}(\boldsymbol{k}) = \frac{\langle \delta \hat{c}_n(\boldsymbol{k},0)\, \delta \hat{c}_p{*}(\boldsymbol{k},0)\rangle}{\langle\delta\hat{c}_p(\boldsymbol{k},0)\, \delta \hat{c}_p{*}(\boldsymbol{k},0)\rangle}. \end{gather}$$
In this linearized theory, the covariance between the initial ionic concentration fluctuations of cation and anion are non-zero due to the charge neutrality constraint. The structure factor terms 
$S_{p n}$ and 
$S_{n p}$ can be obtained from experiment or more detailed simulations. The long-time behaviour of both temporal autocorrelations and cross-correlations are dominated by slowly decaying terms that behave asymptotically as 
$\exp (\lambda _2 t)$, where 
$\lambda _2$ is one of the two (negative) eigenvalues.
3. Numerical results
To validate the closed-form expressions of fluctuation correlations that we derived based on linearized fluctuating hydrodynamic and electrokinetic equations in § 2, we perform both all-atom MD simulations with explicit ions and mesoscale cDPD simulations with semi-implicit ions, as illustrated in figure 1. We first carry out all-atom MD simulations of an aqueous NaCl solution in the bulk, which consists of 400 Na
$^{+}$ ions, 400 Cl
$^{-}$ ions, and 64 000 H
$_2$O water molecules in a periodic cubic computational domain with box size 
$L = 12.8\,\textrm {nm}$. The mass density is 
$\rho ^{0} = 0.616\,\textrm {amu}$ Å
$^{-3}$, and the ion concentration is 
$c_{\pm } = 0.3168\,\textrm {M}$. The SPC/E model (Berendsen, Grigera & Straatsma Reference Berendsen, Grigera and Straatsma1987) is used for water molecules. The ionic force field terms are adopted from published work by Smith & Dang (Reference Smith and Dang1994) and Yoshida et al. (Reference Yoshida, Mizuno, Kinjo, Washizu and Barrat2014). The particle–particle particle-mesh (PPPM) method (Hockney & Eastwood Reference Hockney and Eastwood1988) is used for computing electrostatic interactions with vacuum periodic boundary conditions. The PPPM method maps atom charge to a three-dimensional mesh and uses fast Fourier transforms to solve Poisson's equation on the mesh, then interpolates electric fields on the mesh points back to the atoms, which is implemented in a built-in PPPM solver in LAMMPS (Thompson et al. Reference Thompson2022). A direct space cutoff of 9.8 Å, which is three times larger than the size of the oxygen atoms, is used to ensure the accuracy of the PPPM solver (Isele-Holder, Mitchell & Ismail Reference Isele-Holder, Mitchell and Ismail2012). The bond length of water molecules is constrained using the SHAKE algorithm (Ryckaert, Ciccotti & Berendsen Reference Ryckaert, Ciccotti and Berendsen1977) to allow a time step of 
$1\,\textrm {fs}$ in the velocity Verlet integrator (Allen & Tildesley Reference Allen and Tildesley2017). A constant number–volume–temperature system (NVT ensemble) is simulated at 
$T = 300\,\textrm {K}$ with a Nosé–Hoover thermostat. The system is relaxed for 
$0.5\,\textrm {ns}$ to achieve a thermal equilibrium state, and then for up to 
$10\,\textrm {ns}$ for statistics and sampling.
Figure 1. Typical snapshots of the aqueous NaCl solution from (a) all-atom MD simulations with explicit ions (box size 
$L=12.8\,\textrm {nm}$), and (b) mesoscale cDPD simulations with semi-implicit ions (box size 
$L=106.8\,\textrm {nm}$).
In § 2, we assumed that the mutual diffusion is small compared to ionic self-diffusion, thus it can be ignored in the derivation. To confirm this assumption, we compute both self-diffusion and mutual diffusion coefficients of ionic species from the MD velocity correlation functions via the Green–Kubo relationship (Zhou & Miller Reference Zhou and Miller1996; Wheeler & Newman Reference Wheeler and Newman2004). The computed self-diffusion coefficient is 
$1.3\times 10^{-9}\,\textrm {m}^{2}\,\textrm {s}^{-1}$ for Na
$^{+}$, and 
$2.1 \times 10^{-9}\,\textrm {m}^{2}\,\textrm {s}^{-1}$ for Cl
$^{-}$, while their mutual diffusion coefficient is 
$1.2 \times 10^{-10}\,\textrm {m}^{2}\,\textrm {s}^{-1}$. These results are consistent with previous simulation and experimental results for NaCl solutions (Wheeler & Newman Reference Wheeler and Newman2004), and also confirm that the mutual diffusion coefficient is much smaller than the self-diffusion coefficients in this electrolyte solution, thus it can be ignored in the analytical derivations of current correlation functions presented in § 2.
Alternative to all-atom MD simulation, at the mesoscopic scale, a cDPD model is used to tackle the challenge of simulating coupled fluctuating hydrodynamics and electrostatics with long-range Coulomb interactions. The cDPD model is an extension of the classical DPD model to solve numerically the fluctuating hydrodynamic and electrokinetic equations in the Lagrangian framework. In classical DPD models, ions can be represented by explicit charged particles with electrostatic interactions between explicit ions. Groot (Reference Groot2003) introduced a lattice to the DPD system to spread out the charges over the lattice nodes. Then the long-range portion of the interaction potential was calculated by solving the Poisson equation on the grid based on a PPPM algorithm (Hockney & Eastwood Reference Hockney and Eastwood1988) by transferring quantities (charges and forces) from the particles to the mesh, and vice versa. Because the mesh defines a coarse-grained length for electrostatic interactions, correlation effects on length scales shorter than the mesh size cannot be properly accounted for. Although the particle-to-mesh then mesh-to-particle mapping/redistribution can solve the Poisson equation for particle-based systems, its dependence on a grid may contradict the original motivation for using a Lagrangian method, and additional computational complexity and inefficiencies are introduced. To abandon grids and use a unifying Lagrangian description for mesoscopic electrokinetic phenomena, we developed the cDPD model (Deng et al. Reference Deng, Li, Borodin and Karniadakis2016) to solve the Poisson equation on moving cDPD particles rather than grids. More specifically, cDPD describes the solvent explicitly in a coarse-grained sense as cDPD particles, while the ion species are described semi-implicitly, i.e. using a Lagrangian description of ionic concentration fields, associated with each moving cDPD particle, as shown in figure 1(b), which provides a natural coupling between fluctuating electrostatics and hydrodynamics. The validation of the cDPD model has been confirmed by Deng et al. (Reference Deng, Li, Borodin and Karniadakis2016) in different problems, including electrostatic double-layer and electro-osmotic flows in micro-channels.
The state vector of a cDPD particle can be written as 
$(\boldsymbol {r}, \boldsymbol {v}, c_{\alpha }, \phi )$, which is characterized not only by its position 
$\boldsymbol {r}$ and velocity 
$\boldsymbol {v}$ as in the classical DPD model, but also by ionic species concentration 
$c_{\alpha }$ (with 
$\alpha$ representing the 
$\alpha$th ion type) and electrostatic potential 
$\phi$ on the particle. A cDPD system is simulated in a cubic periodic computational box with length 
$106.8\,\textrm {nm}$, where a cDPD particle is viewed as a coarse-grained fluid volume that contains the solvent and other charged species. Exchange of the concentration flux of charged species occurs between neighbouring cDPD particles, much like the momentum exchange in the classical DPD model. The governing equations of cDPD and model parameters are summarized in Appendix B. We will perform cDPD simulations of electrolyte solutions and compute current correlation functions of fluctuating electrokinetics to compare with theoretical predictions given in § 2.
We first examine the local fluctuations of mass density, momentum density, charge density and ionic concentration. For each of these quantities, we construct the instantaneous function 
$n(\boldsymbol {r})$ on a grid 
$\boldsymbol {r}^{i,j,k}$ by spatial averaging using a step function 
$W(r, h)$, with 
$h$ being the grid size, i.e. 
$W(r, h) = 1$ for 
$|\boldsymbol {r}_i-\boldsymbol {r}|\le 0.5h$, and 
$W(r, h) = 0$ for 
$|\boldsymbol {r}_i-\boldsymbol {r}| > 0.5h$. Then the local quantities can be extracted by 
$n(\boldsymbol {r}) = \sum _{i=1}^{N} W(|\boldsymbol {r} - \boldsymbol {r}_i|, h) n_i$ with 
$h = 3.1$ Å or 
$h = 7.75$ Å for MD data, and 
$h = 4.27\,\textrm {nm}$ for cDPD data. Because all instantaneous physical quantities are computed on a three-dimensional grid, the first point beyond 
$r = 0$ in spatial correlation functions (SCFs) presented in figures 3 and 4 is 
$r = h$, the second point is 
$r =\sqrt {2}h$, the third point is 
$r =\sqrt {3}h$, then 
$r =2h$, and so on. According to the central limit theorem, the local fluctuations in these quantities of interest should be Gaussian in the continuum regime. This is observed for mass density, momentum density and charge density from both MD and cDPD results. However, the local concentration fluctuations for ionic species in the MD system follow gamma distributions and converge to Gaussian with large spacing 
$h$, as shown in figure 2(a). Because the cDPD model assumes that the random fluxes of ionic concentration between neighbouring cDPD particles can be modelled by Gaussian white noises, the MD results indicate that the cDPD model is valid for length scales above a few nanometres, but has a lower bound for its length scale. Figure 2(b) presents the probability distribution functions of local ionic concentration in cDPD simulations, which follow a Gaussian distribution as expected.
Figure 2. Local ionic concentration probability distribution functions from (a) all-atom MD and (b) mesoscale cDPD simulations. The symbols show simulation data, while the lines show fits to gamma (MD) and Gaussian (inset of (a), and cDPD) distributions.
The SCFs for physical quantities of interest 
$n(\boldsymbol {r})$ are computed on the grid 
$\boldsymbol {r}^{i,j,k}$ according to
\begin{equation} {\rm SCF}(\boldsymbol{r}) = \frac{1}{N(\boldsymbol{r})} \sum_{i=1}^{N(\boldsymbol{r})} \delta n(\boldsymbol{r}_i)\,\delta n(\boldsymbol{r}_i + \boldsymbol{r}), \end{equation}
where 
$N(\boldsymbol {r})$ is defined as a normalization coefficient. We have assumed that the spatial correlations are isotropic in weakly charged electrolyte bulk solutions. We observe that the spatial correlations of mass density, momentum density, charge density and ion concentration are all very short ranged (approximately delta-correlated for large systems), which is in agreement with continuum fluctuating hydrodynamics results at equilibrium. However, as shown in figure 3, there is small-scale anti-correlation structure for the charge densities obtained in the MD simulations. We also compute the SCF of charge density in mesoscale cDPD simulations, and find similar small-scale anti-correlation structures for both positive and negative ions, as shown in figure 4. These anti-correlation structures are observed on length scales comparable to the Bjerrum and Debye lengths for electrolyte solution systems, suggesting that this could be related to the ionic screening effect.
Figure 3. Spatial correlation function (SCF) of charge density from full-atom MD simulations for grid distances (a) 
$0.31\,\textrm {nm}$ and (b) 
$0.775\,\textrm {nm}$.
Figure 4. SCF of charge density from mesoscale cDPD simulations.
The temporal correlation function (TCF) between two physical quantities 
$u$ and 
$w$ in Fourier space is defined by
\begin{equation} {\rm TCF}_k(t)=\langle u(\boldsymbol{k},t)\,w(\boldsymbol{k}, 0)\rangle = \frac{1}{N(t)} \sum_{s=1}^{N(t)} \hat{u}(\boldsymbol{k}, t)\,\hat{w}(\boldsymbol{k}, 0), \end{equation}
where 
$\boldsymbol {k}$ is the wave vector in Fourier mode, and 
$\hat {u}, \hat {w}$ are the Fourier components computed directly from particle trajectories as
\begin{equation} \hat{u}(\boldsymbol{k},t) = \frac{1}{N_P} \sum_{i=1}^{N_P} u_i(\boldsymbol{r}_i, t) \exp(-{\rm i} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{r}_i(t)). \end{equation}
We use a wavelength equal to the computational box size 
$L$ by setting 
$k=2{\rm \pi} /L$, which is 
$12.8\,\textrm {nm}$ for the MD system, and 
$106.79\,\textrm {nm}$ for the cDPD system.
Figures 5(a) and 5(b) present the temporal velocity autocorrelation functions from MD and cDPD simulations, respectively. The theoretical predication based on (2.14) gives 
$\textrm {CF}_v(t)=\exp (-0.1616t)$ (shear mode) and 
$\textrm {CF}_v(t)=\exp (-0.2403t)\cos (0.7352t)$ (sound mode) for the MD system, and 
$\textrm {CF}_v(t)=\exp (-0.0216t)$ (shear mode) and 
$\textrm {CF}_v(t)=\exp (-0.0238t)\cos (0.2177t)$ (sound mode) for the cDPD system. It can be observed from figure 5 that the computed results from both MD and DPD simulations are in agreement with the classical linearized theory of fluctuating hydrodynamics, i.e. the transverse autocorrelation function (shear mode) decays as 
$\exp (-\nu k^{2} t)$, with 
$\nu$ the kinematic shear viscosity, while the longitudinal autocorrelation function (sound mode) decays as 
$\exp (-\varGamma _T k^{2} t)\cos (c_s k t)$, with 
$\varGamma _T = 2\nu /3 + \zeta /2\rho$ the sound absorption coefficient and 
$c_s$ the isothermal sound speed. We conclude that the fluctuating hydrodynamics model still follows the classical linearized theory, and is not affected explicitly by the electrolyte bulk solutions. In an electrolyte solution in thermodynamic equilibrium, the fluctuating hydrodynamics and electrokinetics are decoupled explicitly.
Figure 5. Temporal transverse (shear mode) and longitudinal (sound mode) velocity autocorrelation functions in Fourier space from (a) MD and (b) cDPD simulations represented by open symbols, with comparison against the theoretical predictions in the form of (2.14) in solid and dashed lines.
Next, we compare the autocorrelation and cross-correlation functions of charge concentration with the linearized theory derived above, as shown in figure 6. We fit the MD and DPD results with the linearized theory in the form of (2.29) with proper structure factors 
$S_{p n}$ and 
$S_{n p }$, which are around 
$1.0$ and system-dependent, i.e.
\begin{equation} \left.\begin{gathered} {\varPsi}_{pp}(t)={\varPsi}_{nn}(t)={-}0.0490\exp(\lambda_1 t) + 1.0049\exp(\lambda_2 t), \\ {\varPsi}_{pn}(t)={\varPsi}_{np}(t)=0.0276\exp(\lambda_1 t) + 0.9724\exp(\lambda_2 t), \end{gathered}\right\} \end{equation}
with 
$\lambda _1 = -7.4187\times 10^{-3}$ and 
$\lambda _2 = - 3.9067\times 10^{-4}$ for the MD system, and
\begin{equation} \left.\begin{gathered} {\varPsi}_{pp}(t)={\varPsi}_{nn}(t)={-}0.2401\exp(\lambda_1 t) + 1.2401\exp(\lambda_2 t), \\ {\varPsi}_{pn}(t)={\varPsi}_{np}(t)=0.2458\exp(\lambda_1 t) + 0.7542\exp(\lambda_2 t), \end{gathered}\right\} \end{equation}
with 
$\lambda _1 = -3.4446\times 10^{-3}$ and 
$\lambda _2 = - 8.10\times 10^{-4}$ for the cDPD system. We observe in figure 6 that the computed results of temporal correlation functions from both MD and DPD simulations are in good agreement with the theoretical predictions derived from linearized theory of fluctuating hydrodynamics. Compared to the temporal correlation functions of hydrodynamic fluctuations presented in figure 5, the temporal correlations of charge concentration fluctuations shown in figure 6 decay extremely slowly over time. The long-time behaviour of both autocorrelation and cross-correlation follow the slow decay dominated by 
$\propto \exp (\lambda _2 t)$, which is significantly different from the decorrelation processes of hydrodynamic fluctuations.
Figure 6. Temporal cation and anion concentration autocorrelation (
$\varPsi _{pp}$ and 
$\varPsi _{nn}$) and cross-correlation (
$\varPsi _{pn}$ and 
$\varPsi _{np}$) functions in Fourier space from (a) MD and (b) cDPD simulations represented by open symbols, with comparison against the predictions from linearized fluctuating hydrodynamics theory in the form of (2.29) in solid and dashed lines.
4. Summary and discussion
We have presented theoretical closed-form expressions for the fluctuations of electrolyte bulk solution close to thermodynamic equilibrium, with an emphasis on mesoscopic spatiotemporal scales. In particular, we started with the Landau–Lifshitz theory and linearized the fluctuating hydrodynamic and electrokinetic equations to derive analytical solutions for current correlation functions using perturbation theory. To validate these theoretical expressions obtained based on the linearized theory, we performed numerical experiments of fluctuating hydrodynamics and electrokinetics for electrolyte solutions using both all-atom molecular dynamics (MD) and a mesoscale charged dissipative particle dynamics (cDPD) methods. We presented both MD and DPD simulation results of bulk electrolyte solutions in about 10 and 100 nm scales, which are compared directly with the predictions from the linearized continuum theory.
The current correlation functions computed from both MD and cDPD trajectories indicate that the temporal correlations of fluctuations from electrokinetics decay much more slowly than those from hydrodynamics, which agree well with the predictions made by the theoretical closed-form expressions of temporal correlation functions (hydrodynamics and charge) as well as the vast differences in decay time – varying by a few orders of magnitude. In a bulk electrolyte solution close to thermodynamic equilibrium, the fluctuating hydrodynamics and electrokinetics are decoupled explicitly because of the zero mean velocity, thus their behaviour is not dependent explicitly on the electrolyte concentrations in the dilute regime. At length scales above 10 nm, the results obtained from both MD and cDPD simulations are in good agreement with the continuum-limit linearized theories. The good agreement between the theoretical predictions and the particle-based mesoscopic simulations can also be interpreted as the capability of the mesoscopic cDPD model in bridging nano-to-continuum scales. Spatial correlations of charge density demonstrate finite range and non-trivial structure at nanometre length scales, but can be viewed as the delta function in the continuum limit. Simulation results also show that the fluctuations of local ionic concentration follow a gamma distribution at small length scales, while converging to a Gaussian distribution in the continuum limit, which suggests the existence of a lower bound of length scale for mesoscale models using Gaussian fluctuations for electrolyte solutions. Because the probability distribution functions of local ionic concentrations are computed directly from MD simulations before the continuum hypothesis is applied, a lower bound of length scale could exist in more general mesoscopic descriptions, including the Landau–Lifshitz PDE-based approach where stochastic fluxes are modelled as Gaussian processes.
It is worth noting that the present work focused on dilute electrolyte solutions close to thermodynamic equilibrium, where the mutual diffusion of ions is small compared to their self-diffusion, and the fluctuating electrokinetics is decoupled from zero-mean hydrodynamic fluxes. When the mutual diffusion coefficients of ions become comparable to the self-diffusion terms in concentrated electrolyte solutions (Dufrêche, Bernard & Turq Reference Dufrêche, Bernard and Turq2002; Galindres et al. Reference Galindres, Eslava, Ribeiro, Esteso, Vargas and Leaist2021), the contribution of mutual diffusion on ionic transport should be considered correctly in the fluctuating hydrodynamic and electrokinetic equations. Moreover, the linearized fluctuating hydrodynamics and electrokinetics system was solved in the 
$k$-space by a spatial Fourier transform, where the periodic condition is used in derivation of analytic solutions. Although both MD and cDPD simulations can be applied easily to electrokinetic problems involving fixed or induced charges on walls, the current theoretical framework cannot be extended to such problems directly. It is also interesting to consider how the fluctuating hydrodynamics is coupled with mesoscale electrokinetics in shear flows, where the current correction functions can be affected by the coupling between fluctuating terms and advection processes (Bian et al. Reference Bian, Deng and Karniadakis2018), leading to an orientation-dependent decorrelation process of fluctuating variables in the fluid system.
Funding
This work was supported by an ARO/MURI grant W911NF-15-1-0562 and the US Army Research Laboratory under Cooperative Agreement no. W911NF-12-2-0023. F.T. and Z.L. acknowledge research support from the NSF grant OAC-2103967 and the NASA project 80NSSC21M0153, and also Clemson University with generous allotment of computing time on the Palmetto cluster.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Derivation of mass–momentum correlations
The local hydrodynamic field can be expressed as the perturbations around the bulk state. We can define the instantaneous density, velocity and momentum as
$$\begin{gather} \rho(\boldsymbol{r},t) = \rho_0 + \delta \rho(\boldsymbol{r},t), \end{gather}$$
$$\begin{gather}\boldsymbol{u}(\boldsymbol{r},t) = \boldsymbol{u}_0 + \delta \boldsymbol{u}(\boldsymbol{r},t), \end{gather}$$
$$\begin{gather}\boldsymbol{g}(\boldsymbol{r},t) = \delta \boldsymbol{g} (\boldsymbol{r},t) = \rho\,\boldsymbol{u}(\boldsymbol{r},t). \end{gather}$$The linearized equations of fluctuating hydrodynamics are given by
$$\begin{gather} \frac{\partial [\delta \rho(\boldsymbol r,t)]}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} \delta \boldsymbol g(\boldsymbol r,t) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial [{\delta \boldsymbol g(\boldsymbol r,t)]}}{\partial t} + c_s^{2}\, \boldsymbol{\nabla} \delta \rho(\boldsymbol r,t) - \frac{\eta}{\rho_0}\, \nabla^{2} \delta \boldsymbol g(\boldsymbol r,t) - \frac{\dfrac{\eta}{3} + \zeta}{\rho_0}\,\boldsymbol{\nabla} (\boldsymbol{\nabla}\boldsymbol{\cdot} \delta\boldsymbol g(\boldsymbol r,t)) = 0. \end{gather}$$
The above equations can be transformed into 
$k$-space by a spatial Fourier transform using the function
\begin{equation} \hat{f_k}(t) = \int_V f(t) \exp(-{\rm i}\boldsymbol k \boldsymbol{\cdot} \boldsymbol r)\,{{\rm d}} \boldsymbol r. \end{equation}
To make the derivation simple, we set the wave vector 
$\boldsymbol k = (k, 0, 0)$ as one-dimensional along an arbitrary 
$x$-direction, leading to
$$\begin{gather} \frac{\partial \widehat{\delta \rho}_k(t)}{\partial t} + {\rm i} \boldsymbol k \boldsymbol{\cdot} \widehat{\delta \boldsymbol g}_k (t) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \widehat{\delta \boldsymbol g}_k(t)}{\partial t} + {\rm i}c_s^{2} \boldsymbol k\,\widehat{\delta \rho}_k(t) + \frac{\eta}{\rho_0}\,k^{2}\,\widehat{\delta \boldsymbol g}_k(t) + \frac{\dfrac{\eta}{3} + \zeta}{\rho_0}\,\boldsymbol k \boldsymbol k\,\widehat{\delta \boldsymbol g}_k(t) = 0. \end{gather}$$
We can divide the above equations into different components (
$x$-, 
$y$- and 
$z$-directions). The 
$x$-direction equations are
$$\begin{gather} \frac{\partial \widehat{\delta \rho}_k(t)}{\partial t} + {\rm i}k \cdot \widehat{\delta g}_k^{x} (t) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \widehat{\delta g}_k^{x}(t)}{\partial t} + {\rm i}c_s^{2} k\, \widehat{\delta \rho}_k(t) + \frac{\eta}{\rho_0}\,k^{2} \widehat{\delta g}_k^{x}(t) + \frac{\dfrac{\eta}{3} + \zeta}{\rho_0}\,k^{2}\,\widehat{\delta g}_k^{x}(t) = 0. \end{gather}$$
The 
$y$- and 
$z$-components are given by
$$\begin{gather} \frac{\partial \widehat{\delta g}_k^{y}(t)}{\partial t} + \frac{\eta}{\rho_0}\,k^{2}\,\widehat{\delta g}_k^{y}(t) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \widehat{\delta g}_k^{z}(t)}{\partial t} + \frac{\eta}{\rho_0}\, k^{2}\,\widehat{\delta g}_k^{z}(t) = 0. \end{gather}$$
Let 
$\nu = {\eta }/{\rho _0}$ and 
$\nu _L = (\frac {4}{3}\eta + \zeta )/\rho _0$. The above equations can be rewritten into the forms
$$\begin{gather} \frac{\partial \widehat{\delta \rho}_k(t)}{\partial t} + {\rm i}k \cdot \widehat{\delta g}_k^{x} (t) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \widehat{\delta g}_k^{x}(t)}{\partial t} + {\rm i}c_s^{2} k\,\widehat{\delta \rho}_k(t) + \nu_L k^{2}\,\widehat{\delta g}_k^{x}(t) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \widehat{\delta g}_k^{y}(t)}{\partial t} + \nu k^{2}\,\widehat{\delta g}_k^{y}(t) = 0, \end{gather}$$
$$\begin{gather}\frac{\partial \widehat{\delta g}_k^{z}(t)}{\partial t} + \nu k^{2}\,\widehat{\delta g}_k^{z}(t) = 0. \end{gather}$$The last two equations are first-order linear ODEs, which are solved easily as
$$\begin{gather} \widehat{\delta g}_k^{y}(t) = \widehat{\delta g}_k^{y}(0) \exp (-\nu k^{2} t), \end{gather}$$
$$\begin{gather}\widehat{\delta g}_k^{z}(t) = \widehat{\delta g}_k^{z}(0) \exp(-\nu k^{2} t). \end{gather}$$The first two equations of (A7) are two coupled ODEs and can be written in matrix form as
\begin{equation} \frac{{{\rm d}} \boldsymbol{a}_k(t)}{{{\rm d}}t} = \boldsymbol{\mathsf{H}}\,\boldsymbol{a}_k(t), \end{equation}where
\begin{equation}
\boldsymbol a_k(t) = \left[
\begin{array}{@{}c@{}} \widehat{\delta \rho}_k(t)\\
\widehat{\delta g}_k^{x}(t) \end{array}\right]\quad {\rm
and}\quad \boldsymbol{\mathsf{H}} = \left[ \begin{array}{@{}cc@{}} 0 &
-{\rm i}k \\ -{\rm i}c_s^{2}k & -\nu_Lk^{2} \end{array}
\right]. \end{equation}
The solutions of (A9) are determined by the eigenvalues of 
$\boldsymbol{\mathsf{H}}$, which can be obtained by solving the equation
\begin{equation} {\rm det} (\boldsymbol{\mathsf{H}} - \lambda \boldsymbol{\mathsf{I}}) = 0. \end{equation}The eigenvalues are given by
\begin{equation} \lambda_1 ={-}\varGamma_Tk^{2} + {\rm i}s_Tk\quad {\rm and}\quad \lambda_2 ={-}\varGamma_Tk^{2} - {\rm i}s_Tk, \end{equation}where
\begin{equation} \varGamma_T = \frac{\nu_L}{2} \quad {\rm and}\quad s_T = \frac{\sqrt{4c_s^{2} - \nu_L^{2}k^{2}}}{2}. \end{equation}
Here, 
$\varGamma _T$ is the sound absorption coefficient. Next, we consider an under-damped solution, where 
$s_T$ is real, that is, 
$k < 2c_s / \nu _L$. In particular, we consider the continuum limit, where 
$k \ll 2c_s / \nu _L$ so that 
$s_T \approx c_s$, then we arrive at the solutions (De Fabritiis et al. Reference De Fabritiis, Serrano, Delgado-Buscalioni and Coveney2007; Bian et al. Reference Bian, Deng and Karniadakis2018)
$$\begin{gather} \widehat{\delta \rho}_k(t) = {\rm e}^{-\varGamma_T k^{2} t}\left[ \cos(c_skt)\, \widehat{\delta \rho}_k(0) - \frac{{\rm i}}{c_s} \sin(c_skt)\,\widehat{\delta g}_k^{x}(0) \right], \end{gather}$$
$$\begin{gather}\widehat{\delta g}_k^{x}(t) = {\rm e}^{-\varGamma_T k^{2} t}\left[ \cos(c_skt)\, \widehat{\delta g}_k^{x}(0) - {\rm i}c_s \sin(c_skt)\,\widehat{\delta \rho}_k(0) \right]. \end{gather}$$
Therefore, the normalized temporal correlation functions of the mass–momentum fluctuations in 
$k$-space are given by
$$\begin{gather} \frac{\langle \widehat{\delta \rho}_k(t)\,\widehat{\delta \rho}_k(0) \rangle}{\langle \widehat{\delta \rho}_k(0)\,\widehat{\delta \rho}_{k}(0) \rangle} = {\rm e}^{-\varGamma_T k^{2} t} \cos (c_s kt), \end{gather}$$
$$\begin{gather}\frac{\langle \widehat{\delta g}_k^{x}(t)\,\widehat{\delta g}_k^{x}(0) \rangle}{\langle \widehat{\delta g}_k^{x}(0)\,\widehat{\delta g}_k^{x}(0) \rangle} = {\rm e}^{-\varGamma_T k^{2} t} \cos (c_s kt), \end{gather}$$
$$\begin{gather}\frac{\langle \widehat{\delta \rho}_k(t)\,{\rm i}\,\widehat{\delta g}_k^{x}(0) \rangle}{\langle \widehat{\delta \rho}_k(0)\,{\rm i}\,\widehat{\delta g}_k^{x}(0) \rangle} = {\rm e}^{-\varGamma_T k^{2} t} \sin (c_s kt), \end{gather}$$
where we have assumed that the initial cross-correlation is 
$\langle \widehat {\delta \rho }_k(0)\,\widehat {\delta g}_k^{x}(0) \rangle = 0$.
Appendix B. Charged dissipative particle dynamics (cDPD)
We consider a constant number–volume–temperature system (NVT ensemble) consisting of 
$N$ cDPD particles, with the state of each cDPD particle defined by its position 
$\boldsymbol {r}$, velocity 
$\boldsymbol {v}$, and ionic concentration 
$c_{\alpha }$ (with 
$\alpha$ representing the 
$\alpha$th ion species). The time evolution of the 
$i$th particle state with unit mass is governed by Newton's law and a transport equation (Deng et al. Reference Deng, Li, Borodin and Karniadakis2016):
$$\begin{gather} \frac{{{\rm d}}^{2} \boldsymbol{r}_i}{{{\rm d}}t^{2}} = \frac{{{\rm d}}\boldsymbol{v}_i}{{{\rm d}}t} = \boldsymbol{F}_i = \sum_{i \neq j}(\boldsymbol{F}^{C}_{ij}+\boldsymbol{F}^{D}_{ij}+\boldsymbol{F}^{R}_{ij}+\boldsymbol{F}^{E}_{ij}), \end{gather}$$
$$\begin{gather}\frac{{{\rm d}}c_{\alpha i}}{{{\rm d}}t} = q_{\alpha i} = \sum_{i \neq j}(q^{D}_{\alpha ij} + q^{E}_{\alpha ij} + q^{R}_{\alpha ij}), \end{gather}$$
where 
$\boldsymbol {F}_i$ denotes the total force exerted on the 
$i$th particle, which consists of the conservative, dissipative and random forces. Additionally, the electrostatic force 
$\boldsymbol {F}_i^{E}$ is introduced to couple the hydrodynamics and electrokinetics within the DPD framework. In particular,
$$\begin{gather} \boldsymbol{F}^{C}_{ij} = a_{ij}\,\omega_C(r_{ij})\,\hat{\boldsymbol{r}}_{ij}, \end{gather}$$
$$\begin{gather}\boldsymbol{F}^{D}_{ij} ={-} \gamma_{ij}\, \omega_D(r_{ij})\,(\hat{\boldsymbol{r}}_{ij} \boldsymbol{\cdot}\boldsymbol{v}_{ij})\,\hat{\boldsymbol{r}}_{ij}, \end{gather}$$
$$\begin{gather}\boldsymbol{F}^{R}_{ij} = \sigma_{ij}\,\omega_R(r_{ij})\,\theta_{ij}\,\delta t^{{-}1/2}\,\hat{\boldsymbol{r}}_{ij}, \end{gather}$$
$$\begin{gather}\boldsymbol{F}^{E}_{ij} = \lambda_{ij} \left(\sum_{\alpha=1}^{S} z_\alpha c_{\alpha i}\right)\boldsymbol{E}_{ij}, \end{gather}$$
where 
$r_{ij}=|\boldsymbol {r}_{ij}|=|\boldsymbol {r}_i-\boldsymbol {r}_j|$, 
$\hat {\boldsymbol {r}}_{ij}=\boldsymbol {r}_{ij}/r_{ij}$ and 
$\boldsymbol {v}_{ij}=\boldsymbol {v}_i-\boldsymbol {v}_j$. The conservative, dissipative and random forces are pairwise forces with weighting functions 
$\omega _{C}(r_{ij})$, 
$\omega _D(r_{ij})$, 
$\omega _R(r_{ij})$, and corresponding strengths 
$a_{ij}$, 
$\gamma _{ij}$, 
$\sigma _{ij}$, respectively. The 
$\theta _{ij}$ are symmetric Gaussian random variables with zero means and unit variances; these variables are independent for different pairs of particles at different times; 
$\theta _{ij} = \theta _{ji}$ is enforced to satisfy momentum conservation. The dissipative and random forces together act as a thermostat, with their coefficients and weighting functions satisfying the fluctuation–dissipation theorem (FDT) (Español & Warren Reference Español and Warren1995)
\begin{equation} \sigma_{ij}^{2}=2k_BT\gamma_{ij},\quad \omega_D(r_{ij})=\omega_R^{2}(r_{ij}), \end{equation}
where 
$k_B$ is the Boltzmann constant and 
$T$ is the temperature. The coupling parameter 
$\lambda _{ij}$ in electrostatic force is introduced by rescaling the PNP equations with DPD units, which are related linearly to the macroscopic dimensionless coupling parameter 
$\varLambda = c_0^{*} \cdot k_BT \tau ^{2}/({\rho _0 r_0^{5}})$ with 
$c_0^{*} = c_0 r_0^{3}$, the reference concentration in DPD units (usually 
$c_0$ as bulk concentration, 
$r_0$ as unit DPD length). Here, 
$\boldsymbol {E}_{ij}$ is the relative electric fields difference between particles 
$i$ and 
$j$, which is determined by the electrostatic potential field 
$\phi$:
\begin{equation} \boldsymbol{E}_{ij} = (\phi_i - \phi_j)\,\omega_E(r_{ij})\,\boldsymbol{r}_{ij}, \end{equation}
with 
${\omega }_E(r)$ a weighting function. It is very important to note that the electrostatic forces here are not pairwise additive, i.e. 
$\boldsymbol {F}^{E}_{ij} \neq \boldsymbol {F}^{E}_{ji}$; however, we have 
$\sum _{i,j} \boldsymbol {F}^{E}_{ij} = 0$ to guarantee the global momentum conservation when there is no external electrostatic field. In the cDPD framework, the ionic concentration (rescaled by the reference concentration 
$c_0$) evolution is driven by three pairwise flux terms, i.e. the Fickian flux 
$q_{\alpha ij}^{D}$, electrostatic flux 
$q_{\alpha ij}^{E}$ and random flux 
$q_{\alpha ij}^{R}$, induced by concentration gradient, electrostatic potential gradient and thermal fluctuations, respectively. Specifically,
$$\begin{gather} q_{\alpha ij}^{D} ={-}\kappa_{\alpha ij}(c_{\alpha i}-c_{\alpha j})\,\omega_{DC}(r_{ij}), \end{gather}$$
$$\begin{gather}q_{\alpha ij}^{E} ={-}\tfrac{1}{2}\kappa_{\alpha ij}z_\alpha (c_{\alpha i} + c_{\alpha j})(\phi_i - \phi_j)\,\omega_{DC}(r_{ij}), \end{gather}$$
$$\begin{gather}q_{\alpha ij}^{R} = \xi_{\alpha ij}\,\omega_{RC} (r_{ij})\,\theta_{ij}\, \delta t^{{-}1/2}, \end{gather}$$
where 
$\kappa _{ij}$ are the diffusion coefficients, and 
$\omega _{DC}(r_{ij})$ and 
$\omega _{RC}$ are weighting functions. The coefficients and weighting functions of the random flux are determined via the generalized FDT as
\begin{equation} {\xi_{\alpha ij}}^{2} = \frac{\kappa_{\alpha ij}}{c_0^{*}}\,(c_{\alpha i}+c_{\alpha j}),\quad \omega_{DC}(r)=\omega_{RC}^{2}(r). \end{equation}
The electrostatic potential 
$\phi$ on each cDPD particle is determined by solving a modified Poisson equation at every DPD time step. In cDPD, we consider the dimensionless modified Poisson equation rescaled by DPD units,
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\epsilon (\boldsymbol{r})\,\boldsymbol{\nabla} (\phi(\boldsymbol{r})\right) ={-}\varGamma\, \rho_e(\boldsymbol{r}), \end{equation}
with 
$\varGamma =e^{2}c_0^{*}r_0^{2}/\epsilon _0 k_BT$, and 
$\rho _e = \sum _{\alpha =1}^{S} z_\alpha c_\alpha$ is the charge density. The electrostatic potential 
$\phi _i$ on the 
$i$th particle is obtained together via a successive over-relaxation iteration scheme as
\begin{equation} \phi_i^{k} = \phi_i^{k-1} + \vartheta\left[\sum_{\alpha=1}^{S} \varGamma z_\alpha c_\alpha - \sum_{j \neq i}\bar{\epsilon}_{ij} \phi_{ij}^{k}\,\omega_{\phi}(r_{ij}) \right], \end{equation}
where 
$\omega _{\phi }(r)$ is the weight function, 
$k$ represents an iteration step, 
$\vartheta$ is the relaxation factor, and 
$\epsilon _i$ and 
$\epsilon _j$ are the permittivities of the 
$i$th and 
$j$th cDPD particles, respectively. Here, 
$\bar {\epsilon }_{ij}=(\epsilon _i + \epsilon _j)/2$, where 
$\epsilon _i$ and 
$\epsilon _j$ can be different values to model mixtures of heterogeneous solvents.
The initial guesses for 
$\phi _i^{k-1}$ take the value of 
$\phi _i$ from the previous time step. The iteration stops when the absolute differences 
$|\phi _i^{k}-\phi _i^{k-1}|$ are smaller than a tolerance, i.e. 
$10^{-3}$ for all cDPD particles. The relaxation factor 
$\vartheta$ is selected adaptively during the iteration to optimize for faster convergence. For more details on the derivations leading to (B1)–(B8), we refer interested readers to the previous work Deng et al. (Reference Deng, Li, Borodin and Karniadakis2016).
Throughout this paper, the weighting functions are chosen as 
$\omega _C(r) = (1-r/r_c)$, 
$\omega _D(r) = \omega _R^{2}(r) = (1-r/r_c)$, 
$\omega _{qD}(r) = \omega _{qR}^{2}(r) = (1-r/r_{cc})^{2}$, 
$\omega _\phi (r) = (1-r/r_{ec})^{2}$ and 
$\omega _E(r) = 0.5(1-r/r_{ec})^{2} r$. The cDPD parameters used for the simulation are summarized in table 1. The basic DPD units according to these parameters are 
$r_0=21.36$ nm for the length unit, 
$\tau =3.28$ ns for the time unit, 
$k_BT=4.14\times 10^{-21}$ J for the energy units, and 
$c_0=4.08\times 10^{-3}$ M for the concentration unit.
Table 1. Parameter list for cDPD simulations. The values are provided in reduced DPD units. Length unit 
$r_0=21.36\,\textrm {nm}$, time unit 
$\tau =3.28\,\textrm {ns}$, energy unit 
$k_BT=4.14\times 10^{-21}\,\textrm {J}$ and concentration unit 
$c_0=4.08\,\textrm {mM}$ are used for mapping to physical units.
Open access
