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

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.


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 counterions condensed on the plate. Our main result in this study is the quantitative evaluation of l c , yielding l c ∼ a for Ξ ∼ 10 3 . Section 5 contains a summary and conclusions.

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], where q is the valence of counterions, l B ≡ e 2 /4π 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 (πa 2 )σ = q using the surface number density σ 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 Γ increases, the OCP shows a transition from a weakly coupled gaseous regime (Γ 1) to a strongly coupled fluid regime (Γ 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 Ξ defined by [6][7][8][9][10][11] using the Gouy-Chapman length λ = 1/(2πql B σ), a characteristic length of the electric double layer. Inserting the relation σ = q/(πa 2 ) into the definition of λ, we have It follows from Equations (3) and (4) that When macroion-counterion attractions are weak, the structure of such ionic cloud, characterized by λ, will be dispersed instead of forming the 2D OCP. The dispersed electric double layer is thus represented by the weak coupling parameter of Ξ 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 Ξ increases while reducing λ, 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 λ in the strong coupling regime of Ξ = 2Γ 2 1 (i.e., a λ), as found from Equations (4) and (5). A field theoretical treatment provides counterion density distribution, ρ ∞ (z 0 ), in the ground state of the strong coupling limit (Ξ, Γ → ∞) [6,8]: where J(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(r 0 ) is expressed as J(r 0 ) = z 0 /λ with z 0 denoting the distance between the position 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 oneand 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 > λ; see [8] for a recent review.
2.2. Imposing a Given Density Distribution ρ on the Grand Potential Ω Letρ(r) be an instantaneous density of counterions located at r i (i = 1, · · · , N), where the counterion system is rescaled as r = (x, y, z) = (x 0 /a, y 0 /a, z 0 /a) = r 0 /a. Figure 1 shows a schematic of the rescaled system. The instantaneous density is expressed aŝ 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 which is a functional ofρ (see Appendix B for the details). Two schematics illustrating the side view of one charged plate system that consists of a positively charged plate carrying a surface charge density +σe and negatively charged counterions with valence q. Our scaling (z = z 0 /a) of an actual system depicted on the left side implies a coarse-grained system on the right hand side, where a denotes a mean separation between counterions, provided that all of the counterions are condensed on the oppositely charged plate uniformly: πa 2 σ = q. The Gouy-Chapman length λ ≡ 1/(2πql B σ) is also indicated. This paper adopts the coupling constant Γ, defined by Γ = q 2 l B /a that applies to the 2D one-component plasma (OCP), instead of the conventional one, Ξ = q 2 l B /λ, used for the counterion system.
We can impose a given density distribution ρ(r) on the counterion system via the following delta functional [14,[19][20][21]: where the potential field ψ(r) has been introduced in the Fourier transform of the delta functional.
Multiplying the configurational integral representation of the grand potential Ω[J] (see Appendix B for the definition) by the constraint in Equation (9), the formal expression of F[ρ] is obtained: where J(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) where the saddle-point potential field ψ * satisfies the following relation: The first Legendre transform of Ω[ψ * ] provides the Hohenberg-Kohn free energy A[ρ] defined by Equation (12), the central functional of the equilibrium density functional theory [28,29].

Stochastic Density Dynamics Obeying the Dean-Kawasaki Equation
Here we focus on the stochastic dynamics of a density field at time t, ρ(r, t), whose spatially varying distribution is the same as the coarse-grained variation ofρ(r). What matters in terms of the stochastic density dynamics is the free energy functional F[ρ] of a given density field ρ, rather than the grand potential Ω in equilibrium. For the density functional, we have provided the approximate form of F[ρ]. The driving force due to F[ρ] and the density-dependent multiplicative noise ζ[ρ, η] creates the stochastic dynamics that obeys the DK equation [13][14][15][16][17][18][19][20]: 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 with the vectorial white noise field η(r, t) that has the correlation η l (r, t)η m (r , t ) = δ lm δ(r − r )δ(t − t ). Equation (15) can read [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] In general, the stochastic Equation (14) includes not only the multiplicative noise term, but also the nonlinear term associated with F[ρ]. 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 ρ ∞ .

Linearizing the Stochastic Dean-Kawasaki Equation (14)
Expanding the first derivative of due to See Appendix D for the details. It follows from Equations (17) and (19) that the right-hand side (rhs) of Equation (14) reads when n/ρ ∞ 1. It is to be noted that the ideal gas distribution ρ ∞ (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(r, t) is governed by the strong Coulomb interactions as represented by the contribution, − dr c(r − r )n(r , t), on the rhs of Equation (20). In Equations (17) and (20), the direct correlation function c(r − r ) appears because F[ρ] is expressed using the Hohenberg-Kohn free energy functional A[ρ] as given by Equation (12) [28,29]. We also have due to ∂ t ρ ∞ = 0, on the left hand side of the DK Equation (14). Combining Equations (14), (20), and (21), we obtain the linear DK equation: where ψ n denotes a fluctuating Coulomb potential defined by and the longitudinal current j ⊥ (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 ∇ρ ∞ . Incidentally, in the rescaled system of ρ ∞ (r, t) = a 3 ρ ∞ (r 0 , t), Equation (6) is rewritten as in the rescaled system. In the next subsection, we will discuss the details of j ⊥ (r, t) associated with the presence of ∇ρ ∞ in the longitudinal direction. Equation (22) appears to be a simple extension of the Poisson-Nernst-Planck equation [30] when the longitudinal current j ⊥ (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]: which gives where n(r, t) denotes a coarse-grained density that is smeared by the Gaussian distribution function g a (r − r ) 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 Γ 1. In Equation (26), the bare electrostatic potential (∼ 1/r) is modified using the Gaussian distribution function g a (r), and the second equation of Equation (26) introduces the function of v L (r) = 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 v L (k) of v L (r 0 /a) is given by Equation (29) is rewritten in the original coordinate of r 0 as using the Fourier transforms of ψ n (r 0 , t) and n(r 0 , t): ψ n (k, t) and n(k, t). In the limit of a → 0, Equation (29) is reduced to the conventional Poisson equation. It follows from Equations (29) and (31) that the Fourier transform n(k, t) of the coarse-grained density n(r 0 , t) reads n(k, t) = e −k 2 (ma 2 )/4 n(k, t), implying the cut-off of n(k, t) at a high wavenumber in correspondence with the coarse-graining of n(r, t).

Implications of Longitudinal Contributions Given by Equation (33)
As described in Section 3.1, the gradient of ρ ∞ (r) has only the z-component as found from Equation (25). The z-component given by ∂ z ρ ∞ = −2Γρ ∞ yields the longitudinal contribution, −∇ · j ⊥ (r, t), to the rhs of Equation (22): where ΓE z (r, t) denotes the z-component of fluctuating electric field E(r, t) defined by It is found from Equations (23), (33), and (34) that −∇ · j ⊥ (r, t) disappears in the absence of the longitudinal variance in n(r, t) (i.e., ∂ z n(r, t) = 0) as well as ρ ∞ (r), which is the reason why j ⊥ (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: Equation (35) represents a vertical advection of fluctuating density field n(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 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 /λ, which increases as the Gouy-Chapman length λ, 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(r, t), like a Goldstone-mode, due to the advection flow. Figure 2 shows the case of ∂ z 0 n < 0 where the increase from ρ ∞ (0) (i.e., n(r 0 , t) > 0 at z 0 = 0) is lowered when ∂ 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 ΓE z 0 of a fluctuating electric field E(r 0 , t); however, this term reduces the third term on the rhs of Equation (37), or the smeared density n(r 0 , t) induced by E(r 0 , t) itself. We should remember that, in the limit of a → 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 ρ ∞ (r 0 , t)4πn(r 0 , t) with n(r 0 , t) replaced by the original density n(r 0 , t), indicating the electrostatically restoring term to n(r 0 , t) → 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 E z 0 > 0, or ∂ z 0 n(r 0 , t) < 0, which means that the counterions have become more condensed. Meanwhile, the negative sign of E z 0 < 0 yields the restoring contribution to n(r 0 , t) → 0 when accumulated counterions leave the charged plate: ∂ z 0 n(r 0 , t) > 0. In both cases of E z 0 > 0 and E z 0 < 0, the latter longitudinal enhances counterion condensation, which is represented as electrical reverse-flow for ∂ z 0 n(r 0 , t) < 0 in Figure 2. Here we consider the case that the fluctuating density n(r 0 , t) decreases with z 0 , and n(0.t) and n(λ, t) are abbreviations of n(r 0 , t)| z 0 =0 and n(r 0 , t)| z 0 =λ , respectively. The advective flow, or migration of density fluctuations as a whole, is always in the negative direction along the z 0 -axis. Meanwhile, the electrical reverse-flow is also in the negative direction, irrespective of the sign of ∂ z 0 n, or E z . Accordingly, the latter flow acts as a positive feedback of density fluctuations for ∂ z 0 n < 0.

Density-Density Correlations Due to Transverse Dynamics along the Plate Surface
Supposing that ∂ z 0 n(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(r 0 , t), Equation (37) becomes given that ∂ z 0 n(r 0 , t) = 0. In Equation (39), the conventional coupling constant Ξ = q 2 l B /λ given by Equation (2) appears using ρ ∞ (0) = σ/λ. In both of the strong coupling regime (Ξ 1) and the low wavenumber region (ka 1), the propagator G(k) is approximated by using the relation σ = q/(πa 2 ). Now we have the analytical solution to Equation (38) in the following form [25]: Considering the real space representation that D 0 ΞG(r 0 ) ≈ 4πD 0 Ξσ in the above approximation of Equation (40), the above exponential factors, e −tD 0 Ξ G and e −(t−s)D 0 Ξ G (s < t), are negligible due to Ξ 1. Hence, Equation (41) is reduced to n(r 0 , t) = ζ[ρ ∞ (z 0 ), η(r 0 , s)], thereby providing 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π/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 Ξ 1 and ka 1. While increasing the wavenumber and maintaining the strong coupling of Ξ 1, we arrive at the crossover wavenumber k c that is defined by the following relation: 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 based on another expression of Equation (43) as follows: The approximate form of W 0 (x) ≈ ln x for x 1 [34] applies to Equation (44) because of qmΞ 1. It follows that Equation (44) reads 2π a l c = k c a ≈ 2 m 1/2 ln(qmΞ) = 2.16 ln(qmΞ), 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Ξ = 10 4 (or Ξ ∼ 10 3 for q ∼ 10) as an example of the strong coupling regime, Equation (46) provides k c a ≈ 6.56, l c a ≈ 2π 6.56 .
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 ∼ a), which is physically plausible (see also Figure 3).

Summary and Conclusions
We have investigated stochastic density fluctuations n = ρ − ρ ∞ around the ground state distribution (ρ ∞ ∝ e −z 0 /λ ) 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, (δF[ρ]/δρ) around the ground state density ρ ∞ . 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 Ξ ∼ 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.

Conflicts of Interest:
The author declares no conflict of interest.

Appendix A. Electrostatic Interaction Energies: General Forms When Rescaled by the Wigner-Seitz Radius a
The counterion number N is imposed on the global electroneutrality condition, −qN + σΣ = 0 (Σ-the surface area of a macroion). The system has electrostatic interaction energy U = U cc + U cm + U mm , i.e., the sum of the counterion-counterion interaction energy U cc , the macroion-counterion interaction energy U cm , and the self-energy energy of macroion U mm due to electrostatic interactions between charged groups fixed on the macroion surface. In the thermal energy unit, we have where the interaction potential v(r 0 ) in the k B T unit is represented by an ion-ion separation vector, r 0 = (x 0 , y 0 , z 0 ), as v(r 0 ) ≡ 1/r 0 with r 0 = |r 0 | , and dS denotes that the integration of the macroion charge position R 0 is restricted to the surface of the macroion. Let us rescale the system as r = (x, y, z) = (x 0 /a, y 0 /a, z 0 /a) = r 0 /a, using the mean counterion-counterion separation a (or the Wigner-Seitz radius a) that is evaluated from the electrical neutrality condition, πa 2 σ = q, supposing that all counterions located on the macroion surface are uniformly distributed. Since the coupling constant Γ defined by Equation (1) gives ql B σ = q 2 l B /(πa 2 ) = Γ/(πa), Equation (A1) reads where we have used the relation dS = dS 0 /a 2 . Equation (A2) implies that the electrostatic interaction energies become extremely large in the strong coupling limit of Γ → ∞. The 2D OCP formed on the macroion surface has been often regarded as a ground state of strongly coupled counterions. For incorporating this analogy to the ground state of the 2D OCP into our formulation, we set the base energy U{R 1 , · · · , R N } given by U{R 1 , · · · , R N } = U cc {R 1 , · · · , R N } + U cm {R 1 , · · · , R N } + U mm . (A3) We write U mm , a constant interaction energy, as U mm = Nu m using the interaction energy u mm of one charged group located at a reference position of R ref , which is expressed as The base energy vanishes (U{R 1 , · · · , R N } ≈ 0) in the approximation that U cc {R 1 , · · · , R N } ≈ Nu mm U cm {R 1 , · · · , R N } ≈ −2Nu mm .
It is also convenient to introduce the reference energy U m associated with the charged groups on macroion surface: U{R 1 , · · · , R N } = U cc {R 1 , · · · , R N } + U m U m ≡ U cm {R 1 , · · · , R N } + U mm ≈ −Nu mm , where the reference energy U m can be regarded as a constant energy, as found from the above approximation in Equation (A5). In the Gaussian approximation, ∆F[ρ] corresponds to the logarithmic correction term [21], which we have neglected in Equation (11) as the first approximation of the strong coupling regime of Γ 1.