Transverse Density Fluctuations around the Ground State Distribution of Counterions near One Charged Plate: Stochastic Density Functional View

Abstract

We consider the Dean–Kawasaki (DK) equation of overdamped Brownian particles that forms the basis of the stochastic density functional theory. Recently, the linearized DK equation has successfully reproduced the full Onsager theory of symmetric electrolyte conductivity. In this paper, the linear DK equation is applied to investigate density fluctuations around the ground state distribution of strongly coupled counterions near a charged plate, focusing especially on the transverse dynamics along the plate surface. Consequently, we find a crossover scale above which the transverse density dynamics appears frozen and below which diffusive behavior of counterions can be observed on the charged plate. The linear DK equation provides a characteristic length of the dynamical crossover that is similar to the Wigner–Seitz radius used in equilibrium theory for the 2D one-component plasma, which is our main result. Incidentally, general representations of longitudinal dynamics vertical to the plate further suggest the existence of advective and electrical reverse-flows; these effects remain to be quantitatively investigated.

1. Introduction

Water-soluble materials often have surface chemical groups that are dissociated in a polar solvent. Examples of such materials include not only mesoscopic particles, such as viruses, proteins, polyelectrolytes, membranes, and micelles, but also macroscopic objects like a glass plate of sample cell [1,2]. Both mesoscopic and macroscopic particles will be referred to here as ‘macroions’. The macroions are likely to carry a total surface charge exceeding thousands of elementary charges e, surrounded by oppositely charged counterions that are dissociated from the macroions [1,2]; counterions are electrostatically bound around macroions, due to the high asymmetry between counterions and macroions in valence of charges.

Focusing on the counterions, the macroion systems can be rephrased as inhomogeneous one-component ionic fluids in the presence of external fields—the one-component fluids of counterions can be regarded as the one-component plasma (OCP) [3,4,5]. Systems of charged particles immersed in a smooth neutralizing medium are commonly observed in nature, such as a suspension of dust grains in plasmas, as well as colloidal solutions, which can be modeled by the OCP in the unscreened limit of Yukawa fluids [3,4,5]. The 2D OCP has been used for the description of dusty plasmas confined by external fields, where the motion of particles interacting via 3D electrostatic interaction potential is restricted to a 2D surface [3,4,5]. There is, however, a crucial difference between counterion systems and the OCP, due to the localization of the electrically neutralizing background. While the whole space in the OCP is filled with a smooth background in either the 2D or 3D systems, counterions form a 3D electric double layer and are not neutralized unless they are localized on the macroion surfaces [6,7,8].

This paper will address the strong coupling systems of counterions in the presence of one charged plate, focusing especially on the transverse dynamics of density fluctuations around the ground state that will be specified in Section 2.1 [8,9,10,11]. Turning our attention to the dynamics, the strongly coupled counterion system is distinguished from the 2D OCP by one extra dimension, vertical to the macroion surface. Accordingly, density fluctuations occur not only along the 2D plane parallel to the macroion surface, but also along the one extra dimension. Furthermore, the strong coupling regime in the counterion systems may be realized at room temperature, and therefore dynamics due to thermal fluctuations need to be considered; however, there are few studies on the counterion dynamics in the ground state.

Thus, the main aim of this paper is to investigate the anisotropic fluctuation field $n(\mathbf{r},t)$ of counterion density due to the coarse-grained dynamics of counterion density $\rho (\mathbf{r},t)={\rho }_{\infty }(\mathbf{r})+n(\mathbf{r},t)$ around a ground state distribution ${\rho }_{\infty }(\mathbf{r})$, using the stochastic density functional equation. The stochastic density functional theory is based on the so-called Dean–Kawasaki (DK) equation that describes the evolution of the instantaneous microscopic density field of overdamped Brownian particles [12,13,14,15,16,17,18,19,20]. The stochastic density functional theory has been used as one of the most powerful tools for describing slowly fluctuating and/or intermittent phenomena [16,17,18,19,20], such as glassy dynamics, nucleation or pattern formation of colloidal particles, stochastic thermodynamics of colloidal suspensions, dielectric relaxation of Brownian dipoles, and even tumor growth.

The original DK equation includes nonlinear terms of dynamic origin due to the kinetic coefficient that is proportional to fluctuating field $\rho (\mathbf{r},t)$ [12,13,14,15,16,17,18,19,20]. While the nonlinearility of the original DK equation leads to the above successful descriptions of various phenomena [16,17,18,19,20], a more tractable form is required. It has been recently demonstrated that the DK equation can be linearized with respect to $n(\mathbf{r},t)$ when $n(\mathbf{r},t)/{\rho }_{\infty }(\mathbf{r})\ll 1$, and that the linear stochastic equation of density fluctuations is of great practical use [21,22,23,24,25,26,27]. The density fluctuations of fluids near equilibrium are surprisingly well described by model-B dynamics of a Gaussian field theory whose effective quadratic Hamiltonian for the density fluctuation field is constructed to yield the exact form of the static density-density correlation function [25]. Furthermore, we have demonstrated that the DK equation can be directly linearized in the first approximation of the driving force due to the free energy functional $F\left[\rho \right]$ of an instantaneous density distribution $\rho$, when small density fluctuations around a metastable state are considered [21].

The stochastic thermodynamics around a metastable state has been investigated using the stochastic density functional equation (the DK equation), showing that the heat dissipated into the reservoir is generally negligible [21]. The linear stochastic density functional theory has also been found relevant to investigate out-of-equilibrium phenomena, including the formulations of the full Onsager theory of electrolyte conductivity [22].

The remainder of this paper is organized as follows: Section 2 provides formal background in the case of a single charged plate system. We give the linear DK equation as a stochastic density functional equation, after specifying a general form of the free energy functional $F\left[\rho \right]$ of a given density. In Section 3, the linear DK equation is applied to the strongly coupled counterion system by considering density fluctuations $n=\rho -{\rho }_{\infty }$ around ${\rho }_{\infty }$. We can verify that the first derivative of $F\left[\rho \right]$ in the ground state (i.e., ${\delta F\left[\rho \right]/\delta \rho |}_{\rho ={\rho }_{\infty }}$) produces a constant, similar to the above metastable state. Accordingly, the DK equation of the counterion system can be linearized around the ground state. We will also see the underlying physics of anisotropic fluctuations (vertical to the plate) in terms of the general form of the linear DK equation. In Section 4, we focus on the transverse dynamics along the plate surface, assuming the absence of the gradient of the fluctuating density field vertical to the plate. First, we derive the frozen dynamics over a long-range scale beyond the Wigner–Seitz cell, reflecting the formation of the Wigner crystal on the charged plate. Furthermore, the linear DK equation determines a crossover length $l_c$, below which we can observe diffusive behaviors of counterions condensed on the plate. Our main result in this study is the quantitative evaluation of $l_c$, yielding $l_c\sim a$ for $\Xi \sim {10}^3$. Section 5 contains a summary and conclusions.

2. Formal Background

2.1. Ground State of Counterion System in the Strong Coupling Limit

Let us briefly summarize what has been achieved by theoretical and simulation studies on the OCP and the counterion systems in the strong coupling regime.

The thermodynamics of the OCP system is characterized by the coupling parameter [3,4,5],

$$\Gamma =\frac{q^2{l}_B}{a},$$

(1)

where q is the valence of counterions, $l_B\equiv e^2/4\pi ϵk_{B}T$ is the distance (the so-called Bjerrum length) at which two elementary charges interact electrostatically with thermal energy $k_{B}T$, when they are surrounded by a polar solvent with its dielectric permittivity and temperature being $ϵ$ and T, and the Wigner–Seitz cell radius a defined by $(\pi a^2)\sigma =q$ using the surface number density $\sigma$ of macroion charges [3,4,5,7]. Thermodynamic properties of OCP systems have been extensively studied over decades and accurate numerical results as well as their fits are available in the literature [3,4,5]. As $\Gamma$ increases, the OCP shows a transition from a weakly coupled gaseous regime ($\Gamma \ll 1$) to a strongly coupled fluid regime ($\Gamma \gg 1$), and it eventually crystallizes. The concept of the Wigner crystallization due to long-range electrostatic interactions underlies the formation of colloidal crystals, or photonic crystals with large lattice constant, comparable in magnitude to the wavelength of visible light [2,3,4,5].

Meanwhile, in counterion systems, a Wigner–Seitz radius a has not been adopted in rescaling the Bjerrum length as $l_B/a$. We have used another coupling parameter $\Xi$ defined by [6,7,8,9,10,11]

$$\begin{array}{cc}\hfill & \Xi =\frac{q^2{l}_B}{\lambda },\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & \lambda =\frac{1}{2\pi q{l}_B\sigma },\hfill \end{array}$$

(3)

using the Gouy–Chapman length $\lambda =1/(2\pi q{l}_B\sigma )$, a characteristic length of the electric double layer. Inserting the relation $\sigma =q/(\pi a^2)$ into the definition of $\lambda$, we have

$$\frac{a}{\lambda }=2\pi q\sigma l_{B}a=\frac{2q^2{l}_B}{a}=2\Gamma .$$

(4)

It follows from Equations (3) and (4) that

$$\Xi =2{\Gamma }^2.$$

(5)

When macroion-counterion attractions are weak, the structure of such ionic cloud, characterized by $\lambda$, will be dispersed instead of forming the 2D OCP. The dispersed electric double layer is thus represented by the weak coupling parameter of $\Xi \ll 1$, where the Poisson–Boltzmann approach and its systematic improvements via the loop expansion has been found relevant [6,7,8]. On the other hand, as $\Xi$ increases while reducing $\lambda$, the electric double layer thins and the coarse-grained distribution of counterions becomes two-dimensional [6,7,8,9,10,11]. Correspondingly, a becomes much larger than $\lambda$ in the strong coupling regime of $\Xi =2{\Gamma }^2\gg 1$ (i.e., $a\gg \lambda$), as found from Equations (4) and (5).

A field theoretical treatment provides counterion density distribution, ${\rho }_{\infty }(z_0)$, in the ground state of the strong coupling limit ($\Xi ,\Gamma \to \infty$) [6,8]:

$$\begin{array}{c}\hfill {\rho }_{\infty }(z_0)=\frac{\sigma }{\lambda }exp\{-J({\mathbf{r}}_0)\},\end{array}$$

(6)

where $J({\mathbf{r}}_0)$ denotes the external electrostatic potential in the $k_{B}T$-unit due to macroion-counterion interactions. In a single charged plate system, $J({\mathbf{r}}_0)$ is expressed as $J({\mathbf{r}}_0)=z_0/\lambda$ with $z_0$ denoting the distance between the position ${\mathbf{r}}_0$ and the charged plate; therefore, Equation (6) implies that a large portion of the counterions are condensed within the thin electric double layer, which supports the observation that strongly coupled counterions behave like the 2D OCP. Extensive Monte Carlo simulations have been performed on the strong coupling regimes of counterions, especially for one- and two-plate systems [8,9,10,11]. Accordingly, the asymptotic behavior given by Equation (6) has been corroborated by simulation results on the counterion distributions. The correction to Equation (6) has also been evaluated in detail, based on the simulation results. For about two decades, a variety of strong coupling theories have been developed to explain the above simulation results in the strong coupling regime, focusing not only on the validation of the longitudinal distribution mimicked by the ideal gas behavior (i.e., Equation (6)) in the vicinity of the charged plate, but also on the deviations from Equation (6) for $z_0>\lambda$; see [8] for a recent review.

2.2. Imposing a Given Density Distribution $\rho$ on the Grand Potential $\Omega$

Let $\widehat{\rho }(\mathbf{r})$ be an instantaneous density of counterions located at ${\mathbf{r}}_i(i=1,\cdots ,N)$, where the counterion system is rescaled as $\mathbf{r}=(x,y,z)=(x_0/a,y_0/a,z_0/a)={\mathbf{r}}_0/a$. Figure 1 shows a schematic of the rescaled system. The instantaneous density is expressed as

$$\widehat{\rho }(\mathbf{r})=a^3\sum _{i=1}^N\delta (\mathbf{r}-{\mathbf{r}}_i),$$

(7)

the use of which macroion-counterion interaction energy for a one-plate system transforms the configurational representation, given by Equation (A7) of Appendix B, to

$$\begin{array}{cc}\hfill \Delta U_{cm}\left\{\widehat{\rho }\right\}& =\int d\mathbf{r}J(\mathbf{r})\widehat{\rho }(\mathbf{r}),\hfill \\ \hfill J(\mathbf{r})& =2\Gamma z,\hfill \end{array}$$

(8)

which is a functional of $\widehat{\rho }$ (see Appendix B for the details).

We can impose a given density distribution $\rho (\mathbf{r})$ on the counterion system via the following delta functional [14,19,20,21]:

$$\prod _{\mathbf{r}}\delta \left[\widehat{\rho }(\mathbf{r})-\rho (\mathbf{r})\right]=\int D\psi e^{\int d\mathbf{r}i\psi (\mathbf{r})\{\widehat{\rho }(\mathbf{r})-\rho (\mathbf{r})\}},$$

(9)

where the potential field $\psi (\mathbf{r})$ has been introduced in the Fourier transform of the delta functional. Multiplying the configurational integral representation of the grand potential $\Omega \left[J\right]$ (see Appendix B for the definition) by the constraint in Equation (9), the formal expression of $F\left[\rho \right]$ is obtained:

$$\begin{array}{cc}{e}^{-F\left[\rho \right]}\hfill & ={\displaystyle \prod _{\mathbf{r}}}\delta \left[\widehat{\rho }(\mathbf{r})-\rho (\mathbf{r})\right]e^{-\Omega \left[J\right]}\hfill \\ \hfill & =e^{-\int d\mathbf{r}J(\mathbf{r})\rho (\mathbf{r})}\int D\psi e^{-\Omega [-i\psi ]-\int d\mathbf{r}i\psi (\mathbf{r})\rho (\mathbf{r})},\hfill \end{array}$$

(10)

where $J(\mathbf{r})$ given in Equation (8) represents the external potential created by the charged plate. In the mean-field approximation, we obtain (see also Appendix C)

$$\begin{array}{cc}\hfill F\left[\rho \right]& =\mathcal{A}[\rho ]+\int d\mathbf{r}J(\mathbf{r})\rho (\mathbf{r}),\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \mathcal{A}\left[\rho \right]& =\Omega \left[{\psi }^{*}\right]-\int d\mathbf{r}{\psi }^{*}(\mathbf{r})\rho (\mathbf{r}),\hfill \end{array}$$

(12)

where the saddle-point potential field ${\psi }^{*}$ satisfies the following relation:

$${\left.\frac{\delta \left(\beta \Omega [-i\psi ]\right)}{\delta \psi (\mathbf{r})}\right|}_{\psi =i{\psi }^{*}}=-i\rho (\mathbf{r}).$$

(13)

The first Legendre transform of $\Omega \left[{\psi }^{*}\right]$ provides the Hohenberg–Kohn free energy $\mathcal{A}\left[\rho \right]$ defined by Equation (12), the central functional of the equilibrium density functional theory [28,29].

2.3. Stochastic Density Dynamics Obeying the Dean–Kawasaki Equation

Here we focus on the stochastic dynamics of a density field at time t, $\rho (\mathbf{r},t)$, whose spatially varying distribution is the same as the coarse-grained variation of $\widehat{\rho }(\mathbf{r})$. What matters in terms of the stochastic density dynamics is the free energy functional $F\left[\rho \right]$ of a given density field $\rho$, rather than the grand potential $\Omega$ in equilibrium. For the density functional, we have provided the approximate form of $F\left[\rho \right]$. The driving force due to $F\left[\rho \right]$ and the density-dependent multiplicative noise $\zeta [\rho ,\overrightarrow{\eta }]$ creates the stochastic dynamics that obeys the DK equation [13,14,15,16,17,18,19,20]:

$${\partial }_t\rho =\nabla ·D\rho \nabla \frac{\delta F\left[\rho \right]}{\delta \rho }+\zeta [\rho ,\overrightarrow{\eta }],$$

(14)

where we have introduced a scaled diffusion constant, $D=D_0/a^2$, using the bare diffusion constant $D_0$ and the spatio-temporal average of the multiplicative noise correlations is given by

$$\begin{array}{c}〈\zeta [\rho (\mathbf{r},t),\overrightarrow{\eta }(\mathbf{r},t)]\zeta [\rho ({\mathbf{r}}^{\prime },t^{\prime }),\overrightarrow{\eta }({\mathbf{r}}^{\prime },t^{\prime })]〉=-2D\delta (t-t^{\prime }){\nabla }_{\mathbf{r}}·\rho (\mathbf{r},t){\nabla }_{\mathbf{r}}\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),\hfill \end{array}$$

(15)

with the vectorial white noise field $\overrightarrow{\eta }(\mathbf{r},t)$ that has the correlation $〈{\eta }_l(\mathbf{r},t){\eta }_m({\mathbf{r}}^{\prime },t^{\prime })〉={\delta }_{lm}\delta (\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-t^{\prime })$. Equation (15) can read [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27]

$$\zeta [\rho ,\eta ]=\nabla ·\sqrt{2D\rho (\mathbf{r},t)}\overrightarrow{\eta }(\mathbf{r},t).$$

(16)

In general, the stochastic Equation (14) includes not only the multiplicative noise term, but also the nonlinear term associated with $F\left[\rho \right]$. As shown below, however, a linear DK equation may be used to investigate the stochastic density dynamics due to fluctuations of counterions in the ground state by expanding Equation (14) around ${\rho }_{\infty }$.

3. Stochastic Density Functional Equation for Fluctuations round the Ground State Distribution ${\rho }_{\infty }$

3.1. Linearizing the Stochastic Dean–Kawasaki Equation (14)

Expanding the first derivative of $F\left[\rho \right]$ around ${\rho }_{\infty }$, we have

$$\begin{array}{cc}\nabla \frac{\delta F\left[\rho \right]}{\delta \rho (\mathbf{r},t)}\hfill & =\nabla \left[{\left.\frac{\delta F\left[\rho \right]}{\delta \rho (\mathbf{r},t)}\right|}_{\rho ={\rho }_{\infty }}+\int d{\mathbf{r}}^{\prime }{\left.\frac{{\delta }^{2}F\left[\rho \right]}{\delta \rho (\mathbf{r},t)\delta \rho ({\mathbf{r}}^{\prime },t)}\right|}_{\rho ={\rho }_{\infty }}n({\mathbf{r}}^{\prime })\right]\hfill \\ \hfill & =\nabla \left[-\int d{\mathbf{r}}^{\prime }c(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime },t)+\frac{n(\mathbf{r},t)}{{\rho }_{\infty }(z)}\right],\hfill \end{array}$$

(17)

$$\begin{array}{c}\rho (\mathbf{r},t)={\rho }_{\infty }(\mathbf{r})+n(\mathbf{r},t),\end{array}$$

(18)

due to

$$\nabla ·D\rho (\mathbf{r}){\left.\nabla \frac{F\left[\rho \right]}{\delta \rho }\right|}_{\rho ={\rho }_{\infty }}=0;$$

(19)

See Appendix D for the details. It follows from Equations (17) and (19) that the right-hand side (rhs) of Equation (14) reads

$$\begin{array}{cc}\nabla ·D\rho \nabla \frac{\delta F\left[\rho \right]}{\delta \rho (\mathbf{r},t)}+\zeta [\rho ,\overrightarrow{\eta }]\hfill & =\nabla ·D{\rho }_{\infty }\left(1+\frac{n(\mathbf{r},t)}{{\rho }_{\infty }}\right)\nabla \left\{-\int d{\mathbf{r}}^{\prime }c(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime },t)+\frac{n(\mathbf{r},t)}{{\rho }_{\infty }(z)}\right\}+\zeta \left[{\rho }_{\infty }\left(1+\frac{n(\mathbf{r},t)}{{\rho }_{\infty }}\right),\overrightarrow{\eta }\right]\hfill \\ \hfill & \approx \nabla ·D{\rho }_{\infty }\nabla \left\{-\int d{\mathbf{r}}^{\prime }c(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime },t)+\frac{n(\mathbf{r},t)}{{\rho }_{\infty }(z)}\right\}+\zeta [{\rho }_{\infty },\overrightarrow{\eta }],\hfill \end{array}$$

(20)

when $n/{\rho }_{\infty }\ll 1$. It is to be noted that the ideal gas distribution ${\rho }_{\infty }(\mathbf{r})$ given by Equation (6) reproduces only the longitudinal distribution of simulation results in the vicinity of a highly charged plate, resulting from attractive and repulsive Coulomb interactions in the strong coupling limit (see also Appendix A). Correspondingly, Equation (20) reveals that the dynamics of fluctuating density $n(\mathbf{r},t)$ is governed by the strong Coulomb interactions as represented by the contribution, $-\int d{\mathbf{r}}^{\prime }c(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime },t)$, on the rhs of Equation (20). In Equations (17) and (20), the direct correlation function $c(\mathbf{r}-{\mathbf{r}}^{\prime })$ appears because $F\left[\rho \right]$ is expressed using the Hohenberg–Kohn free energy functional $\mathcal{A}\left[\rho \right]$ as given by Equation (12) [28,29]. We also have

$${\partial }_t\{{\rho }_{\infty }(z)+n(\mathbf{r},t)\}={\partial }_{t}n(\mathbf{r},t),$$

(21)

due to ${\partial }_t{\rho }_{\infty }=0$, on the left hand side of the DK Equation (14).

Combining Equations (14), (20), and (21), we obtain the linear DK equation:

$$\begin{array}{cc}\hfill {\partial }_{t}n(\mathbf{r},t)& =D{\nabla }^{2}n(\mathbf{r},t)+D{\rho }_{\infty }(z){\nabla }^2{\psi }_n(\mathbf{r},t)-\nabla ·{\mathbf{j}}_{\perp }(\mathbf{r},t)+\zeta [{\rho }_{\infty },\overrightarrow{\eta }],\hfill \end{array}$$

(22)

where ${\psi }_n$ denotes a fluctuating Coulomb potential defined by

$${\psi }_n(\mathbf{r},t)\equiv -\int d{\mathbf{r}}^{\prime }c(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime }),$$

(23)

and the longitudinal current ${\mathbf{j}}_{\perp }(\mathbf{r},t)=(0,0,j_z)$, which is along the z-axis vertical to the charged plate, arises from the gradient of the ground density distribution $\nabla {\rho }_{\infty }$. Incidentally, in the rescaled system of ${\rho }_{\infty }(\mathbf{r},t)=a^3{\rho }_{\infty }({\mathbf{r}}_0,t)$, Equation (6) is rewritten as

$${\rho }_{\infty }(z)=\frac{\sigma a^3}{\lambda }{e}^{-2\Gamma z},$$

(24)

thereby providing

$$\nabla {\rho }_{\infty }=(0,0,-2\Gamma {\rho }_{\infty }),$$

(25)

in the rescaled system. In the next subsection, we will discuss the details of ${\mathbf{j}}_{\perp }(\mathbf{r},t)$ associated with the presence of $\nabla {\rho }_{\infty }$ in the longitudinal direction.

Equation (22) appears to be a simple extension of the Poisson–Nernst–Planck equation [30] when the longitudinal current ${\mathbf{j}}_{\perp }(\mathbf{r},t)$ disappears. However, in actuality, the Poisson-like equation in the second term on the rhs of Equation (22) differs from the conventional Poisson equation because the interaction potential is replaced by the direct correlation function. In this paper, we adopt the direct correlation function, being of the following form [31,32,33]:

$$\begin{array}{cc}\hfill -c(\mathbf{r})& =\int d{\mathbf{r}}^{\prime }\frac{\Gamma }{|\mathbf{r}-{\mathbf{r}}^{\prime }|}{g}_a({\mathbf{r}}^{\prime })=\Gamma {\tilde{v}}_L(r),\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill g_a(\mathbf{r}-{\mathbf{r}}^{\prime })& =\frac{e^{-|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^2/m}}{{(m\pi )}^{3/2}}=\frac{e^{-|{\mathbf{r}}_0-{\mathbf{r}}_0^{\prime }{|}^2/(m{a}^2)}}{{(m{a}^2\pi )}^{3/2}}(m={1.08}^{-2}),\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\tilde{v}}_L(r)& =\frac{\mathrm{erf}(1.08r)}{r}(r\equiv |\mathbf{r}|=|{\mathbf{r}}_0|/a),\hfill \end{array}$$

(28)

which gives

$$\begin{array}{cc}{\nabla }^2{\psi }_n(\mathbf{r},t)\hfill & \equiv -\int d{\mathbf{r}}^{\prime }{\nabla }^{2}c(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime },t)\hfill \\ \hfill & =-4\pi \Gamma \int d{\mathbf{r}}^{\prime }{g}_a(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime },t)\hfill \\ \hfill & =-4\pi \Gamma \tilde{n}(\mathbf{r},t),\hfill \end{array}$$

(29)

where $\tilde{n}(\mathbf{r},t)$ denotes a coarse-grained density that is smeared by the Gaussian distribution function $g_a(\mathbf{r}-{\mathbf{r}}^{\prime })$ over a range of the Wigner–Seitz radius a. The above form of the direct correlation function has been demonstrated to be available for the OCP in the strong coupling regime of $\Gamma \gg 1$. In Equation (26), the bare electrostatic potential ($\sim 1/r$) is modified using the Gaussian distribution function $g_a(\mathbf{r})$, and the second equation of Equation (26) introduces the function of ${\tilde{v}}_L(r)=\mathrm{erf}(1.08r)/r$ that represents the long-range part of the Coulomb interaction potential [31,32]. It is to be noted that the internal energy of the OCP, obtained using this direct correlation function, or Equation (26), exhibits an error of less than $0.8\%$ in the strong coupling regime [31].

Considering that the Fourier transform ${\tilde{v}}_L(k)$ of ${\tilde{v}}_L(r_0/a)$ is given by

$${\tilde{v}}_L(k)=\frac{4\pi }{k^2}{e}^{-k^2(m{a}^2)/4},$$

(30)

Equation (29) is rewritten in the original coordinate of ${\mathbf{r}}_0$ as

$$\begin{array}{cc}{\nabla }^2{\psi }_n({\mathbf{r}}_0,t)\hfill & =-k^2{\psi }_n(k,t)\hfill \\ \hfill & =-4\pi q^2{l}_B{e}^{-k^2(m{a}^2)/4}n(k,t),\hfill \end{array}$$

(31)

using the Fourier transforms of ${\psi }_n({\mathbf{r}}_0,t)$ and $n({\mathbf{r}}_0,t)$: ${\psi }_n(k,t)$ and $n(k,t)$. In the limit of $a\to 0$, Equation (29) is reduced to the conventional Poisson equation. It follows from Equations (29) and (31) that the Fourier transform $\tilde{n}(k,t)$ of the coarse-grained density $\tilde{n}({\mathbf{r}}_0,t)$ reads

$$\begin{array}{c}\hfill \tilde{n}(k,t)=e^{-k^2(m{a}^2)/4}n(k,t),\end{array}$$

(32)

implying the cut-off of $\tilde{n}(k,t)$ at a high wavenumber in correspondence with the coarse-graining of $n(\mathbf{r},t)$.

3.2. Implications of Longitudinal Contributions Given by Equation (33)

As described in Section 3.1, the gradient of ${\rho }_{\infty }(\mathbf{r})$ has only the z-component as found from Equation (25). The z-component given by ${\partial }_z{\rho }_{\infty }=-2\Gamma {\rho }_{\infty }$ yields the longitudinal contribution, $-\nabla ·{\mathbf{j}}_{\perp }(\mathbf{r},t)$, to the rhs of Equation (22):

$$\begin{array}{cc}-\nabla ·{\mathbf{j}}_{\perp }(\mathbf{r},t)\hfill & =-{\partial }_z{j}_z(\mathbf{r},t)\hfill \\ \hfill & =2\Gamma D{\partial }_{z}n(\mathbf{r},t)-D{\partial }_z{\rho }_{\infty }(z)\Gamma {\mathcal{E}}_z(\mathbf{r},t)\hfill \\ \hfill & =2\Gamma D\left\{{\partial }_{z}n(\mathbf{r},t)+{\rho }_{\infty }(z)\Gamma {\mathcal{E}}_z(\mathbf{r},t)\right\},\hfill \end{array}$$

(33)

where $\Gamma {\mathcal{E}}_z(\mathbf{r},t)$ denotes the z-component of fluctuating electric field $\mathbf{E}(\mathbf{r},t)$ defined by

$$\mathbf{E}(\mathbf{r},t)=-\nabla {\psi }_n(\mathbf{r},t)=\Gamma \left[\begin{array}{c}{\mathcal{E}}_x(\mathbf{r},t)\\ {\mathcal{E}}_y(\mathbf{r},t)\\ {\mathcal{E}}_z(\mathbf{r},t)\end{array}\right].$$

(34)

It is found from Equations (23), (33), and (34) that $-\nabla ·{\mathbf{j}}_{\perp }(\mathbf{r},t)$ disappears in the absence of the longitudinal variance in $n(\mathbf{r},t)$ (i.e., ${\partial }_{z}n(\mathbf{r},t)=0$) as well as ${\rho }_{\infty }(\mathbf{r})$, which is the reason why ${\mathbf{j}}_{\perp }(\mathbf{r},t)$ has been referred to as the longitudinal current. The two terms on the rhs of Equation (33) are expressed in the original coordinate as follows:

$$\begin{array}{cc}\hfill & 2\Gamma D{\partial }_{z}n(\mathbf{r},t)=\frac{a}{\lambda }\left(\frac{D_0}{a^2}\right){\partial }_{z}n(\mathbf{r},t)=\frac{D_0}{\lambda }{\partial }_{z_0}n({\mathbf{r}}_0,t),\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & 2\Gamma D{\rho }_{\infty }(z)\Gamma {\mathcal{E}}_z(\mathbf{r},t)=\frac{a}{\lambda }\left(\frac{D_0}{a^2}\right){\rho }_{\infty }(z)q^2{l}_B{\mathcal{E}}_z({\mathbf{r}}_0,t)=D_0{\rho }_{\infty }(z)\left\{\frac{q^2{l}_B{\mathcal{E}}_{z_0}({\mathbf{r}}_0,t)}{\lambda }\right\}.\hfill \end{array}$$

(36)

Equation (35) represents a vertical advection of fluctuating density field $n({\mathbf{r}}_0,t)$. Putting this advection term given by Equation (35) on the left hand side of Equation (22), the combination of Equations (33)–(36) transforms Equation (22) to

$$\begin{array}{c}\hfill {\partial }_{t}n({\mathbf{r}}_0,t)\underline{-\frac{D_0}{\lambda }{\partial }_{z_0}n({\mathbf{r}}_0,t)}=D_0{\nabla }^{2}n({\mathbf{r}}_0,t)-D_0{q}^2{l}_B{\rho }_{\infty }({\mathbf{r}}_0,t)\left\{4\pi \tilde{n}({\mathbf{r}}_0,t)\underline{-\frac{{\mathcal{E}}_{z_0}({\mathbf{r}}_0,t)}{\lambda }}\right\}+\zeta [{\rho }_{\infty },\overrightarrow{\eta }]\end{array}$$

(37)

in the original coordinate representation, where the underlined terms corresponding to the longitudinal contributions. The former contribution, the second term on the left hand side of Equation (37), suppresses density fluctuations, whereas the latter, the third term on the rhs of Equation (37), acts as a positive feedback to enhance counterion condensation in proximity to the charged plate. Figure 2 is a schematic of such opposite roles of longitudinal contributions from the above underlined terms.

On the one hand, the underlined term on the left hand side of Equation (37) represents the advective flow term. The advection velocity is given by $D_0/\lambda$, which increases as the Gouy–Chapman length $\lambda$, a characteristic length of the electric double layer, is shorter. The negative sign of this term indicates that the flow direction is always in the opposite direction to the z-axis. Figure 2 illustrates translation of whole fluctuating density field $n(\mathbf{r},t)$, like a Goldstone-mode, due to the advection flow. Figure 2 shows the case of ${\partial }_{z_0}n<0$ where the increase from ${\rho }_{\infty }(0)$ (i.e., $n({\mathbf{r}}_0,t)>0$ at $z_0=0$) is lowered when ${\partial }_{z_0}n<0$, and vice versa because of the fixed direction of the advection flow. In other words, density fluctuations are suppressed due to the former longitudinal contribution associated with advection flow.

On the other hand, the underlined contribution of the third term on the rhs of Equation (37) arises from the $z_0$-component $\Gamma {\mathcal{E}}_{z_0}$ of a fluctuating electric field $\mathbf{E}({\mathbf{r}}_0,t)$; however, this term reduces the third term on the rhs of Equation (37), or the smeared density $\tilde{n}({\mathbf{r}}_0,t)$ induced by $\mathbf{E}({\mathbf{r}}_0,t)$ itself. We should remember that, in the limit of $a\to 0$, the second term on the rhs of Equation (37) corresponds to the electrostatic term that is associated with the Poisson equation: we have $-D_0{q}^2{l}_B{\rho }_{\infty }({\mathbf{r}}_0,t)4\pi n({\mathbf{r}}_0,t)$ with $\tilde{n}({\mathbf{r}}_0,t)$ replaced by the original density $n({\mathbf{r}}_0,t)$, indicating the electrostatically restoring term to $n({\mathbf{r}}_0,t)\to 0$ as represented by the negative sign. Accordingly, the latter longitudinal contribution plays a role of positive feedback for counterion condensation, as opposed to the above electrostatic suppression. In actuality, the latter term increases the strength of longitudinal fluctuating field when ${\mathcal{E}}_{z_0}>0$, or ${\partial }_{z_0}n({\mathbf{r}}_0,t)<0$, which means that the counterions have become more condensed. Meanwhile, the negative sign of ${\mathcal{E}}_{z_0}<0$ yields the restoring contribution to $n({\mathbf{r}}_0,t)\to 0$ when accumulated counterions leave the charged plate: ${\partial }_{z_0}n({\mathbf{r}}_0,t)>0$. In both cases of ${\mathcal{E}}_{z_0}>0$ and ${\mathcal{E}}_{z_0}<0$, the latter longitudinal enhances counterion condensation, which is represented as electrical reverse-flow for ${\partial }_{z_0}n({\mathbf{r}}_0,t)<0$ in Figure 2.

4. Density-Density Correlations Due to Transverse Dynamics Along the Plate Surface

Supposing that ${\partial }_{z_0}n({\mathbf{r}}_0,t)=0$, we can focus on the transverse dynamics parallel to the charged plate at $z_0=0$, which we will investigate quantitatively. With the use of the Fourier transform $n(k,t)$ of $n({\mathbf{r}}_0,t)$, Equation (37) becomes

$$\begin{array}{cc}{\partial }_{t}n(k,t)\hfill & =-D_0\left\{k^2+4\pi q^2{l}_B{\rho }_{\infty }(0)e^{-k^2(m{a}^2)/4}\right\}n(k,t)+\zeta [{\rho }_{\infty },\overrightarrow{\eta }]\hfill \\ \hfill & =-D_0\Xi \mathcal{G}(k)n(k,t)+\zeta [{\rho }_{\infty },\overrightarrow{\eta }],\hfill \end{array}$$

(38)

$$\begin{array}{cc}\mathcal{G}(k)\hfill & =\frac{k^2}{\Xi }+4\pi \sigma e^{-k^2(m{a}^2)/4},\hfill \end{array}$$

(39)

given that ${\partial }_{z_0}n({\mathbf{r}}_0,t)=0$. In Equation (39), the conventional coupling constant $\Xi =q^2{l}_B/\lambda$ given by Equation (2) appears using ${\rho }_{\infty }(0)=\sigma /\lambda$. In both of the strong coupling regime ($\Xi \gg 1$) and the low wavenumber region ($ka\ll 1$), the propagator $\mathcal{G}(k)$ is approximated by

$$\begin{array}{c}\hfill \mathcal{G}(k)\approx \frac{1}{a^2}\left\{\frac{{(ka)}^2}{\Xi }+4q\right\}\approx \frac{4q}{a^2}=4\pi \sigma ,\end{array}$$

(40)

using the relation $\sigma =q/(\pi a^2)$.

Now we have the analytical solution to Equation (38) in the following form [25]:

$$\begin{array}{cc}\hfill n({\mathbf{r}}_0,t)& =e^{-t{D}_0\Xi \mathcal{G}({\mathbf{r}}_0)}n({\mathbf{r}}_0,0)+{\int }_0^{t}ds{e}^{-(t-s)D_0\Xi \mathcal{G}({\mathbf{r}}_0)}\zeta [{\rho }_{\infty }(z_0),\eta ({\mathbf{r}}_0,s)].\hfill \end{array}$$

(41)

Considering the real space representation that $D_0\Xi \mathcal{G}({\mathbf{r}}_0)\approx 4\pi D_0\Xi \sigma$ in the above approximation of Equation (40), the above exponential factors, $e^{-t{D}_0\Xi \mathcal{G}}$ and $e^{-(t-s)D_0\Xi \mathcal{G}}(s<t)$, are negligible due to $\Xi \gg 1$. Hence, Equation (41) is reduced to $n({\mathbf{r}}_0,t)=\zeta [{\rho }_{\infty }(z_0),\eta ({\mathbf{r}}_0,s)]$, thereby providing

$$\begin{array}{cc}〈n({\mathbf{r}}_0,t)n({\mathbf{r}}_0^{\prime },t^{\prime })〉\hfill & =〈\zeta [{\rho }_{\infty }(z_0),\overrightarrow{\eta }({\mathbf{r}}_0,t)]\zeta [{\rho }_{\infty }({\mathbf{r}}_0^{\prime }),\overrightarrow{\eta }({\mathbf{r}}_0^{\prime },t^{\prime })]〉\hfill \\ \hfill & =2D_0{\rho }_{\infty }(z_0)\delta ({\mathbf{r}}_0-{\mathbf{r}}_0^{\prime })\delta (t-t^{\prime }).\hfill \end{array}$$

(42)

It is found from Equation (42) that there is no correlation of transverse density fluctuations, which represents a coarse-grained frozen dynamics in the strong coupling regime of coarse-grained 2D OCP.

We can determine a crossover scale $l_c$. or associated crossover wavenumber $k_c=2\pi /l_c$ by comparing two terms on the rhs of Equation (39). In the above approximation, the first term on the rhs of Equation (39) has been neglected based on the condition that $\Xi \gg 1$ and $ka\ll 1$. While increasing the wavenumber and maintaining the strong coupling of $\Xi \gg 1$, we arrive at the crossover wavenumber $k_c$ that is defined by the following relation:

$$\begin{array}{c}\hfill \frac{k_c^2}{\Xi }=4\pi \sigma e^{-k_c^2(m{a}^2)/4},\end{array}$$

(43)

stating that the two terms on the rhs of Equation (39) are comparable to each other. We now introduce the main branch $W_0(x)$ of the Lambert W-function [34], so that Equation (43) is converted to

$$\begin{array}{cc}\nu \hfill & =W_0(\nu e^{\nu })=W_0(qm\Xi ),\hfill \\ \nu \hfill & \equiv k_c^2(m{a}^2)/4,\hfill \end{array}$$

(44)

based on another expression of Equation (43) as follows:

$$\begin{array}{c}\hfill \nu e^{\nu }=\pi \sigma \Xi (m{a}^2)=qm\Xi .\end{array}$$

(45)

The approximate form of $W_0(x)\approx lnx$ for $x\gg 1$ [34] applies to Equation (44) because of $qm\Xi \gg 1$. It follows that Equation (44) reads

$$\begin{array}{c}\hfill 2\pi \left(\frac{a}{l_c}\right)=k_{c}a\approx \frac{2}{m^{1/2}}\sqrt{ln(qm\Xi )}=2.16\sqrt{ln(qm\Xi )},\end{array}$$

(46)

which is our main result in this study. Below this scale specified by the crossover length $l_c$, we can observe the diffusive behavior of counterions, instead of frozen correlations represented by Equation (42). Taking $qm\Xi ={10}^4$ (or $\Xi \sim {10}^3$ for $q\sim 10$) as an example of the strong coupling regime, Equation (46) provides

$$\begin{array}{cc}{k}_{c}a\hfill & \approx 6.56,\hfill \\ \frac{l_c}{a}\hfill & \approx \frac{2\pi }{6.56}.\hfill \end{array}$$

(47)

The latter relation implies that the transverse dynamics of strongly-coupled counterions still retain diffusive behaviors within each Wigner–Seitz cell (i.e., $l_c\sim a$), which is physically plausible (see also Figure 3).

5. Summary and Conclusions

We have investigated stochastic density fluctuations $n=\rho -{\rho }_{\infty }$ around the ground state distribution (${\rho }_{\infty }\propto e^{-z_0/\lambda }$) of strongly coupled counterions near a single charged plate, focusing especially on the transverse dynamics parallel to the charged plate at $z_0=0$. The key to treating the stochastic dynamics is to use the DK equation of overdamped Brownian particles [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27] that is linearized by expanding the first derivative of a free energy functional of given density, ($\delta F\left[\rho \right]/\delta \rho$) around the ground state density ${\rho }_{\infty }$. As a result, we have obtained the linear DK Equation (37), which is applicable to the longitudinal and transverse dynamics of counterions in the strong coupling regime where the stationary density distribution has been investigated using Monte Carlo simulations [8,9,10,11].

The linear DK equation allows us to quantitatively investigate the dynamical crossover of transverse density fluctuations along the plate surface. Accordingly, we have found a crossover scale, given by Equations (46) and (47), above which the transverse density dynamics appear frozen, generating white noise that is uncorrelated with respect to time and space. Below the crossover scale, on the other hand, diffusive behavior of counterions can be observed along the plate surface, as illustrated in Figure 3. For instance, the crossover length $l_c$ is of the order of the Wigner–Seitz radius a when $\Xi \sim {10}^3$. Furthermore, the longitudinal dynamics vertical to the plate arises from the gradient of a fluctuating density field along the z-axis, producing additional contributions to the transverse dynamics, such as electrical reverse-flow as well as advective flow (see Figure 2). The electrical reverse-flow would be crucial in experimental situations where mobile ions, including not only counterions but also added salt, are affected considerably by the longitudinal dynamics. This remains to be addressed in a quantitative manner, by extending the present formulation to multi-component systems.