# Dynamic correlations in a polar fluid: confronting stochastic density functional theory to simulations

Sleeba Varghese, Pierre Illien, Benjamin Rotenberg

► To cite this version:

Sleeba Varghese, Pierre Illien, Benjamin Rotenberg. Dynamic correlations in a polar fluid: confronting stochastic density functional theory to simulations. The Journal of Chemical Physics, 2025, 163, pp.124107. 10.1063/5.0292306. hal-05272762

HAL Id: hal-05272762

https://hal.sorbonne-universite.fr/hal-05272762v1

Submitted on 22 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic correlations in a polar fluid: confronting stochastic density functional theory to simulations

Sleeba Varghese,¹ Pierre Illien,¹ and Benjamin Rotenberg¹, ², a)

¹) Sorbonne Université, CNRS, Laboratoire PHENIX (Physicochimie des Electrolytes et Nanosystèmes Interfaciaux), 4 place Jussieu, 75005 Paris, France

²) Réseau sur le Stockage Electrochimique de l'Energie (RS2E), FR CNRS 3459, 80039 Amiens Cedex, France

(Dated: 8 September 2025)

Understanding the dynamic behavior of polar fluids is essential for modeling complex systems such as electrolytes and biological media. In this work, we develop and apply a Stochastic Density Functional Theory (SDFT) framework to describe the polarization dynamics in the Stockmayer fluid, a prototypical model of dipolar liquids consisting of Lennard-Jones particles with embedded point dipoles. Starting from the overdamped Langevin dynamics of dipolar particles, we derive analytical expressions for the intermediate scattering functions and dynamic structure factors of the longitudinal and transverse components of the polarization field, within linearized SDFT. To assess the theory's validity, we compare its predictions with results from Brownian Dynamics simulations of the Stockmayer fluid. We find that SDFT captures the longitudinal polarization fluctuations accurately, while transverse fluctuations are underestimated due to the neglect of dipolar correlations. By incorporating the Kirkwood factor into a modified SDFT, we recover quantitative agreement for both components across a range of dipole strengths. This study highlights the utility of SDFT as a coarse-grained description of polar fluid dynamics and provides insights into the role of collective effects in polarization relaxation.

## I. INTRODUCTION

Understanding the microscopic behavior and organization of polar solvents, such as water, is a fundamental question in the field of chemical physics¹⁻³. This challenge has direct relevance to a broad spectrum of applications, ranging from the design of advanced batteries and fuel cells in electrochemistry⁴ to the behavior of charged soft matter systems such as colloids and biological fluids⁵⁻⁷. Accurate prediction of solvent structure and dynamics from first principles is especially difficult due to the complex interplay between thermal fluctuations, electrostatic and steric interactions, and the challenge of accounting for the orientation of the molecules.

From a theoretical point of view, two complementary approaches have emerged to address these questions. On the one hand, particle-based simulation methods, such as molecular dynamics, offer detailed insights into the time evolution of individual solvent molecules⁸⁻¹⁴. For instance, proposing accurate force fields to describe the structure and dynamics of water has been a central goal of computational physical chemistry during the past decades, resulting in a wealth of well-calibrated options¹⁵⁻¹⁹. On the other hand, analytical frameworks enable a more general, and often more computationally efficient, exploration of such systems. Their static structure can be accurately predicted e.g. using field-theoretical descriptions²⁰⁻²², or (molecular-) density functional theory²³⁻²⁶. Describing their dynamical response requires other tools, such as dynamical extensions of the density functional theory or mode-coupling-like theories.

For example, a theoretical model for the orientational dynamics of a polar solvent was proposed by Chandra

and Bagchi²⁷⁻²⁹. Building on the Smoluchowski-Vlasov equation³⁰, they developed a dynamical density functional theory (DDFT) framework for orientable particles, which under appropriate closure approximations, yields deterministic equations for the evolution of the polarization field. Subsequent studies refined this framework to account for the intermediate molecular interactions and investigated condensed matter phenomena such as polarization relaxation, solvation dynamics, and dielectric friction³¹⁻³⁵. A comprehensive review of classical dynamical density functional theory can be found in Ref. 36. While DDFT successfully establishes a link between microscopic and collective relaxation dynamics³⁷, its deterministic evolution equations allow it to only describe the average dynamical behavior of the system and provides no information on the fluctuations around this average. As a result, the DDFT framework becomes inadequate for probing the microscopic origins of fluctuation-induced phenomena in bulk or confined systems³⁸,³⁹.

A viable solution to the limitation of deterministic DDFT is to use more general theoretical models with intrinsic stochastic descriptions⁴⁰⁻⁴³. Within this category, stochastic density functional theory (SDFT) has become an important tool to investigate the dynamics of interacting entities, such as colloids or molecules, on the microscopic and nanoscopic scales⁴⁴. Building on the foundational works of Kawasaki⁴¹ and Dean⁴³, who introduced a formalism to coarse-grain the dynamics of interacting Langevin particles, SDFT provides a powerful theoretical framework to describe the time evolution of density fields under thermal noise. In physical chemistry, this framework was applied to study electrolyte solutions, and more specifically conductivity under various conditions³⁸,⁴⁵⁻⁵², spatial and temporal ionic correlations⁵³,⁵⁴, and viscosity in charged fluids⁵⁵⁻⁵⁷. This approach was later extended to include orientational degrees of freedom. In this context, the BBGKY hierarchy of equations and its different

<footnote>a) benjamin.rotenberg@sorbonne-universite.fr</footnote>
closures were discussed⁵⁸. Using DDFT-like closure schemes, static as well as dynamic correlations were computed from this approach⁵⁹⁻⁶¹, but were not compared to numerical simulations, therefore limiting the applicability of the approach for describing real systems.

In this context, we recently proposed a SDFT description of water⁶² as a solvent for ions. In this model, water molecules are represented as dipoles with translational and rotational motion that interact with each other through electrostatic interactions. Importantly, this model does not include any short-range repulsion between the dipoles. They are not meant to represent water molecules: instead, they represent 'polarized blobs', and the parameters of the model can be chosen in such a way that the target properties of water are reproduced. The main outcome of this work are the intermediate scattering functions associated with the longitudinal and transverse polarizations. Derived quantities, such as the power spectral density, are obtained straightforwardly. We showed that this simple representation of the solvent was valid in the limit of low frequencies and small wavevectors.

The goal of the present work is to investigate dynamical correlations in a standard model for dipolar liquids (Stockmayer fluid, consisting of particles interacting via point dipoles and Lennard-Jones potentials), and to assess how well these correlations are captured with SDFT in the linearization approximation and neglecting short-range nonelectrostatic interactions, for which analytical results were obtained in Ref. 62.

To this end, we first derive analytical expressions for the time-dependent polarization fluctuations from SDFT, building upon our previous framework for polar solvents. We then perform Brownian Dynamics (BD) simulations of the Stockmayer fluid to serve as a microscopic reference. By systematically comparing SDFT predictions with simulation data across various dipole moments and diffusion coefficients, we identify the strengths and limitations of the theory. We find that while SDFT reproduces the longitudinal polarization dynamics quantitatively, it underestimates transverse fluctuations due to the neglect of dipole-dipole correlations. By incorporating the Kirkwood g-factor into a modified SDFT formulation, we achieve much improved agreement, particularly for the transverse component. These results establish SDFT as a viable coarse-grained model for capturing the essential features of polarization dynamics in dipolar fluids, and provide a framework for extending this approach to more complex polar media.

## II. POLARIZATION FLUCTUATIONS IN POLAR FLUIDS

### A. Microscopic description

The system of interest is made of N solvent molecules in a volume $V$, represented as dipoles of dipolar moment $p$. Their number density is indicated by $C_s = N/V$, and their positions and unit orientations at time $t$ are denoted by $\mathbf{r}_1(t), \dots, \mathbf{r}_N(t)$ and $\hat{\mathbf{u}}_1(t), \dots, \hat{\mathbf{u}}_N(t)$ respectively. The microscopic polarization

density is defined as:

$$P(\mathbf{r},t) = \sum_{i=1}^{N} p \hat{u}_i(t) \delta(\mathbf{r} - \mathbf{r}_i(t)). \quad (1)$$

In Fourier space, the polarization vector field can be decomposed into a longitudinal and transverse parts (with respect to the wavevector $\mathbf{q}$), with Cartesian components

$$\tilde{P}_a(\mathbf{q},t) = \underbrace{\frac{q_a q_b}{q^2} \tilde{P}_b(\mathbf{q},t)}_{\equiv \tilde{P}_{L,a}(\mathbf{q},t)} + \underbrace{\left(\delta_{ab} - \frac{q_a q_b}{q^2}\right) \tilde{P}_b(\mathbf{q},t)}_{\equiv \tilde{P}_{T,a}(\mathbf{q},t)}, \quad (2)$$

where an implicit sum over index $b$ is assumed. Note that, throughout the paper, we will use the following convention for Fourier transformation: $\tilde{f}(\mathbf{q}) = \int d\mathbf{r} e^{-i\mathbf{q} \cdot \mathbf{r}} f(\mathbf{r})$ and $f(\mathbf{r}) = \int \frac{dq}{(2\pi)^3} e^{i\mathbf{q} \cdot \mathbf{r}} \tilde{f}(\mathbf{q})$. Explicitly, the longitudinal component $\tilde{P}_L$ is the projection of $\tilde{P}(\mathbf{q},t)$ in the direction of the wavevector $\mathbf{q}$, and the transverse component is the remainder $\tilde{P}_T = \tilde{P} - \tilde{P}_L$.

The spatio-temporal fluctuations of these components of the polarization are described by the corresponding intermediate scattering functions (ISF),

$$F_{L,T}(\mathbf{q},t) = \frac{1}{N} \langle \tilde{P}_{L,T}(\mathbf{q},t) \cdot \tilde{P}_{L,T}(-\mathbf{q},0) \rangle, \quad (3)$$

where $\langle \dots \rangle$ denotes an ensemble average and we have used the fact that for a bulk fluid the ISF only depends on the magnitude $q$ but not on the orientation. The initial value of the ISF is simply the static polarization structure factor

$$S_{L,T}(\mathbf{q}) = \frac{1}{N} \langle |\tilde{P}_{L,T}(\mathbf{q})|^2 \rangle. \quad (4)$$

The static fluctuations of the total dipole $\mathbf{M} = \tilde{P}(\mathbf{q} = 0) = \sum_i p \hat{u}_i$ are often characterized by the so-called Kirkwood g-factor:

$$g_K = \frac{\langle \mathbf{M}^2 \rangle}{N p^2}, \quad (5)$$

which is a measure of the local orientational correlations between neighboring dipoles⁶³⁻⁶⁵. For uncorrelated dipoles, $\langle \mathbf{M}^2 \rangle = Np^2$ and the Kirkwood factor is equal to 1.

The dynamics is often also considered in the frequency domain, introducing the Laplace transforms with (with Laplace variable $s = -i\omega$), $\tilde{\mathcal{F}}(\omega) = \int_0^\infty f(t) e^{+i\omega t} dt$, or the power spectral densities (PSD) also known as dynamic structure factors:

$$S_{L,T}(\mathbf{q}, \omega) = \int_{-\infty}^{\infty} F_{L,T}(\mathbf{q}, t) e^{+i\omega t} dt = 2\Re[\tilde{F}_{L,T}(\mathbf{q}, \omega)] \quad (6)$$

where the last equality follows from the fact $F_{L,T}(\mathbf{q},t)$ are even functions of time.

### B. Properties of interest

The ISF and PSD, which describe the equilibrium fluctuations of the polarization, are directly related to important
quantities describing the linear response of the polar fluid to an external field $E_{\text{ext}}(q, \omega)$. We refer the reader to Ref. 38 for a more comprehensive discussion of electrical fluctuations in electrolytes and their link between a variety of observable properties, but for the present work we consider the permittivity tensor, which relates the polarization to the local electric field $P(q, \omega) = \epsilon_0 [\epsilon(q, \omega) - \mathbf{I}] \cdot E(q, \omega)$, with $\epsilon_0$ the vacuum permittivity. The link between its longitudinal and transverse components and the corresponding susceptibilities depend on the boundary conditions (in particular in simulations using periodic boundary conditions). For $q \neq 0$, one has$^{66-68}$

$$ \chi_L(q, \omega) = 1 - \frac{1}{\epsilon_L(q, \omega)}, \quad (7) $$

$$ \chi_T(q, \omega) = \epsilon_T(q, \omega) - 1. \quad (8) $$

For non-polarizable systems the longitudinal susceptibility is related to the ISF Eq. 3 and its Laplace transform as$^{3,65,69,70}$

$$ \chi_{L,T}(q, \omega) = \frac{\beta N}{V \epsilon_0 \nu_{L,T}} \left[ S_{L,T}(q) + i \omega \tilde{F}_{L,T}(q, \omega) \right], \quad (9) $$

where $\beta = 1/k_B T$ ($k_B$ is the Boltzmann constant and $T$ is the temperature), and where $\nu_L = 1$ and $\nu_T = 2$. The above results establish a direct link between the microscopic structure and dynamics of the system with key properties of polar fluids. Decades of theoretical approaches with various levels of refinement and computer simulations have already provided many important results.

The static permittivity $\epsilon_r \equiv \epsilon(q=0, \omega=0)$ of the system described in Section II A, which is a material property, is related to the static fluctuations of total dipole (see e.g. Ref. 65). This link is rather subtle, because the polarization fluctuations also depend on the shape of the considered system and the dielectric properties of the medium in which it is “embedded”. This is particularly important to compute the permittivity in simulations under periodic boundary conditions. Specifically, for a spherical sample embedded in a medium with permittivity $\epsilon'$ (see e.g. Refs. 71,72), one has

$$ \frac{(\epsilon_r - 1)(2\epsilon' + 1)}{2\epsilon' + \epsilon_r} = 3yg_K, \quad (10) $$

where $g_K$ defined in Eq. (5) and $y$ is a dimensionless parameter often found in the literature on polar systems

$$ y \equiv \frac{C_s p^2}{9 \epsilon_0 k_B T} = \frac{C_s \alpha_{\text{or}}}{3}. \quad (11) $$

In the second expression, we have introduced the orientational polarizability $\alpha_{\text{or}} = p^2/3\epsilon_0 k_B T$ characterizing the linear response of a single dipole at temperature $T$ to an external field. This quantity is a characteristic volume and its product with the number density $C_s$ in the $y$ parameter controls the importance of dipolar couplings in the system, just as the product $C_s \Lambda^3$, with $\Lambda$ the de Broglie thermal wavelength, controls that of quantum effects, or the product $C_s \sigma^3$, with $\sigma$ a particle diameter, controls that of excluded volume.

When the permittivity of the embedding medium is the same permittivity as the system of interest ($\epsilon' = \epsilon_r$), Eq. (10)

results in the Kirkwood formula: $(\epsilon_r - 1)(2\epsilon_r + 1)/\epsilon_r = 9yg_K$, whereas for $\epsilon' = 1$ one recovers the Clausius-Mosotti equation $(\epsilon_r - 1)/(\epsilon_r + 2) = yg_K$, and in the limit $\epsilon' \to \infty$ (so-called “tin-foil” boundary conditions in Ewald summation used in most simulations), Eq. (10) reduces to

$$ \epsilon_r - 1 = 3yg_K. \quad (12) $$

Another important physical insight is the anisotropy of the dipole fluctuations. For example, the static longitudinal fluctuations of the dipole are smaller than the transverse ones$^{65}$, by a factor $\epsilon_r$ in the $q \to 0$ limit:

$$ \lim_{q \to 0} \frac{S_L(q)}{S_T(q)} = \frac{1}{\epsilon_r}. \quad (13) $$

The dynamics of polarization fluctuations can be characterized by effective relaxation times defined as the integral of the normalized correlation functions,

$$ \tau_{L,T}(q) = \int_0^\infty \frac{F_{L,T}(q,t)}{S_{L,T}(q)} dt = \frac{\tilde{F}_{L,T}(q,0)}{S_{L,T}(q)}. \quad (14) $$

The longitudinal polarization fluctuations decay faster than the transverse ones, with a ratio of characteristic times$^{37,65,73}$

$$ \lim_{q \to 0} \frac{\tau_L(q)}{\tau_T(q)} = \frac{1}{\epsilon_r}. \quad (15) $$

From a microscopic perspective, this can be understood in the limit of a large permittivity: whereas parallel components of the dipoles (i.e. their contribution to the longitudinal quantities) tend to compensate each other and therefore decrease the time it takes to modify $S_L(q=0)$ significantly, their perpendicular components tend to reinforce each other, and therefore make their relaxation slower$^{65,73}$.

Experimentally, one can only measure the $q \to 0$ limits of the longitudinal and transverse components of the permittivity tensor, which coincide. It is customary to analyze the permittivity using the Debye model (or a sum of Debye modes):

$$ \begin{aligned} \epsilon(\omega) &= \lim_{q \to 0} \epsilon_L(q, \omega) = \lim_{q \to 0} \epsilon_T^{\text{SDFT}}(q, \omega) \\ &= 1 + \frac{\epsilon_r - 1}{1 - i\omega\tau_D}, \end{aligned} \quad (16) $$

where the Debye time $\tau_D$ characterizes the dynamics of collective polarization fluctuations. This time can also be related to the individual dipole dynamics (in the presence of the other dipoles), characterized by:

$$ \tau_p = \int_0^\infty \langle \hat{u}_i(t) \cdot \hat{u}_i(0) \rangle dt, \quad (17) $$

where the unit vector describing the orientation of dipole $i$ follows Eq. (19). Note that this is not the definition in Ref. 65, even though it must coincide for exponential relaxation.

### III. STOCHASTIC DENSITY FUNCTIONAL THEORY FOR A POLAR FLUID

In this section, we introduce the SDFT description for the dynamics in polar fluids, corresponding to the ion-free case of
the more general results presented in Ref. 62 for electrolytes. We then summarize the main results for the polarization fluctuations in the absence of ions, which had not been reported in our previous work. Finally, we discuss the relevance of this level of theory to describe polar fluids, in order to clarify the points to be addressed by comparison with microscopic simulations of such systems.

### A. Governing equations

In the overdamped limit, the evolution equations of the positions and orientations read:

$$ \frac{d\mathbf{r}_i}{dt} = -\mu_s(p\hat{\mathbf{u}}_i \cdot \nabla)\nabla\varphi(\mathbf{r}_i) - \mu_s \sum_{j \neq i} \nabla V(\hat{\mathbf{u}}_i, \hat{\mathbf{u}}_j, \mathbf{r}_i - \mathbf{r}_j) + \sqrt{2D_s}\xi_i^t(t) \quad (18) $$

$$ \frac{d\hat{\mathbf{u}}_i}{dt} = \left\{ -[\mu_s^r p \hat{\mathbf{u}}_i \times \nabla \varphi(\mathbf{r}_i)] + \sqrt{2D_s} \xi_i^r(t) \right\} \times \hat{\mathbf{u}}_i \quad (19) $$

where $D_s$ is the translational diffusion coefficient and $D_s^r$ is the rotational diffusion coefficient of the polar molecules, with their associated mobilities defined as $\mu_s = \beta D_s$ and $\mu_s^r = \beta D_s^r$. The translational and orientational noises $\xi_i^t(t)$ and $\xi_i^r(t)$ are uncorrelated Gaussian white noises of zero average and unit variance, i.e. $\langle \xi_{i,n}^\alpha(t) \xi_{j,m}^\gamma(t') \rangle = \delta_{\alpha\gamma}\delta_{mn}\delta_{ij}\delta(t-t')$, where n or m are components of the vectors and Greek letters refer to particle labels. $V(r)$ is the interparticle potential, that may account for short-range steric repulsion – this term is not included in our SDFT analysis, but will be considered in the numerical simulations. $\varphi(r)$ is the electrostatic potential at the position r respectively, and is obtained by solving Poisson’s equation,

$$ -\nabla^2 \varphi(\mathbf{r}) = \frac{\rho_s(\mathbf{r})}{\varepsilon_0}, \quad (20) $$

where $\rho_s$ is the charge density associated with the polar molecules.

### B. Polarization fluctuations

The microscopic polarization density is defined as:

$$ \mathbf{P}(\mathbf{r}, t) = \sum_{i=1}^{N} p \hat{\mathbf{u}}_i(t) \delta(\mathbf{r} - \mathbf{r}_i(t)). \quad (21) $$

The charge density will be considered within the dipolar approximation, so that $\rho_s \simeq -\nabla \cdot \mathbf{P}$. Incorporating the orientational degrees of freedom to the standard SDFT approach introduced by Dean$^{43}$, and after linearizing around a constant, isotropic and homogeneous state, we obtain the general equation for the evolution of the polarization field as$^{62}$

$$ \partial_t \mathbf{P}(\mathbf{r}, t) = D_s \nabla^2 \mathbf{P} - 2D_s^r \mathbf{P} + \frac{1}{3} p^2 C_s \nabla (\mu_s \nabla^2 \varphi - 2\mu_s^r \varphi) + \Xi(\mathbf{r}, t), \quad (22) $$

where $\Xi$ is a Gaussian white noise of zero average and of variance:

$$ \begin{aligned} \langle \Xi_k(\mathbf{r}, t) \Xi_l(\mathbf{r}', t') \rangle &= \frac{2p^2 C_s}{3} \delta_{kl} \delta(t-t') (-D_s \nabla^2 + 2D_s^r) \delta(\mathbf{r}-\mathbf{r}') . \end{aligned} \quad (23) $$

Starting from Eq. (22), we obtain the evolution equations for the longitudinal and transverse components of the polarization field in Fourier space as

$$ \partial_t \tilde{\mathbf{P}}_L(\mathbf{q}, t) = -\frac{1}{\tau_L(q)} \tilde{\mathbf{P}}_L + \Xi_L(\mathbf{q}, t), \quad (24) $$

$$ \partial_t \tilde{\mathbf{P}}_T(\mathbf{q}, t) = -\frac{1}{\tau_T(q)} \tilde{\mathbf{P}}_T + \Xi_T(\mathbf{q}, t), \quad (25) $$

where the longitudinal and transverse relaxation rates are

$$ \frac{1}{\tau_L(q)} = \frac{1 + 3y}{\tau_s^r} (1 + q^2 a^2), \quad (26) $$

$$ \frac{1}{\tau_T(q)} = \frac{1}{\tau_s^r} (1 + q^2 a^2). \quad (27) $$

The dimensionless parameter $y$ is defined in Eq. (11) and we have introduced the rotational time $\tau_s^r = 1/2D_s^r$ for an isolated dipole obeying Eq. (19), and a length scale emerging from the ratio of the translational and diffusion coefficients, $a \equiv \sqrt{D_s/2D_s^r}$. The latter defines a crossover between two regimes: for $qa \ll 1$, one can consider dipole rotation as sufficiently fast to be averaged, whereas for $qa \gg 1$ rotation plays a role in the spatial fluctuations of the polarization. Finally, the noises $\Xi_L(q, t)$ and $\Xi_T(q, t)$ are uncorrelated and have respective variances:

$$ \langle \Xi_{L,i}(\mathbf{q}, t) \Xi_{L,j}(\mathbf{q}', t') \rangle = 2 \frac{q_i q_j}{q^2} D_s k_B T \kappa_s(q)^2 \varepsilon_0 (2\pi)^3 \delta(\mathbf{q} + \mathbf{q}') \delta(t - t') \quad (28) $$

and

$$ \langle \Xi_{T,i}(\mathbf{q}, t) \Xi_{T,j}(\mathbf{q}', t') \rangle = 2 \left( \delta_{ij} - \frac{q_i q_j}{q^2} \right) D_s k_B T \kappa_s(q)^2 \varepsilon_0 (2\pi)^3 \delta(\mathbf{q} + \mathbf{q}') \delta(t - t'). \quad (29) $$

where we have introduced the reciprocal square of a wavenumber-dependent screening length associated with the polarization charge:

$$ \kappa_s(q)^2 = \frac{p^2 C_s}{3\varepsilon_0 k_B T a^2} (1 + q^2 a^2) = \frac{3y}{a^2} (1 + q^2 a^2). \quad (30) $$

From Eqs. (24)-(25), one gets the expressions for correlation functions for the polarization densities in the infinite space limit as

$$ \langle \tilde{\mathbf{P}}_L(\mathbf{q}, t) \cdot \tilde{\mathbf{P}}_L(\mathbf{q}', t') \rangle = \frac{(2\pi)^3 \delta(\mathbf{q} + \mathbf{q}') p^2 C_s}{3(1 + 3y)} e^{-|t-t'|/\tau_L(q)}, \quad (31) $$

$$ \langle \tilde{\mathbf{P}}_T(\mathbf{q}, t) \cdot \tilde{\mathbf{P}}_T(\mathbf{q}', t') \rangle = 2 \frac{(2\pi)^3 \delta(\mathbf{q} + \mathbf{q}') p^2 C_s}{3} e^{-|t-t'|/\tau_T(q)}. \quad (32) $$
For a system with periodic boundary conditions, the corresponding expressions are obtained by changing

$$ (2\pi)^3 \delta(q+q') \rightarrow \gamma \delta_{q,-q'} \quad (33) $$

and $C_s = N/\gamma$, resulting in intermediate scattering functions (see Eq. (3)) decaying as single exponential:

$$ F_L(q,t) = \frac{p^2}{3} \frac{1}{1+3y} e^{-t/\tau_L(q)}, \quad (34) $$

$$ F_T(q,t) = \frac{2p^2}{3} e^{-t/\tau_T(q)}. \quad (35) $$

Note that the characteristic times $\tau_{L,T}(q)$ introduced in Eqs. (26) and (27) are thus consistent with the definition of the effective relaxation times Eq. (14). The dynamic structure factors (see Eq. (6)) are obtained from Eqs. (34) and (35), resulting in:

$$ S_L(q, \omega) = \frac{p^2}{3} \frac{1}{1+3y} \frac{2\tau_L(q)}{1+\omega^2\tau_L(q)^2}, \quad (36) $$

$$ S_T(q, \omega) = \frac{2p^2}{3} \frac{2\tau_T(q)}{1+\omega^2\tau_T(q)^2}. \quad (37) $$

Finally, we note that the solution of the Poisson equation (20) in reciprocal space using Fourier transforms corresponds to Fourier series for a periodic system in the limit of an infinite period. Since we did not consider an additional “surface” term (see e.g. Refs. 71,72), the above results correspond to the tin-foil boundary conditions mentioned in Section II B. In particular, the corresponding permittivity is given by Eq. (12).

### C. Relevance of SDFT for polar fluids

As mentioned in Sections III A and III B, the above analytical results for the polarization fluctuations were obtained after several approximations. In particular, non-electrostatic interactions are neglected and the dipoles are considered to interact only via the mean-field electrostatic potential induced by the other dipoles. In addition, the coupling between dipoles is considered to be sufficiently small that one can linearize the equations around a constant, isotropic, and homogeneous state. Here, we discuss the relevance of these results for the description of polar fluids, before turning to an explicit comparison with particle-based simulations in Section IV.

We first observe that the SDFT predictions for the ISF (Eqs. (34)-(35)) and PSD (Eqs. (36)-(37)) are consistent with the general results for the $q \to 0$ limits of the ratios between longitudinal and transverse static structure factors (Eq. (13)) and relaxation times (Eq. (15)), provided that the permittivity of the fluid is equal to

$$ \varepsilon_r^{\text{SDFT}} \equiv 1 + 3y. \quad (38) $$

Interestingly, Eqs. (34)-(35) and (36)-(37) show that in the SDFT case, the relations (13) and (15) also hold for all $q$. This is likely due to the fact that the approximations leading to

these results neglect short range correlations and correspond to the hydrodynamic regime $q \to 0, \omega \to 0$. In turn, Eq. (12) then leads to

$$ g_K^{\text{SDFT}} = 1. \quad (39) $$

As noted below its definition, Eq. (5), the Kirkwood factor is a measure of the local orientational correlations. The absence of short-range correlations in SDFT is consistent with the mean-field treatment of interactions.

The SDFT predictions Eqs. (34)-(37) depend on only 4 independent physical parameters: the dipole moment $p$; the dimensionless parameter $y$, or equivalently the number density $C_s$ for a given dipole moment and temperature, characterizing the importance of dipolar couplings; the rotational diffusion coefficient $D_s^r$ corresponding to the rotational time for an isolated dipole; the characteristic length $a$, or equivalently for a given $D_s^r$, the translational diffusion coefficient $D_s$ (see below Eqs. (26)-(27)).

Introducing Eqs. (34)-(35) in Eq. (9) provides the longitudinal and transverse susceptibilities. Eqs. (7) and (8) then yield the corresponding components of the wavelength- and frequency-dependent permittivity tensor. The SDFT prediction for the frequency-dependent permittivity is:

$$ \begin{aligned} \varepsilon^{\text{SDFT}}(\omega) &= \lim_{q \to 0} \varepsilon_L^{\text{SDFT}}(q, \omega) = \lim_{q \to 0} \varepsilon_T^{\text{SDFT}}(q, \omega) \\ &= 1 + \frac{3y}{1 - i\omega\tau_D}, \end{aligned} \quad (40) $$

i.e. of the Debye form, with a Debye relaxation time

$$ \tau_D = \tau_s^r = \frac{1}{2D_s^r}. \quad (41) $$

This result is consistent with previous results within the Mean Spherical Approximation³⁷, in which, in contrast with the present approach, short-range interactions between solvent particles are accounted for by setting the pair distribution function to zero below the radius of the particles. This leads to $\tau_D = \tau_T(q \to 0) = \varepsilon_r \tau_L(q \to 0)$, since in the present case $\tau_T(q \to 0) = \varepsilon_r \tau_L(q \to 0) = \tau_s^r$, i.e. the rotational time for an isolated dipole. However, the fact that it corresponds to the individual relaxation time for rotation is a consequence of the approximations made to model the polarization fluctuations. Indeed, the dynamics of collective fluctuations are in general not identical to that of individual dipoles, and a dynamic Kirkwood factor $\dot{g}$ can be introduced to link the Debye and individual relaxation times⁶⁵. We still need to check that SDFT predicts that $\tau_p^{\text{SDFT}}$, defined by Eq. (17), is equal $\tau_s^r$.

In Ref. 62, we used the experimental permittivity spectrum of water, which is well described in the low-frequency domain by a single Debye mode with $\tau_D \approx 10$ ps and a relative permittivity $\varepsilon_r = 78.5$, together with Eqs. (40) and (41), to obtain values of $D_s^r = 0.05$ ps⁻¹ and y. Combining $D_s^r$ with the experimental value of the translational (self-)diffusion of water at 25°C ($D_s = 2.310^{-9}$ m²s⁻¹), we estimated the lengthscale $a \approx 2.14$ Å. Finally, assuming a number density equal to that of water at the same temperature, $C_s = 55$ M (0.033 Å⁻³), then enforces the dipole $p = 1.4$ D consistent with the value
of y determined from the experimental permittivity. While this procedure is consistent with experimental data and results in parameters compatible with physical intuition, a better strategy to assess the range of validity and determine the effective parameters of the coarse-grained theory (SDFT) is to build the latter from a microscopic basis. Therefore, in the following we introduce such a microscopic reference obtained by simulations of a more realistic model of a dipolar fluid.

## IV. BROWNIAN DYNAMICS SIMULATIONS

### A. Stockmayer fluid

In order to obtain 'exact' results for the polarization fluctuations of a polar fluid with microscopic dynamics described by the overdamped Langevin Eqs. (18)-(19) without resorting to the approximations used in SDFT, it is possible to integrate them numerically in Brownian Dynamics simulations. A standard model for polar fluids is the Stockmayer fluid, consisting of point dipoles $p\hat{u}_i$, where $\hat{u}_i$ is the unitary orientation vector of molecule $i$ and $p$ the dipole moment magnitude of the particle, also interacting via Lennard-Jones (LJ) potentials to account for short-range range repulsion as well as dispersion between particles. The pairwise Stockmayer potential between particles $i$ and $j$ separated by a distance $r$ is thus

$$V(\hat{u}_i, \hat{u}_j, r) = U_{\text{LJ}}(r) + U_{\text{dipole}}(\hat{u}_i, \hat{u}_j, r) \quad (42)$$

where

$$U_{\text{LJ}}(r) = \begin{cases} 4\epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right] & \text{if } r \le r_c, \\ 0 & \text{if } r > r_c, \end{cases} \quad (43)$$

with $\sigma$ and $\epsilon$ are the LJ distance and energy, respectively, and $r_c$ is a cut-off radius, and

$$U_{\text{dipole}}(\hat{u}_i, \hat{u}_j, r) = \frac{p^2}{4\pi\epsilon_0} \left[ \frac{1}{r^3}(\hat{u}_i \cdot \hat{u}_j) - \frac{3}{r^5}(\hat{u}_i \cdot r)(\hat{u}_j \cdot r) \right] \quad (44)$$

In MD simulations of the Stockmayer fluid$^{8–14}$, short-range repulsion between particles results in momentum transfer by collision and momentum conservation results in hydrodynamic effects such as long-time tails in the velocity autocorrelation functions$^{74}$. In BD simulations, momentum is not conserved because of the energy exchange with the thermal bath, so that such effects may not be captured. The overdamped dynamics of Eqs. (18) and (19) do not describe short-time features arising from the inertia of particles, but rather introduce the effect of such high-frequency collective modes on the dynamics at intermediate and long times via the friction and stochastic forces acting on each particle.

From the trajectory of Stockmayer particles, we measure the longitudinal and transverse components of the polarization

$$\tilde{P}(q,t) = \sum_{i=1}^{N} p \hat{u}_i(t) e^{-iq \cdot r_i(t)} \quad (45)$$

<table><caption>TABLE I. Parameters common to all simulations of the Stockmayer fluid. See text for the values of the dipole, rotational and translational diffusion coefficients.</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Value</th></tr></thead><tbody><tr><td>σ</td><td>LJ diameter</td><td>3.024 Å</td></tr><tr><td>ε</td><td>LJ energy</td><td>0.4412 kcal mol<sup>-1</sup></td></tr><tr><td>T</td><td>temperature</td><td>298 K</td></tr><tr><td>C<sub>s</sub></td><td>number density</td><td>0.033 Å<sup>-3</sup></td></tr></tbody></table>

for selected wavevectors $q$ compatible with the periodic boundary conditions, from which we compute the corresponding ISF Eqs. (34)-(35) and PSD Eqs. (36)-(37). In order to assess the effects of these physical parameters, we consider several systems with different dipole moments $p$ and rotational/translational diffusion coefficients $D_s^r$ and $D_s$, keeping the number density $C_s$ and temperature $T$ constant.

### B. Simulation Details

Table I summarizes the simulation parameters common to all considered systems. In all case, we consider a fixed temperature $T = 298$ K and number density $C_s = 0.033 \text{ Å}^{-3}$ (experimental value for water at this temperature and a pressure of 1 atm), which corresponds to a reduced density $C_s^* = C_s \sigma^3 \approx 0.9$ for the present choice of LJ diameter. The LJ interaction parameters for the Stockmayer fluid were taken from Ramirez et al.$^{75}$. As a reference system, we consider particles with a dipole moment $p = 1.85$ D (or in reduced units $p^* = p\sqrt{\beta/\epsilon_0\sigma^3} = 6.15$), corresponding to Ref. 26 (comparing MD simulations of the Stockmayer fluid with the same LJ parameters to molecular Density Functional Theory as a simple model for water), translational diffusion coefficient $D_s = 2.310^{-9} \text{ m}^2\text{s}^{-1}$ (experimental value for water) and rotational diffusion coefficient $D_s^r = 0.05 \text{ ps}^{-1}$ (deduced from the experimental Debye relaxation time of water in Ref. 62). This reference system corresponds to a dipolar coupling parameter (see Eq. 11) $y = 3.84$. The corresponding permittivity $\epsilon_r \approx 140$ (see Ref. 26) is much larger than that of water, reflecting the importance of dipolar couplings in this system. In order to probe a wide range of regimes, we also considered systems with dipole moments in the range $p \in [0.148, 1.85]$ D, corresponding to $y \in [0.025, 3.84]$.

Short-range LJ and dipole-dipole interactions are computed using Eq. (43) and (44), respectively, with a real-space cutoff $r_c = 7.56 \text{ Å} \approx 2.5\sigma$. Long-range dipolar interactions are computed using the P3M algorithm$^{76}$ with a relative root mean square error in the per-atom force calculations below $10^{-4}$. Brownian dynamics (BD) simulations were performed using the large atomic/molecular massively parallel simulator$^{77}$ (LAMMPS) package, using the *brownian* integrator$^{78,79}$. Eqs. (18)-(19) are integrated using a timestep $\delta t = 1$ fs, at a temperature $T = 298$ K. Each simulation consists of 1 ns of equilibration followed by 50 ns of production used for data collection, with observables sampled every 40 fs. As a compromise between computational efficiency and ability to probe the limit of small wave vectors, we use a non-cubic sim-
ulation box with dimensions $L_x \times L_y \times L_z = 500 \times 50 \times 50$ Å³, with periodic boundary conditions in all directions. The smallest accessible wave vector is thus $2\pi/L_x \approx 0.0126$ Å⁻¹, which corresponds to $q_{min}\sigma \approx 0.039$. The results presented in the following sections are averaged over two independent initial conditions and the uncertainties are reported as the standard error computed across the independent trajectories. See Appendix A for additional computational details on the computation of the dynamic structure factor.

## V. RESULTS AND DISCUSSION

### A. Static structure factors and relaxation times

We first compute the longitudinal and static structure factors, whose definition is recalled: $S_{L,T}(q) \equiv \frac{1}{N}\langle|\vec{P}_{L,T}(q)|^2\rangle$. Their $q \to 0$ limits (calculated at the smallest considered wavevector $q_{min}\sigma \approx 0.039$) are plotted respectively on Figs. 1(a) and (b), as a function of the reduced dipole moment $p^*$. The BD results for the Stockmayer fluid (green squares) are compared with the prediction from SDFT (i.e. $q \to 0$ limit of Eqs. (34)-(35), blue circles), obtained by neglecting non-electrostatic interactions and treating electrostatic couplings at the mean-field level. While SDFT predictions for the $q \to 0$ limit of the longitudinal polarization structure factor, $S_L(q \to 0) = p^2/3(1+3y)$, are in excellent agreement with the BD results over the whole dipole range considered, those for the transverse one $S_T(q \to 0) = p^2/3$, deviate from the BD results except for small dipoles.

Identical conclusions can be drawn on the $q \to 0$ limit of the corresponding relaxation times. The latter are determined in BD simulations by fitting the dynamic structure factor $S_{L,T}(q_{min}, \omega)$ as $S_{L,T}(q_{min}) \times 2\tau_{L,T}(q_{min})/(1+\omega^2\tau_{L,T}(q_{min})^2)$, where the value of the static structure factor, discussed above, is not considered as a fitting parameter (see Appendix A for the actual form taking into account the finite sampling frequency in BD simulations). The BD results for the longitudinal and transverse relaxation times are shown as green symbols in Figs. 1(c) and (d), together with the corresponding SDFT predictions $\tau_L(q \to 0) = \tau_s^r/(1+3y)$ and $\tau_L(q \to 0) = \tau_s^r$ (see Eqs. (26) and (27), blue symbols). As for the static structure factor, while SDFT predictions for the $q \to 0$ limit of the longitudinal relaxation time are in excellent agreement with the BD results over the whole dipole range considered, those for the transverse one deviate from the BD results except for small dipoles.

### B. Kirkwood factor and permittivity

Since the above predictions were obtained by neglecting non-electrostatic interactions and treating electrostatic couplings at the mean-field level, we can anticipate that at least some of their shortcomings are due to the fact that short-range correlations are not properly taken into account. Such correlations between dipole are usually characterized by the Kirkwood factor $g_K$ defined in Eq. (5) and shown as a function of

<figure><img src="image_2.png" /><figcaption>FIG. 1. Longitudinal (a) and transverse (b) polarization density structure factors and longitudinal (c) and transverse (d) relaxation times, calculated at the smallest considered wavevector $q_{min}\sigma \approx 0.039$, as a function of the reduced dipole moment $p^* \equiv p\sqrt{\beta/\epsilon_0\sigma^3}$. The structure factors and relaxation times are shown normalized by $p^2$ and the rotational diffusion time of a single dipole, $\tau_s^r = 1/2D_s^r$, respectively. Green squares indicate the results from Brownian Dynamics simulations of the Stockmayer fluid, while blue circles are the SDFT predictions. ‘SDFT’ and ‘effective SDFT’ correspond respectively to the predictions obtained without and with considering the effect of dipolar correlations, respectively. In other words, red stars are the SDFT predictions using effective parameters computed via Eqs. (47), (49), (50) and (51), using the Kirkwood factor $g_K$ measured in BD simulations.</figcaption></figure>

the reduced dipole $p^*$ in Fig. 2. Starting from the expected $g_K = 1$ (no correlations) in the limit $p^* \to 0$, it grows significantly over the considered range of dipoles. As shown in the inset of Fig. 2, these correlations have a direct impact on the static permittivity $\epsilon_r = 1+3yg_K$ (see Eq. (12)), which rises up to $\epsilon_r = 146$, which is close to the value reported in Ref. 26 for similar simulation parameters as the ones used in this study and much larger than the value predicted by assuming $g_K = 1$ (see Eq. (38), also shown as dashed line).

Overall, these observations support the claim that the limitations of SDFT to predict the transverse static structure factor and correlation time are related to the neglect of short-range correlations. The effect of such correlations is less pronounced on the longitudinal static structure factor. This can be understood by taking the limit $q \to 0$, $\omega \to 0$ of Eq. (9), which yields $\chi_L(q=0, \omega) = \frac{C_s}{\epsilon_0 k_B T} S_L(q=0) = 1 - \frac{1}{\epsilon_r}$, where the second equality originates from Eq. (7). Using the relation between $g_K$ and $\epsilon_r$ [Eq. (12)], and with a similar reasoning for the transverse component, we get

$$S_L(q \to 0) = \frac{S_T(q \to 0)}{\epsilon_r} = \frac{p^2 g_K}{3\epsilon_r} = \frac{p^2 g_K}{3(1+3yg_K)} \quad (46)$$

While both components are well described by SDFT Eqs. (34)-(35) when $g_K \approx 1$ (i.e. sufficiently small $p^*$), for
<figure><img src="image_3.png" /><figcaption>FIG. 2. Kirkwood factor $g_K$ as a function of reduced dipole moment $p^*$, obtained from BD simulations (see Eq. (5)). The horizontal dashed line corresponds to $g_K = 1$. The inset shows the corresponding permittivity computed by Eq. (12) (full symbols), and the prediction assuming $g_K = 1$ (empty symbols).</figcaption></figure>

$g_K \gg 1$ there is a compensation between the numerator and denominator for $S_L(q) \approx p^2/9y$ (as correctly predicted by SDFT), whereas no such compensation occurs for $S_T(q)$. A similar compensation might be at play to rationalize the fact that the longitudinal relaxation time is well reproduced even when dipolar correlations are important, unlike the transverse one. However, to the best of our knowledge, there are no exact relations analogous to Eq. 46 for the relaxation times (see nevertheless the discussion below of the approximate link between collective and individual dynamics, which would provide a similar argument).

### C. Accounting for short-range correlations: effective SDFT

It should be no surprise that the approximations leading to the analytical predictions of Section III B fail to capture short-range effects. In fact, one should expect these predictions to be limited to small wavevectors and frequencies. In the following, we show how to introduce the effect of short-range correlations into the above theory via effective quantities (indicated by a tilde), in order to describe polarization fluctuations in the hydrodynamic regime $q \to 0$, $\omega \to 0$ without introducing additional complexity to the model.

To that end, we first note that in order to correctly describe the static correlations $S_T(q \to 0)$ by the SDFT prediction, which corresponds to $g_K^{SDFT} = 1$ (see Eq. (39)), with an effective dipole $\tilde{p}$, Eq. (46) implies $\tilde{p}^2 \equiv p^2 g_K$, i.e.

$$ \tilde{p} \equiv p \sqrt{g_K} = p g_K^{1/2} \quad (47) $$

where $g_K$ is the actual Kirkwood factor of the Stockmayer fluid. Previous works have already introduced an effective dipole to account for the effect of short-range correlations on the electrostatic potential or the pair distribution functions at long distance (see e.g. Refs. 80–82), generally leading to an

increase of the “bare” dipole by a factor $g_K$ (the effective interaction between two distant dipoles being further screened by a factor $\epsilon_r$). In the present case, the square root can be understood by considering the fact that the description with effective dipoles aims at reproducing collective fluctuations and not the actual structure around a given dipole. In order to also correctly describe $S_L(q \to 0)$, one should also map the permittivity at the microscopic and effective levels of description, i.e. $\epsilon_r = 1 + 3y g_K = 1 + 3\tilde{y}$, leading to:

$$ \tilde{y} \equiv y g_K, \quad (48) $$

which together with Eq. (47) implies

$$ \tilde{C}_s \equiv C_s = C_s g_K^0. \quad (49) $$

Kournopoulos et al. also proposed recently a scheme to describe the static dielectric constant of dipolar fluids from molecular considerations, by introducing a numerical scaling factor in the Kirkwood-Onsager relation to improve the agreement between a mean-field theory and the static permittivity measured in MD simulations of the Stockmayer fluid$^{83}$. Their augmented mean-field theory for the static dielectric constant of dipolar hard spheres showed better quantitative agreement with the simulation results when the product $y g_K$ was scaled by a factor of 1.3. In contrast, here we introduce a mapping between the two levels of theory (microscopic $\to$ coarse-grained) using the Kirkwood g-factor obtained from simulations, in order to account for the orientational correlations in numerical applications.

Turning now to dynamical properties, in order to correctly describe the relaxation times $\tau_{L,T}(q \to 0)$ (see Eqs. (14)) one should consider the effect of short-range correlations on the collective dynamics of dipoles. This issue has been investigated in particular by Kivelson and Madden$^{65,84,85}$, who showed that by neglecting some dynamical correlations, one obtains a particularly simple relation between the Debye relaxation time $\tau_D = \tau_T(q \to 0) = \epsilon_r \tau_L(q \to 0)$ and the correlation time for the reorientation dynamics of a single dipole (in the presence of the others), $\tau_1$, namely $\tau_D = \tau_1 g_K$. The SDFT result Eq. (41) suggests that one should rescale the individual rotational time $\tilde{\tau}_s^r \equiv \tau_s^r g_K$, or equivalently the rotational diffusion time:

$$ \tilde{D}_s^r \equiv D_s g_K^{-1}. \quad (50) $$

Finally, in order to consistently describe the $q$-dependence of the relaxation times with SDFT, one should keep $a$ unchanged (see Eqs. (26) and (27)) under rescaling of the rotational diffusion coefficient. This implies the same rescaling of the translational diffusion coefficient,

$$ \tilde{D}_s \equiv D_s g_K^{-1}. \quad (51) $$

In the following, we will refer to ‘effective SDFT’ when using the SDFT predictions with the effective dipole $\tilde{p}$, concentration $\tilde{C}_s$, rotational diffusion coefficient $\tilde{D}_s^r$ and translational diffusion coefficient $\tilde{D}_s$, obtained by rescaling the corresponding bare quantities with the Kirkwood factor determined from simulations using the appropriate exponents presented in
Eqs. (47), (49), (50) and (51). To some extent, this effective way to incorporate dipolar correlations through rescaling of the relevant parameters of the model can be seen as analogous to the dressed ion theory for electrolytes, which consists in a self-consistent renormalization of the charge of ions in order to simplify the analytical description of their static correlations$^{86,87}$.

The predictions of this effective DFT for the static struc-
ture factors and relaxation times, also shown in Fig. 1, are in
excellent agreement with the BD results. This suggests that
incorporating the short-range correlations into SDFT by ap-
propriately rescaling the relevant parameters is sufficient to
capture these key observable properties in the hydrodynamic
regime $q \to 0$, $\omega \to 0$. Before turning to a more stringent test
of effective SDFT with the dynamic structure factors in Sec-
tion V E, we first consider the static structure factor beyond
the $q \to 0$ limit.

D. q-dependence of the static structure factors

The SDFT results obtained by neglecting short-range corre-
lations and treating electrostatic interactions at the mean-field
level predict that the static polarization structure factors (cor-
responding to *t* = 0 in Eqs. (34) and (35)) are independent of
the wavevector *q*. Fig. 3 shows that the BD results for a wide
range of dipoles generally satisfy this prediction for wave vec-
tors *qσ* ≲ 1, *i.e.* down to wavelengths comparable to the par-
ticle size. In addition, the magnitude of the two components
in the *q* → 0 limit satisfy the exact results Eq. (46), with the
Kirkwood factors *g*<sub>K</sub> and permittivity *ε*<sub>r</sub> reported in Fig. 2.

However, we note that for the largest considered dipoles,
the transverse components does depend on the wave vector,
in contrast with the simple SDFT prediction. The BD results
in that case are well described by $S_T(q) = p^2 g_K / (3(1 + q^2 \xi_T^2))$
where $\xi_T$ is a correlation length used as a fitting parameter.
The resulting values are $\xi_T \approx 0.603\sigma$ and $\xi_T \approx 1.28\sigma$ for
$p^* = 5.53$ and 6.15, respectively. These values are compa-
rable to the particle size, consistently with Kirkwood's point
that the dipolar correlations remain relatively short-range (see
e.g. Refs. 63,65,88). In order to capture these correlations,
which are also reflected in the Kirkwood factor, one should
go beyond the assumptions leading to the simple analytical
results of the present work, namely the dipolar approxima-
tion and the omission of any short-range steric repulsion in
the SDFT description.

E. Dynamic polarization structure factors

We finally turn to the polarization dynamics, quantified by
the dynamic structure factors S<sub>L,T</sub>(q, ω). The BD simulation
results for the smallest considered wave vector q<sub>min</sub> σ ≈ 0.039
are shown in Fig. 4. They are very well described by the ef-
fect SDFT predictions, i.e. Eqs. (36)-(37) with effective pa-
rameters (see Eqs. (47), (49), (50) and (51)), also shown as
dashed lines, over the whole range of considered frequencies.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Longitudinal (a) and transverse (b) polarization structure factors (see Eq. (4)), normalized by $p^2$ and as a function of reduced wave number $q\sigma$, for Stockmayer fluids with reduced density $C_s^* = C_s\sigma^3 = 0.9$ and various reduced dipoles $p^* \equiv p\sqrt{\beta/\epsilon_0\sigma^3}$. The dashed lines in panels (a) and (b) correspond to $S_L(q)/p^2 = g_K/3\epsilon_r$ and $S_T(q)/p^2 = g_K/3(1+q^2\xi_T^2)$, with $g_K$ and $\epsilon_r$ from Fig. 2, and $\xi_T$ a correlation length used as a fitting parameter. Colors encode dipole moments in the range $p^* \in [0.5, 6.15]$ as indicated by the colorbar.</figcaption></figure>

The plateau deviating from the $\omega^{-2}$ scaling at the largest fre-
quencies is due to the finite sampling frequency in simulations
and can be captured analytically by considering the discrete
Fourier transform (see Appendix A) of the ISF Eqs. (34) and
(35). Overall, the results of Fig. 4 show that $S_{L,T}(q \to 0, \omega)$
for the overdamped Stockmayer fluid correspond to exponen-
tially decaying intermediate scattering functions of the polar-
ization at all times, as predicted by SDFT. Furthermore, they
confirm that the simple rescaling of the physical parameters by
the Kirkwood factor with appropriate exponents is sufficient
to capture the effect of short-range correlations on the dynam-
ical polarization fluctuations. The SDFT approach therefore
extends well beyond the static limit, and correctly captures
the behavior of the dynamics of the system for a large range
of frequencies, when the correlations between dipoles are ac-
counted for in an effective way.

Finally, we investigate the ability of effective SDFT to cap-
ture the BD results beyond the $q \to 0$ limit, as well as the
effect of the translational and rotational diffusion coefficients
on the dynamics of the polarization fluctuations. Fig. 5 shows
<figure><img src="image_5.png" /><figcaption>FIG. 4. Longitudinal (a) and transverse (b) dynamic polarization structure factors (see Eq. (6)), calculated at the smallest considered wave vector $q_{min} \sigma \approx 0.039$, for Stockmayer fluids with reduced density $C_s^* = C_s \sigma^3 = 0.9$ and various reduced dipoles $p^* \equiv p \sqrt{\beta/\epsilon_0 \sigma^3}$ indicated by the color bar. Results are normalized by $p^2 \tau_s^r$ and expressed as a function of reduced frequency $\omega \tau_s^r$. The shaded regions denote the uncertainty estimates corresponding to one standard error for the BD results. Dashed lines indicate the effective SDFT predictions, i.e. Eqs. (36)-(37) with effective parameters (see Eqs. (47), (49), (50) and (51)).</figcaption></figure>

$S_{L,T}(q, \omega)$ for a wave vector $q\sigma = 0.5$ corresponding to distances comparable to a few particle diameters, for the Stockmayer fluid with the largest considered dipole ($p = 1.85$ D, i.e. $p^* = 6.15$). In this system, short-range correlations are significant, with a Kirkwood factor $g_K \approx 12.6$) and a correlation length $\xi_T \approx 1.28\sigma$, so that the static structure factor $S_T(q = 0.5\sigma^{-1})$ differs significantly from $S_T(q \to 0)$ (see Fig. 3). Starting from a reference system with $D_s = 2.3 \times 10^{-9} \text{ m}^2\text{s}^{-1}$ and $D_s^r = 0.05 \text{ ps}^{-1}$, we consider the effect of $D_s^r$ and $D_s$ by multiplying or dividing them by a factor of 4. The BD results show that, for the considered wave vector, the effect of $D_s^r$ is significant, while that of $D_s$ is very limited, as expected in the small $q$ regime (in the reference case, $a \approx 0.5\sigma$, and in the modified ones $a \approx 0.25\sigma$ or $\sigma$, leading to $qa < 1$ in all cases). Out of this regime, $D_s$ might also play a role on the polarization dynamics$^{28}$. Increasing $D_s^r$ results in a decrease of the characteristic time (larger crossover frequency) and of the plateau at low frequency (note that the frequency and dy-

<figure><img src="image_6.png" /><figcaption>FIG. 5. Effect of the translational and rotational diffusion coefficients on the longitudinal [(a) and (c)] and transverse [(b) and (d)] dynamic polarization structure factors (see Eq. (6)), calculated at $q_{min} \sigma \approx 0.5$, for the Stockmayer fluid with reduced density $C_s^* = C_s \sigma^3 = 0.9$ and reduced dipoles $p^* = 6.15$. Starting from a reference system (blue) with $D_s = 2.3 \times 10^{-9} \text{ m}^2\text{s}^{-1}$, and $D_s^r = 0.05 \text{ ps}^{-1}$, we consider the effect of $D_s^r$ [(a) and (b)] and $D_s$ [(c) and (d)] by dividing (red) or multiplying (green) them by a factor of 4. The x and y-axes of all plots are normalized by $1/\tau_s^r$ and $p^2 \tau_s^r$, respectively, with the rotational diffusion time of the reference system. The shaded regions denote the uncertainty estimates corresponding to one standard error for the BD results. Dashed lines indicate the effective SDFT predictions, i.e. Eqs. (36)-(37) with effective parameters (see Eqs. (47), (49), (50) and (51)).</figcaption></figure>

namic structure factors in Fig. 5 are adimensionalized using $\tau_s^r$ of the reference system). Importantly, Fig. 5 shows that the BD results and the effects of $D_s^r$ and $D_s$ are well captured by effective SDFT even for this case beyond the $q \to 0$ limit of Fig. 4.

## VI. CONCLUSION AND PERSPECTIVES

In this work, we developed and validated a stochastic density functional theory (SDFT) framework to describe the polarization dynamics of polar fluids, using the Stockmayer fluid as a prototypical model. Analytical expressions were derived for the intermediate scattering functions and dynamic structure factors of longitudinal and transverse polarization components under linearized SDFT assumptions. These predictions were then systematically compared with Brownian Dynamics simulations, which served as a microscopic reference.

Our results demonstrate that SDFT, which neglects short-range correlations and treats dipole-dipole interactions at the mean-field level in its current form, provides accurate predictions for the longitudinal polarization fluctuations. However, it underestimates the amplitude and relaxation time of transverse fluctuations, especially for fluids with strong dipo-
lar interactions. We show that this discrepancy arises from the neglect of local orientational correlations, which can be effectively captured by introducing the Kirkwood factor into a modified SDFT. This effective theory, based on parameter rescaling guided by simulation data, achieves quantitative agreement with both static and dynamic properties in the long-wavelength limit.

Beyond validating the theoretical framework, this work provides practical guidelines for the use of SDFT to model real polar solvents and their role in complex environments, and several future research directions. (i) In particular, the present SDFT theory sets the stage for investigating fluctuation-induced phenomena in electrolytes,⁸⁹⁻⁹² where polar solvent dynamics could play a critical role, and has been overlooked so far. This could involve the integration of SDFT description of solvent into the fluctuating hydrodynamics model for electrolyte systems⁹³,⁹⁴, thereby providing a more unified theoretical framework that simultaneously accounts for the hydrodynamic interactions and solvent polarization. (ii) Furthermore, the model can be extended to account for spatial confinement — an essential factor in many nanoscale systems, including ionic channels, membranes, and porous materials. Previous studies have already shown that confinement induces anisotropy in the static dielectric properties of polar fluids⁹⁵⁻⁹⁷. Further investigation into the confinement effects on the polarization correlations and collective relaxation dynamics could substantially enhance our understanding of nanoscale transport processes. To conclude, we assert that by providing a consistent and computationally tractable approach, SDFT thus offers a promising path toward bridging molecular-scale interactions with mesoscopic behavior in a wide range of soft and electrochemical systems.

## ACKNOWLEDGMENTS

This project received funding from the European Research Council under the European Union's Horizon 2020 research and innovation program (grant agreement no. 863473). P.I. acknowledges financial support from Agence Nationale de la Recherche through project TraNonEq (ANR-24-CE30-0651).

## AUTHOR DECLARATIONS

### Conflict of interest

There is no conflict of interest to declare.

### Author contributions

**Sleeba Varghese:** Conceptualization (supporting); Formal analysis (equal); Investigation (lead); Methodology (supporting); Writing/Original Draft Preparation (lead); Writing – review & editing (equal). **Pierre Illien:** Conceptualization (equal); Formal analysis (equal); Investigation (supporting);

Methodology (equal); Supervision (supporting); Writing – review & editing (equal). **Benjamin Rotenberg:** Conceptualization (equal); Formal analysis (equal); Funding Acquisition (lead); Investigation (supporting); Methodology (equal); Supervision (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The original data presented in this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.17061930

## Appendix A: Additional computational details

The partial Fourier components of the longitudinal and transverse polarization densities, $\tilde{P}_{L,T}(q,t)$, are sampled every 40 fs during the production run of the simulation, for selected wave vectors in the range $q\sigma \in [0.039, 11]$, satisfying $q = m2\pi/L_x$ with $m \in \mathbb{N}$ to ensure compatibility with periodic boundary conditions. The dynamic polarization structure factor is computed for each wave vector using the Fast Fourier Transform algorithm, applied over $n$ blocks of the time series assumed to be statistically independent ($n = 50$ corresponding to 1 ns for the longitudinal polarization density, $n = 10$ corresponding to 5 ns for the transverse one). The results are further averaged over the two independent initial conditions, with uncertainties estimated as the standard error across the independent trajectories.

To incorporate the effect of finite sampling frequency into the analytical predictions, we express the dynamic polarization structure factors using the discrete Fourier transform of Eqs. (34)-(35). Hence, Eqs. (36)-(37) become

$$S_{L,T}(q, \omega_k) = \frac{2S_{L,T}(q)}{\tau_{L,T}(q)} \left| \frac{1}{1 - r(\omega_k)} \right|^2 (\Delta t)^2, \quad (A1)$$

where for each discrete frequency $\omega_k = 2\pi k/N_b\Delta t$, with $\Delta t = 40\delta t$ the interval between successive samples, $N_b$ the total number of samples in a block and $k \in \{0, 1, 2, ..., N_b\}$, we introduced

$$r(\omega_k) = \exp\left[-\left(\frac{\Delta t}{\tau_{L,T}} + i\omega_k\Delta t\right)\right]. \quad (A2)$$

<footnote n="1">¹C. Gray and K. Gubbins, *Theory of Molecular Liquids* (Oxford University Press, 1984).</footnote>

<footnote n="2">²P. J. Rossky, Annu. Rev. Phys. Chem. **36**, 321 (1985).</footnote>

<footnote n="3">³J. P. Hansen and I. R. McDonald, *Theory of Simple Liquids*, 4th ed. (Elsevier, Amsterdam, 2013).</footnote>

<footnote n="4">⁴G. Jeanmairet, B. Rotenberg, and M. Salanne, Chem. Rev. **122**, 10860 (2022).</footnote>

<footnote n="5">⁵B. Bagchi and A. Chandra, in *Adv. Chem. Phys.*, Vol. 80, edited by I. Prikogine and S. A. Rice (Wiley, 1991) 1st ed., pp. 1–126.</footnote>

<footnote n="6">⁶N. Nandi, K. Bhattacharyya, and B. Bagchi, Chem. Rev. **100**, 2013 (2000).</footnote>

<footnote n="7">⁷B. Bagchi and B. Jana, Chem. Soc. Rev. **39**, 1936 (2010).</footnote>

<footnote n="8">⁸E. L. Pollock and B. J. Alder, Physica A **102A**, 1 (1980).</footnote>

<footnote n="9">⁹M. Neumann, Mol. Phys. **50**, 841 (1983).</footnote>

<footnote n="10">¹⁰C. Hesse-Bezot, G. Bossis, and C. Brot, J. Chem. Phys. **80**, 3399 (1984).</footnote>
<footnote n="11">¹¹M. Neumann, Mol. Phys. **57**, 97 (1986).</footnote>

<footnote n="12">¹²M. Neumann, J. Chem. Phys. **85**, 1567 (1986).</footnote>

<footnote n="13">¹³P. A. Bopp, A. A. Kornyshev, and G. Sutmann, J. Chem. Phys. **109**, 1939 (1998).</footnote>

<footnote n="14">¹⁴J. Bartke and R. Hentschke, Phys. Rev. E **75**, 061503 (2007).</footnote>

<footnote n="15">¹⁵F. H. Stillinger and A. Rahman, J. Chem. Phys. **60**, 1545 (1974).</footnote>

<footnote n="16">¹⁶H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. **91**, 6269 (1987).</footnote>

<footnote n="17">¹⁷M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. **112**, 8910 (2000).</footnote>

<footnote n="18">¹⁸S. W. Rick, J. Chem. Phys. **120**, 6085 (2004).</footnote>

<footnote n="19">¹⁹J. L. F. Abascal and C. Vega, J. Chem. Phys. **123**, 234505 (2005).</footnote>

<footnote n="20">²⁰H. Berthoumieux, J. Chem. Phys. **148**, 104504 (2018).</footnote>

<footnote n="21">²¹H. Berthoumieux and F. Paillusson, J. Chem. Phys. **150**, 094507 (2019).</footnote>

<footnote n="22">²²H. Berthoumieux and A. C. Maggs, J. Chem. Phys. **143**, 104501 (2015).</footnote>

<footnote n="23">²³P. Frodl and S. Dietrich, Phys. Rev. A **45**, 7330 (1992).</footnote>

<footnote n="24">²⁴P. Frodl and S. Dietrich, Phys. Rev. E **48**, 3741 (1993).</footnote>

<footnote n="25">²⁵G. Jeanmairet, M. Levesque, R. Vuilleumier, and D. Borgis, J. Phys. Chem. Letters **4**, 619 (2013).</footnote>

<footnote n="26">²⁶G. Jeanmairet, N. Levy, M. Levesque, and D. Borgis, J. Phys.: Condens. Matter **28**, 244005 (2016).</footnote>

<footnote n="27">²⁷B. Bagchi and A. Chandra, in *Proc. Indian Acad. Sci. Chem. Sci.*, Vol. 100 (Springer, 1988) pp. 353–357.</footnote>

<footnote n="28">²⁸A. Chandra and B. Bagchi, Chem. Phys. Lett. **151**, 47 (1988).</footnote>

<footnote n="29">²⁹A. Chandra and B. Bagchi, J. Chem. Phys. **91**, 1829 (1989).</footnote>

<footnote n="30">³⁰D. F. Calef and P. G. Wolynes, J. Chem. Phys. **78**, 4145 (1983).</footnote>

<footnote n="31">³¹A. Chandra and B. Bagchi, J. Phys. Chem. **93**, 6996 (1989).</footnote>

<footnote n="32">³²G. Vijayadamodar, A. Chandra, and B. Bagchi, Chem. Phys. Lett. **161**, 413 (1989).</footnote>

<footnote n="33">³³G. Vijayadamodar and B. Bagchi, J. Chem. Phys. **95**, 1168 (1991).</footnote>

<footnote n="34">³⁴S. Roy and B. Bagchi, J. Chem. Phys. **99**, 1310 (1993).</footnote>

<footnote n="35">³⁵S. Roy and B. Bagchi, J. Chem. Phys. **99**, 9938 (1993).</footnote>

<footnote n="36">³⁶M. te Vrugt, H. Löwen, and R. Wittkowski, Adv. Phys. **69**, 121 (2020).</footnote>

<footnote n="37">³⁷A. Chandra and B. Bagchi, J. Phys. Chem. **94**, 3152 (1990).</footnote>

<footnote n="38">³⁸M.-T. Hoang Ngoc, J. Kim, G. Pireddu, I. Chubak, S. Nair, and B. Rotenberg, Faraday Discuss. **246**, 198 (2023).</footnote>

<footnote n="39">³⁹A. D. Young, A. L. Thornevwork, and S. Marbach, “Decoding noise in nanofluidic systems: Adsorption versus diffusion signatures in power spectra,” (2025), arXiv:2506.23867.</footnote>

<footnote n="40">⁴⁰T. Munakata, J. Phys. Soc. Jpn. **58**, 2434 (1989).</footnote>

<footnote n="41">⁴¹K. Kawasaki, Physica A **208**, 35 (1994).</footnote>

<footnote n="42">⁴²K. Kawasaki, J. Stat. Phys. **93**, 527 (1998).</footnote>

<footnote n="43">⁴³D. S. Dean, J. Phys. A: Math. Gen. **29**, L613 (1996).</footnote>

<footnote n="44">⁴⁴P. Illien, “The Dean-Kawasaki equation and stochastic density functional theory. arXiv:2411.13467,” (2024).</footnote>

<footnote n="45">⁴⁵V. Démery and D. S. Dean, J. Stat. Mech. **2016**, 023106 (2016).</footnote>

<footnote n="46">⁴⁶J.-P. Péraud, A. Nonaka, A. Chaudhri, J. B. Bell, A. Donev, and A. L. Garcia, Phys. Rev. Fluids **1**, 074103 (2016).</footnote>

<footnote n="47">⁴⁷A. Donev, A. L. Garcia, J. P. Péraud, A. J. Nonaka, and J. B. Bell, Curr. Opin. Electrochem. **13**, 1 (2019).</footnote>

<footnote n="48">⁴⁸Y. Avni, D. Andelman, and H. Orland, J. Chem. Phys. **157**, 154502 (2022).</footnote>

<footnote n="49">⁴⁹Y. Avni, R. M. Adar, D. Andelman, and H. Orland, Phys. Rev. Lett. **128**, 098002 (2022).</footnote>

<footnote n="50">⁵⁰O. Bernard, M. Jardat, B. Rotenberg, and P. Illien, J. Chem. Phys. **159**, 164105 (2023).</footnote>

<footnote n="51">⁵¹H. Bonneau, V. Démery, and É. Raphaël, J. Stat. Mech. **2023**, 073205 (2023).</footnote>

<footnote n="52">⁵²H. Berthoumieux, V. Démery, and A. C. Maggs, J. Chem. Phys. **161**, 184504 (2024).</footnote>

<footnote n="53">⁵³H. Frusawa, Entropy **22**, 34 (2019).</footnote>

<footnote n="54">⁵⁴H. Frusawa, Soft Matter **18**, 4280 (2022).</footnote>

<footnote n="55">⁵⁵H. Wada, J. Stat. Mech. **2005**, P01001 (2005).</footnote>

<footnote n="56">⁵⁶R. Okamoto, J. Stat. Mech. **2022**, 093203 (2022).</footnote>

<footnote n="57">⁵⁷P. Robin, J. Chem. Phys. **160**, 064503 (2024).</footnote>

<footnote n="58">⁵⁸L. F. Cugliandolo, P.-M. Déjardin, G. S. Lozano, and F. Van Wijland, Phys. Rev. E **91**, 032139 (2015).</footnote>

<footnote n="59">⁵⁹P. M. Déjardin, Y. Cornaton, P. Ghesquière, C. Caliot, and R. Brouzet, J. Chem. Phys. **148**, 044504 (2018).</footnote>

<footnote n="60">⁶⁰P.-M. Déjardin, S. V. Titov, and Y. Cornaton, Phys. Rev. B **99**, 024304 (2019).</footnote>

<footnote n="61">⁶¹P.-M. Déjardin and Y. Cornaton, J. Phys.: Conf. Ser. **1322**, 012039 (2019).</footnote>

<footnote n="62">⁶²P. Illien, A. Carof, and B. Rotenberg, Phys. Rev. Lett. **133**, 268002 (2024).</footnote>

<footnote n="63">⁶³J. G. Kirkwood, J. Chem. Phys. **7**, 911 (1939).</footnote>

<footnote n="64">⁶⁴R. L. Fulton, Mol. Phys. **29**, 405 (1975).</footnote>

<footnote n="65">⁶⁵P. Madden and D. Kivelson, in *Adv. Chem. Phys.* (John Wiley & Sons, Ltd, 1984) pp. 467–566.</footnote>

<footnote n="66">⁶⁶B. U. Felderhof, Physica A **101**, 275 (1980).</footnote>

<footnote n="67">⁶⁷E. L. Pollock and B. J. Alder, Phys. Rev. Lett. **46**, 950 (1981).</footnote>

<footnote n="68">⁶⁸J. M. Caillol, Europhys. Lett. **4**, 159 (1987).</footnote>

<footnote n="69">⁶⁹P. V. Giaquinta, M. Parrinello, and M. P. Tosi, Physica A **92**, 185 (1978).</footnote>

<footnote n="70">⁷⁰B. M. Ladanyi and B.-C. Perng, AIP Conf. Proc. **492**, 250 (1999).</footnote>

<footnote n="71">⁷¹S. W. de Leeuw, J. W. Perram, E. R. Smith, and J. S. Rowlinson, Proc. R. Soc. London, Ser. A **373**, 27 (1980).</footnote>

<footnote n="72">⁷²J. M. Caillol, J. Chem. Phys. **96**, 7039 (1992).</footnote>

<footnote n="73">⁷³D. Kivelson and H. Friedman, J. Phys. Chem. **93**, 7026 (1989).</footnote>

<footnote n="74">⁷⁴G. D. Harp and B. J. Berne, Phys. Rev. A **2**, 975 (1970).</footnote>

<footnote n="75">⁷⁵R. Ramirez, R. Gebauer, M. Mareschal, and D. Borgis, Phys. Rev. E **66**, 031206 (2002).</footnote>

<footnote n="76">⁷⁶J. J. Cerdà, V. Ballenegger, O. Lenz, and C. Holm, J. Chem. Phys. **129** (2008).</footnote>

<footnote n="77">⁷⁷A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in 't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, Comput. Phys. Commun. **271**, 108171 (2022).</footnote>

<footnote n="78">⁷⁸S. Delong, F. B. Usabiaga, and A. Donev, J. Chem. Phys. **143**, 144107 (2015).</footnote>

<footnote n="79">⁷⁹I. M. Ilie, W. J. Briels, and W. K. Den Otter, J. Chem. Phys. **142**, 114103 (2015).</footnote>

<footnote n="80">⁸⁰G. Nienhuis and J. M. Deutch, J. Chem. Phys. **56**, 5511 (1972).</footnote>

<footnote n="81">⁸¹R. Finken, V. Ballenegger, and J.-P. Hansen, Mol. Phys. **101**, 2559 (2003).</footnote>

<footnote n="82">⁸²V. Ballenegger and J.-P. Hansen, Mol. Phys. **102**, 599 (2004).</footnote>

<footnote n="83">⁸³S. Kournopoulos, A. J. Haslam, G. Jackson, A. Galindo, and M. Schoen, J. Chem. Phys. **156**, 154111 (2022).</footnote>

<footnote n="84">⁸⁴D. Kivelson and P. Madden, Mol. Phys. **30**, 1749 (1975).</footnote>

<footnote n="85">⁸⁵T. Samanta and D. V. Matyushov, J. Mol. Liq. **364**, 119935 (2022).</footnote>

<footnote n="86">⁸⁶R. Kjellander and D. J. Mitchell, J. Chem. Phys. **101**, 603 (1994).</footnote>

<footnote n="87">⁸⁷R. Kjellander, J. Chem. Phys. **145**, 124503 (2016).</footnote>

<footnote n="88">⁸⁸C. Zhang, J. Hutter, and M. Sprik, J. Phys. Chem. Lett. **7**, 2696 (2016).</footnote>

<footnote n="89">⁸⁹S. Mahdisoltani and R. Golestanian, Phys. Rev. Lett. **126**, 158002 (2021).</footnote>

<footnote n="90">⁹⁰S. Mahdisoltani and R. Golestanian, New J. Phys. **23**, 073034 (2021).</footnote>

<footnote n="91">⁹¹G. Du, D. S. Dean, B. Miao, and R. Podgornik, Phys. Rev. Lett. **133**, 238002 (2024).</footnote>

<footnote n="92">⁹²G. Du, D. S. Dean, B. Miao, and R. Podgornik, Phys. Rev. E **111**, 044108 (2025).</footnote>

<footnote n="93">⁹³J.-P. Péraud, A. Nonaka, A. Chaudhri, J. B. Bell, A. Donev, and A. L. Garcia, Phys. Rev. Fluids **1**, 074103 (2016).</footnote>

<footnote n="94">⁹⁴A. Donev, A. L. Garcia, J.-P. Péraud, A. J. Nonaka, and J. B. Bell, Curr. Opin. Electrochem. **13**, 1 (2019).</footnote>

<footnote n="95">⁹⁵I. Szalai and S. Dietrich, Eur. Phys. J. E **28**, 347 (2009).</footnote>

<footnote n="96">⁹⁶A. Schlaich, E. W. Knapp, and R. R. Netz, Phys. Rev. Lett. **117**, 048001 (2016).</footnote>

<footnote n="97">⁹⁷L. Fumagalli, A. Esfandiar, R. Fabregas, S. Hu, P. Ares, A. Janardanan, Q. Yang, B. Radha, T. Taniguchi, K. Watanabe, et al., *Science* **360**, 1339 (2018).</footnote>