Memory effects in fluctuating dynamic density-functional theory: theory and simulations

Many applications, from chemistry to biology and engineering, exploit complex fluids consisting of colloidal particles. A number of challenging issues in applications requires chemically reacting particles while colloidal chemistry is often the preferred path for particulate systems and formulated products. Examples range from the design of micro-/nanofluidic devices [1], lab-on-chip applications [2] and analysis of chemical reactors [3] to biological macromolecules in solutions [4], dispersion of carbon nanotubes in solvents [5] and preparation–formulation of organic nanoparticles [6].

Physical properties at the nanoscale may be investigated in a 'brute-force' manner by solving the Newton's equations of motion for each individual particle, i.e. by means of molecular dynamics or Monte Carlo simulations. However, many interesting dynamical events, such as protein folding [7, 8] and phase nucleation [9, 10], occur at time- and lengthscales several orders of magnitude larger than those accessible by atomistic simulations [11]. Due also to the prohibitively high computational cost associated with such simulations, coarse-graining remains a highly relevant instrument in the study of reacting systems. Precursory model reduction techniques go back to the phenomenological description of pollen particles in water by Brown [12]. However, it was not until the work of Langevin [13] that a low-dimensional (stochastic) model for such phenomena was formulated. Since then, several reduced models have been proposed [14–26]. A standard way to coarse-grain a system is by using projection-operator techniques to decouple the relevant and the non-relevant degrees of freedom in the system [16, 20, 25, 27]. This leads to a system of stochastic generalized Langevin equations (GLEs) for the time evolution of the resolved variables (relevant degrees of freedom). The effects of the unresolved variables (the irrelevant degrees of freedom) are included as the noisy part of the GLEs. Such GLEs typically have non-Markovian terms involving time convolutions which carry information about the history of the dynamical system at hand. Most existing studies focus on the so-called Markovian limit of the full GLEs where the memory effects are ignored [28, 29]. Importantly, the Markovian approximation can be inaccurate when the system does not have a clear scale separation between the resolved and the unresolved variables. It is worth underlying that, in general, GLEs are used to describe the stochastic time-evolution of a set of observables of interest of a dynamical system. These observables can be either some particle coordinates (microscopic observables) or continuum field variables (macroscopic observables). In this work, GLEs are used to obtain the time evolution of some particle coordinates of interest (or resolved particles). It follows that this Langevin approach still gives a description of the system at particle-level, and thus it is essentially equivalent in computational complexity to atomistic simulations.

A continuum/mean-field framework capable of describing small-scale systems is offered by density-functional theory (DFT) or its dynamical counterparts known as dynamical DFT (DDFT), and fluctuating hydrodynamics (FH). Unlike the atomistic GLEs, both DDFT and FH operate with local averaged densities and are thus computationally tractable. DDFT is usually obtained by deriving the Fokker–Planck equation governing the evolution of the system probability density function, which is then averaged over all but one degrees of freedom [22, 30–32]. Recently, several extensions of DDFT were proposed for multi-component systems [30] to include hydrodynamic interactions [31–33] and effects of particle orientability [34]. On the other hand, FH for simple liquids was initially stipulated by Landau et al from empirical arguments [35]. The connection between FH and the density-functional framework for small-scale systems was later highlighted in the works of Kawasaki and Dean (DK) [28, 36] (with Dean's equation being the overdamped counterpart of Kawasaki's). In the DK approach, a heuristic energy functional accounting for small-scale features was introduced. Recently, a rigorous and systematic derivation of FH and fluctuating DDFT (FDDFT) for the general case of arbitrarily shaped thermalized particles was offered by Durán-Olivencia et al [34].

It should be emphasised, however, that the DK model describes the evolution of instantaneous microscopic fields, thus containing the same physical information as the original set of GLEs, and hence leaves the equations computationally intractable and at the same time disconnected from the original Landau et al theory. Even more importantly, this disconnection has often led to the misconception that the DK model describes the evolution of macroscopic observables. The formalisation of the connection between DK and FH to enable the description of the evolution of observable fields and remove the dependence on the microscopic GLEs has lately received a great deal of attention. Indeed, recent efforts have achieved the desired formal connection under the local-equilibrium approximation [29, 37]. Finally, applications of FH to multi-component reacting systems have been studied in many works, e.g. references [38–40].

Thus far, the inclusion of fluctuations in the DDFT framework has posed many questions [41]. For instance, in some studies such as the work by Elder et al [42], an additional noise term was artificially included, but the physical meaning of this term remains questionable. At the same time, the derivation of the Fokker–Planck equation, the first step in the derivation of DDFT, requires that we already averaged out the noise. The recent effort by Durán-Olivencia et al [29] showed that DDFT can be understood as the most-likely realisation of FDDFT, thus providing closure to a long standing debate in the density-functional theory community about the inclusion of fluctuations establishing also rigorously the connection between FDDFT and FH from first principles.

In the present work we start with microscopic GLEs of reacting systems to derive the corresponding FH equations. The underlying GLEs allow us to obtain the proper noise correlations in FH, as well as include memory effects. We derive the non-Markovian FH governing the time evolution of the one-body density and momentum fields in systems of several reacting species. Adopting a density-functional formulation, we transform the FH equation to an equivalent FDDFT equation. In the appropriate overdamped limit, the latter reduces to a single equation for the time-dependent one-body density field.

In general, numerically solving the non-Markovian FH is quite challenging due to the memory term and the time-correlated noise. For a memory term in the form of decaying exponential series, the computational issues are overcome by introducing an additional field, whose time evolution accounts for both the memory term and noise time correlations. Such an approach will be denoted as 'extended field dynamics'. The resulting non-Markovian DDFT equations are then solved numerically using a finite-volume approach. Our numerical scheme is exemplified with several applications to multi-component reacting systems leading to the emergence of stochastic Turing patterns [43–51].

In section 2 we introduce the GLEs for an isothermal multi-component system of particles, which provide the starting point for our derivation. Section 3 deals with the inclusion of microscopic density and momentum fields, and the derivation of the FH equations for a multi-species reacting system. A non-trivial deconvolution operator has to be defined to take the overdamped limit of the previously derived FH. This leads to the formulation of a non-Markovian overdamped FDDFT equation for the density field in section 4. The computational challenges due to both the memory term and the correlated noise are overcome by introducing an extended-field variable dynamics framework in section 5. The finite-volume method is then adopted for the spatiotemporal discretization of the extended variable dynamics in section 6. Section 7 reports several computational experiments including a comparison between molecular dynamics simulations and FH that validates our main theoretical results, and a dedicated study of memory-driven Turing patterns. Conclusions and outlook for possible future work are offered in section 8.

In figure 1, we show a flowchart summarizing the main equations in this work and their connections to previous formulations. Moreover, for the sake of clarity, we include in table 1 a summary of the abbreviations used in the text and their definition.

Zoom In Zoom Out Reset image size

Consider an isolated system containing a total of N_tot,s (t) particles at time t, with s = 1, ..., K denoting the different coexisting species. The phase-space coordinates of a given species s are noted as z_s = {r_s , p_s }, with ${\mathbf{r}}_{s}=\left\{{\mathbf{r}}_{1,s}\dots {\mathbf{r}}_{{N}_{\text{tot},s}}\right\}$ and ${\mathbf{p}}_{s}=\left\{{\mathbf{p}}_{1,s}\dots {\mathbf{p}}_{{N}_{\text{tot},s}}\right\}$ being the particle positions and momenta, respectively. In many cases, such as colloids in suspension, one is interested in the dynamical evolution of only the first N_s of the original N_tot,s particles belonging to the species s. The coordinates of the particles of interest are also known as resolved variables, while the remaining ones are denoted as unresolved variables or thermal bath. Using time-scale separation (or weak coupling) between resolved and unresolved variables, the resolved ones can be shown to obey the following GLE [16, 20, 25, 56]:

$$\begin{equation}\begin{aligned}\hfill {\dot {\boldsymbol{r}}}_{i,s}\left(t\right)& ={m}_{s}^{-1}{\boldsymbol{p}}_{i,s}\left(t\right),\hfill \\ \hfill {\dot {\boldsymbol{p}}}_{i,s}\left(t\right)& ={\boldsymbol{F}}_{i,s}\left(t\right)-{\int }_{{t}_{0}}^{t}\mathrm{d}{t}^{\prime }{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right){\boldsymbol{p}}_{i,s}\left({t}^{\prime }\right)+{\boldsymbol{R}}_{i,s}\left(t\right),\hfill \end{aligned}\end{equation} \tag{ 1 }$$

where ${\boldsymbol{F}}_{i,s}=-{\nabla }_{{\boldsymbol{r}}_{i,s}}U\left({\boldsymbol{r}}_{i,s}\right)$ is the effective force experienced by each particle, with U( r _i,s ) being the potential of mean force in the form:

$$\begin{equation}U\left({\boldsymbol{r}}_{i,s},{\boldsymbol{r}}^{N}\right)=\sum _{\alpha =1,\alpha \ne s}^{K}\sum _{j=1}^{{N}_{\alpha }}{V}_{\alpha }\left({\boldsymbol{r}}_{i,s}-{\boldsymbol{r}}_{j,\alpha }\right)+\sum _{j=1,j\ne i}^{{N}_{s}\left(t\right)}{V}_{s}\left({\boldsymbol{r}}_{i,s}-{\boldsymbol{r}}_{j,s}\right)+{U}_{\text{ext}}\left({\boldsymbol{r}}_{i,s}\right),\end{equation} \tag{ 2 }$$

where V_α is the two-body interacting potential acting between particles of the species s and α, and U_ext is the external potential acting on the system. It is worth highlighting that the underlying assumption in equation (2) is that the inter-particle interactions are governed by one- and two-body interactions. While this assumption allows us to significantly simplify the derivation in what follows, it also poses some limitations. In fact, the energy functional that we derive is only able to account for the two-body correlations. For applications requiring higher-order interactions, the assumption could be relaxed, and our derivation could be easily extended in order to account for any N-body interactions. Should this be the case, the structure of the equations derived hereinafter [specifically equation (26)] would lead to the exact first equation of the YBG hierarchy, which involves all N-body densities, and hence correlations.

The second term on the right-hand side of equation (1) is widely known as the damping force, representing the deterministic drag coming out from the random collisions of the target particles with the surrounding bath particles. Such a drag force involves a memory-kernel tensor θ _s (t) which makes the force dependent on the entire history of the particle dynamics [25, 26, 52–55]. It is worth noticing that if the thermal bath is characterized by multiple timescales, the dynamical effects of the bath particles of the original system may not be fully embedded in the memory kernel θ _s (t) of equation (1). But these applications are beyond the scope of this work.

The stochastic term R _i,s (t) in equation (1) accounts for the volatility which the just mentioned collisions with the bath induces on the particle's acceleration. This emerges as a Gaussian distribution with zero mean and a correlation function satisfying the fluctuation dissipation theorem:

$$\begin{equation}\begin{aligned}\hfill \langle {\boldsymbol{R}}_{i,s}\left(t\right)\rangle & =0,\hfill \\ \hfill \langle {\boldsymbol{R}}_{i,s}\left(t\right){\boldsymbol{R}}_{i,s}\left({t}^{\prime }\right)\rangle & ={m}_{s}{k}_{\text{B}}T{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right),\hfill \end{aligned}\end{equation} \tag{ 3 }$$

where k_B is Boltzmann's constant, T is the absolute temperature imposed by the system of bath particles and the brackets ⟨⋅⟩ denote the ensemble average over multiple realizations of the noise, i.e. ⟨A⟩ = ∫ ρ(z)A(z) dz, with ρ(z) being the probability density function defined in the phase space of the original system. Despite being derived for equilibrium conditions, in this work equation (3) is also applied to systems in local-equilibrium. This extension is rigorously valid if the thermal bath is assumed to relax instantly toward local-equilibrium conditions. For the sake of clarity, we need to highlight that all particles of the same species, say s, have the same mass and memory kernel and interact only via the pairwise potential V_s , as specified in equation (2).

Some systems are characterized by a time-scale separation between resolved and unresolved particle dynamics, due for instance to the much higher inertia of the first compared to the latter particles. In these cases, the memory kernel tensor in equation (1) can be approximated as θ _s (t) = 2 θ _s,0 δ(t), where δ(t) is a Dirac's delta function and θ _s,0 is a constant. This approximation allows one to obtain a memory-less (or Markovian) formulation, with a much simpler structure. In fact, in the Markovian limit the memory term in equation (1) loses its convolution structure, i.e. ${\int }_{{t}_{0}}^{t}\mathrm{d}{t}^{\prime }{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right){\boldsymbol{p}}_{i,s}\left({t}^{\prime }\right)={\boldsymbol{\theta }}_{s,0}{p}_{i,s}\left(t\right)$, while the stochastic term becomes delta-correlated in time, i.e. ⟨ R _i,s (t) R _i,s (t')⟩ = m_s k_B T θ _s,0 δ(t − t') [26]. However, when no clear scale separation between colloidal and bath particles is present, the Markovian hypothesis is not applicable and one needs to seek a way of approximating the memory kernel to numerically track the dynamics [21, 23, 25, 56].

In a recent study on the parameterization of GLEs, Lei et al used a data-driven approach to approximate the memory kernel [56]. Following this study, one could consider atomistic observations of a system at equilibrium, then right-multiply equation (1) by the initial momentum p _i (0) and, finally, take the ensemble average with respect to the equilibrium distribution function, obtaining the following equation:

$$\begin{equation}{\boldsymbol{g}}_{i,s}\left(t\right)={\int }_{0}^{t}{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right){\boldsymbol{h}}_{i,s}\left({t}^{\prime }\right)\enspace \mathrm{d}{t}^{\prime },\end{equation} \tag{ 4 }$$

where we have defined the matrices ${\boldsymbol{g}}_{i,s}\left(t\right)=\langle \left({\boldsymbol{F}}_{i,s}-\dot {{\boldsymbol{p}}_{i,s}}\left(t\right)\right)\cdot {\boldsymbol{p}}_{i,s}^{T}\left(0\right)\rangle$ and ${\boldsymbol{h}}_{i,s}\left(t\right)=\langle {\boldsymbol{p}}_{i,s}\left(t\right)\cdot {\boldsymbol{p}}_{i,s}^{T}\left(0\right)\rangle$ and made use of the orthogonality between the stochastic term and the particle momentum, i.e. ⟨ R _i,s (t) ⋅ p _i,s (0)^T ⟩ = 0 [25]. Taking the Laplace transform ¹ of equation (4), we get

$$\begin{equation}{\tilde {\boldsymbol{\theta }}}_{s}\left(\lambda \right)={\tilde {\boldsymbol{g}}}_{i,s}\left(\lambda \right){\tilde {\boldsymbol{h}}}_{i,s}^{-1}\left(\lambda \right),\end{equation} \tag{ 5 }$$

where ${\tilde {\boldsymbol{g}}}_{i,s}\left(\lambda \right)$, ${\tilde {\boldsymbol{h}}}_{i,s}\left(\lambda \right)$ and ${\tilde {\boldsymbol{\theta }}}_{s}\left(\lambda \right)$ denote the Laplace transforms of h _i,s (t), g _i,s (t) and θ _s (t), respectively. Assuming that the memory kernel θ(t) is isotropic and represented in the form of decaying exponential series, ${\boldsymbol{\theta }}_{s}\left(t\right)={\sum }_{k}{A}_{s,k}{\mathrm{e}}^{{B}_{s,k}t}\mathbf{I}$, we can express its Laplace transform as:

$$\begin{equation}{\tilde {\boldsymbol{\theta }}}_{s}\left(\lambda \right)=\sum _{k}^{d}\frac{{A}_{s,k}\lambda }{1-{B}_{s,k}\lambda }\mathbf{I},\end{equation} \tag{ 6 }$$

where I is the identity matrix with dimension 1, 2 or 3 depending on the dimensionality of the system, A_s,k and B_s,k are real-valued and B_s,k ⩽ 0. In a data-driven approach, if the matrices h _i,s (t) and g _i,s (t) are computed from historical data, the sets of coefficients {A_s,1, ..., A_s,k , ..., A_s,d } and {B_s,1, ..., B_s,k , ..., B_s,d } can be fitted using equations (5) and (6).

The goal of this section is to go from the non-Markovian GLEs to a more hydrodynamic picture of the equations. For this, we introduce the microscopic number and momentum density fields of species s, containing a number of particles N_s (t) evolving in time according to the particle reactions:

$$\begin{equation}\begin{aligned}\hfill {\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)& =\sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right),\hfill \\ \hfill {\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)& =\sum _{i=1}^{{N}_{s}\left(t\right)}{\boldsymbol{p}}_{i,s}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right),\hfill \end{aligned}\end{equation} \tag{ 7 }$$

where N_s is the number of particles in species s.

To obtain the time-evolution equations governing the dynamics of the microscopic densities just introduced, equation (7) must be differentiated with respect to time, which can be done by applying the discrete Leibniz integral rule:

$$\begin{equation}{\partial }_{t}\sum _{i=1}^{{N}_{s}\left(t\right)}\boldsymbol{f}\left(i;t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)=\sum _{i=1}^{{N}_{s}\left(t\right)}{\partial }_{t}\left[\boldsymbol{f}\left(i;t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)\right]+\boldsymbol{f}\left({N}_{s};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{{N}_{s},s}\right)\frac{\mathrm{d}{N}_{s}\left(t\right)}{\mathrm{d}t}.\end{equation} \tag{ 8 }$$

It is convenient to express the derivative of the discrete upper summation limit in the distribution form, $\frac{\mathrm{d}{N}_{s}\left(t\right)}{\mathrm{d}t}=\kappa \delta \left(t\right)$, where κ indicates the reaction rate. The second term in the right-hand side of equation (8), which in our case will be stochastic in nature, can then be determined empirically [57].

Applying the time derivative defined in equation (8) to the equation (7), and using equation (1), we obtain the time evolution of the number and momentum density fields of each species s:

$$\begin{equation}{\partial }_{t}{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)=-{\nabla }_{\boldsymbol{r}}\cdot \left({m}_{s}^{-1}{\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)\right)+{\hat{w}}_{\rho ,s}\left(\boldsymbol{r};t\right),\end{equation} \tag{ 9 }$$

$$\begin{align}\hfill {\partial }_{t}{\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)& =-{\nabla }_{\boldsymbol{r}}\cdot \sum _{i=1}^{{N}_{s}\left(t\right)}\left({m}_{s}^{-1}{\boldsymbol{p}}_{i,s}\right)\otimes {\boldsymbol{p}}_{i,s}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)+\sum _{i=1}^{{N}_{s}\left(t\right)}{\boldsymbol{F}}_{i,s}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)\hfill \\ \hfill & \quad -{\int }_{{t}_{0}}^{t}\enspace \mathrm{d}t\prime {\boldsymbol{\theta }}_{s}\left(t-t\prime \right){\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)+{\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right)+{\hat{\boldsymbol{w}}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right),\hfill \end{align} \tag{ 10 }$$

where ${\hat{w}}_{\rho ,s}$ and ${\hat{\boldsymbol{w}}}_{\boldsymbol{j},s}$ are the contributions to density and momentum due to chemical reactions, and the local random force η _s is:

$$\begin{equation}{\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right)=\sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right){\boldsymbol{R}}_{i,s}\left(t\right).\end{equation} \tag{ 11 }$$

Notice that the terms ${\hat{w}}_{\rho ,s}$ and ${\hat{\boldsymbol{w}}}_{\boldsymbol{j},s}$ in general allow for reactions between the non-resolved variables and particle sources. However, if reactions can occur between resolved particles only, the conservation of mass and momentum of the total system must also be imposed:

$$\begin{equation}\begin{aligned}\hfill \sum _{\alpha =1}^{K}{\hat{w}}_{\rho ,\alpha }\left(\boldsymbol{r};t\right)=0,\\ \hfill \sum _{\alpha =1}^{K}{\hat{\boldsymbol{w}}}_{\boldsymbol{j},\alpha }\left(\boldsymbol{r};t\right)=\boldsymbol{0},\end{aligned}\end{equation} \tag{ 12 }$$

Finally, we can rewrite the noise term without changing its statistical properties (see appendix A) as [28, 29, 58–60]

$$\begin{equation}{\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right)=\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right),\end{equation} \tag{ 13 }$$

where R _s denotes the spatial-temporal noise satisfying the fluctuation-dissipation theorem

$$\begin{equation}\langle {\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right){\boldsymbol{R}}_{s}\left({\boldsymbol{r}}^{\prime };{t}^{\prime }\right)\rangle ={m}_{s}{k}_{\text{B}}T{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right),\end{equation} \tag{ 14 }$$

where k_B denotes the Boltzmann constant.

The first and the second terms in the right-hand side of equation (10) can be rewritten now in terms of ${\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)$ and ${\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)$. The second term with the force F _i,s , experienced by the particle can be further rewritten as

$$\begin{equation}\begin{aligned}\hfill \sum _{i=1}^{{N}_{s}\left(t\right)}{\nabla }_{{\boldsymbol{r}}_{i,s}}U\left({\boldsymbol{r}}_{i,s},{\boldsymbol{r}}^{N}\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)=& {\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)\left[{\nabla }_{\boldsymbol{r}}{U}_{\text{ext}}\left(\boldsymbol{r};t\right)\right.\hfill \\ \hfill & +\left.\sum _{\alpha =1}^{K}\int \mathrm{d}{\boldsymbol{r}}^{\prime }{\nabla }_{\boldsymbol{r}}{V}_{\alpha }\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\left({\hat{\rho }}_{\alpha }\left({\boldsymbol{r}}^{\prime };t\right)-\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\right)\right].\hfill \end{aligned}\end{equation} \tag{ 15 }$$

The convective flux then becomes [29, 61]

$$\begin{align}\hfill \sum _{i=1}^{{N}_{s}\left(t\right)}\left({m}_{s}^{-1}{\boldsymbol{p}}_{i,s}\right)\otimes {\boldsymbol{p}}_{i,s}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\right)& =\frac{{\sum }_{i=1}^{{N}_{s}\left(t\right)}{\sum }_{j=1}^{{N}_{s}\left(t\right)}\left({m}_{s}^{-1}{\boldsymbol{p}}_{i,s}\right)\otimes {\boldsymbol{p}}_{i,s}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{j,s}\right)}{{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)}\hfill \\ \hfill & =\frac{\left({m}_{s}^{-1}\hat{{\boldsymbol{j}}_{s}}\left(\boldsymbol{r};t\right)\right)\otimes \hat{{\boldsymbol{j}}_{s}}\left(\boldsymbol{r};t\right)}{{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)}.\hfill \end{align} \tag{ 16 }$$

It is worth underlying that, as discussed extensively in references [29, 61], the approach followed to obtain the convective flux in equation (16) is mathematically rigorous only if a physical interpretation of the Dirac delta function is accepted. Specifically, if the phase space is partitioned in cells each containing a single particle, the Dirac delta can be approximated as a Kronecker delta normalized with respect to the volume of the cell. As highlighted in [29], such approximation is strictly valid only for particle systems with a repulsive hard core, which represent by far the most common kind of systems studied in classical particle physics.

Thus, the fluctuating momentum-moment equation reads

$$\begin{align}\hfill {\partial }_{t}{\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)=& -{\nabla }_{\boldsymbol{r}}\cdot \left(\frac{\left({m}_{s}^{-1}{\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)\right)\otimes {\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)}{{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)}\right)-{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)\left({\nabla }_{\boldsymbol{r}}{U}_{\text{ext}}\left(\boldsymbol{r};t\right)\right.\hfill \\ \hfill & -{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)\sum _{\alpha =1}^{K}\int \mathrm{d}{\boldsymbol{r}}^{\prime }{\nabla }_{\boldsymbol{r}}{V}_{\alpha }\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\left({\hat{\rho }}_{\alpha }\left({\boldsymbol{r}}^{\prime };t\right)-\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\right)\hfill \\ \hfill & -{\int }_{{t}_{0}}^{t}\enspace \mathrm{d}{t}^{\prime }{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right){\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)+\sqrt{{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)+{\hat{\boldsymbol{w}}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right).\hfill \end{align} \tag{ 17 }$$

The following definitions of the kinetic, internal and external energy functionals

$$\begin{align}\hfill {H}_{K}\left[\phi ,\boldsymbol{\pi }\right]& =\frac{1}{2}\int \mathrm{d}\boldsymbol{r}\frac{\boldsymbol{\pi }{\left(\boldsymbol{r}\right)}^{2}}{\phi \left(\boldsymbol{r}\right)},\hfill \\ \hfill {H}_{U}\left[{\phi }_{s}\right]& =\frac{1}{2}\sum _{\alpha =1}^{K}\iint \mathrm{d}\boldsymbol{r}\enspace \mathrm{d}{\boldsymbol{r}}^{\prime }{V}_{\alpha }\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\left({\phi }_{s}\left(\boldsymbol{r}\right){\phi }_{\alpha }\left({\boldsymbol{r}}^{\prime }\right)-\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right){\phi }_{s}\left(\boldsymbol{r}\right)\right),\hfill \\ \hfill {H}_{U}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left[\phi \right]& =\frac{1}{2}\int \mathrm{d}\boldsymbol{r}\enspace {U}_{\text{ext}}\left(\boldsymbol{r}\right)\phi \left(\boldsymbol{r}\right),\hfill \end{align} \tag{ 18 }$$

allow us to rewrite equation (17) in terms of the functional derivative of equation (18) as

$$\begin{align}\hfill {\partial }_{t}{\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)=& -{\nabla }_{\boldsymbol{r}}\cdot \left(\frac{\left({m}_{s}^{-1}{\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)\right)\otimes {\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)}{{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)}\right)-{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)\left[{\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{U}\left[{\hat{\rho }}_{s}\right]}{\delta {\hat{\rho }}_{s}\left(\boldsymbol{r},t\right)}+{\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{U}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left[{\hat{\rho }}_{s}\right]}{\delta {\hat{\rho }}_{s}\left(\boldsymbol{r},t\right)}\right]\hfill \\ \hfill & -{\int }_{{t}_{0}}^{t}\enspace \mathrm{d}{t}^{\prime }{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right){\hat{\rho }}_{s}\left(\boldsymbol{r},t\right)\left[{\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{K}\left[{\hat{\rho }}_{s},{\hat{\boldsymbol{j}}}_{s}\right]}{\delta {\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r},t\right)}\right]+\sqrt{{\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)+{\hat{\boldsymbol{w}}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right).\hfill \end{align} \tag{ 19 }$$

Notice that equation (19) still depends on all particle positions and momenta through the terms ${\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)$ and ${\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)$. Thus, a proper ensemble average is required to obtain a mesoscopic description of the system [29]. With this aim, Durán-Olivencia et al extended the concept of equilibrium canonical distribution, well known from classical statistical mechanical arguments, to local-equilibrium conditions. Local-equilibrium conditions are out-of-equilibrium configurations of the system under which the particles, despite having small deviations from their equilibrium velocities and energies, behave locally as in equilibrium. This approach allows us to take the ensemble average of the system out-of-equilibrium, under a local-equilibrium approximation. To this end, we introduce a probability density function f_le( r ^N , p ^N ; t) for the system in local-equilibrium condition:

$$\begin{equation}{f}_{\text{le}}\left({\boldsymbol{r}}^{N},{\boldsymbol{p}}^{N};t\right)={\mathcal{Z}}^{-1}{\mathrm{e}}^{-\beta \mathcal{H}\left({\boldsymbol{r}}^{N},{\boldsymbol{p}}^{N};t\right)}\end{equation} \tag{ 20 }$$

where $\mathcal{Z}$ is the partition function, $\beta ={\left({k}_{\text{B}}T\right)}^{-1}$ and $\mathcal{H}\left({\boldsymbol{r}}^{N},{\boldsymbol{p}}^{N};t\right)$ is the Hamiltonian given by:

$$\begin{align}\hfill \mathcal{H}=& \frac{1}{2}\int \mathrm{d}\boldsymbol{r}\sum _{i=1}^{{N}_{s}\left(t\right)}\left[{m}^{-1}{\boldsymbol{p}}_{i}\cdot {\boldsymbol{p}}_{i}+\sum _{\alpha =1,\alpha \ne s}^{K}\sum _{j=1}^{{N}_{\alpha }\left(t\right)}{V}_{\alpha }\left({\boldsymbol{r}}_{i,s}-{\boldsymbol{r}}_{j}\right)\right.\hfill \\ \hfill & +\left.\sum _{j=1,j\ne i}^{{N}_{s}\left(t\right)}{V}_{s}\left({\boldsymbol{r}}_{i,s}-{\boldsymbol{r}}_{j}\right)+{U}_{\text{ext}}\left({\boldsymbol{r}}_{i,s}\right)\right]\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\right).\hfill \end{align} \tag{ 21 }$$

Hence, the ensemble average of density, momenta and reaction sources become fields:

$$\begin{equation}\left(\begin{matrix}\hfill {\rho }_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {w}_{\rho ,s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\boldsymbol{w}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right)\hfill \end{matrix}\right)=\int {f}_{\text{le}}\left({\boldsymbol{r}}^{N},{\boldsymbol{p}}^{N};t\right)\left(\begin{matrix}\hfill {\hat{\rho }}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\hat{w}}_{\rho ,s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\hat{\boldsymbol{j}}}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\hat{\boldsymbol{w}}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right)\hfill \end{matrix}\right)\mathrm{d}{\boldsymbol{r}}^{N}\mathrm{d}{\boldsymbol{p}}^{N}.\end{equation} \tag{ 22 }$$

Ensemble averaging equation (19) yields the mesoscopic FDDFT for species s:

$$\begin{equation}{\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)=-{\nabla }_{\boldsymbol{r}}\cdot \left({m}_{s}^{-1}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\right)+{w}_{\rho ,s}\left(\boldsymbol{r};t\right)\end{equation} \tag{ 23 }$$

$$\begin{align}\hfill {\partial }_{t}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)& =-{\nabla }_{\boldsymbol{r}}\cdot \left(\frac{\left({m}_{s}^{-1}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\right)\otimes {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)}{{\rho }_{s}\left(\boldsymbol{r};t\right)}\right)-{\rho }_{s}\left(\boldsymbol{r};t\right)\left[{\nabla }_{\boldsymbol{r}}\frac{\delta {\mathcal{F}}_{s}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}+{\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{U}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}\right]\hfill \\ \hfill & \quad -{\int }_{{t}_{0}}^{t}\enspace \mathrm{d}t\prime {\boldsymbol{\theta }}_{s}\left(t-t\prime \right){\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\prime \right)+\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)+{\boldsymbol{w}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right),\hfill \end{align} \tag{ 24 }$$

where we introduced the free-energy functional, defined as sum of ideal and excess-over-ideal contributions ${\mathcal{F}}_{s}\left[\rho \right]={\mathcal{F}}_{s}^{\mathrm{i}\mathrm{d}}\left[{\rho }_{s}\right]+{\mathcal{F}}_{s}{\left[{\rho }_{s}\right]}^{\mathrm{e}\mathrm{x}}$ [29]:

$$\begin{equation}{\mathcal{F}}_{s}^{\mathrm{i}\mathrm{d}}\left[{\rho }_{s}\right]={k}_{\text{B}}T\int \mathrm{d}\boldsymbol{r}\enspace {\rho }_{s}\left(\boldsymbol{r};t\right)\left(\mathrm{log}\left({{\Lambda}}^{3}{\rho }_{s}\left(\boldsymbol{r};t\right)\right)-1\right)\end{equation} \tag{ 25 }$$

with Λ the thermal de Broglie wavelength. The excess part of the free energy is defined through the relation:

$$\begin{equation}{\rho }_{s}\left(\boldsymbol{r};t\right){\nabla }_{\boldsymbol{r}}\frac{\delta {\mathcal{F}}_{s}^{\mathrm{e}\mathrm{x}}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}=\left\langle {\hat{\rho }}_{s}\left(\boldsymbol{r};t\right){\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{U}\left[{\hat{\rho }}_{s}\right]}{\delta {\hat{\rho }}_{s}\left(\boldsymbol{r},t\right)}\right\rangle =\sum _{\alpha =1}^{K}\int \mathrm{d}{\boldsymbol{r}}^{\prime }{\nabla }_{\boldsymbol{r}}{V}_{\alpha }\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right){\rho }_{s,\alpha }^{\left(2\right)}\left(\boldsymbol{r},{\boldsymbol{r}}^{\prime };t\right)\end{equation} \tag{ 26 }$$

where ${\rho }_{s,\alpha }^{\left(2\right)}\left(\boldsymbol{r},{\boldsymbol{r}}^{\prime };t\right)=\left\langle {\rho }_{s}\left(\boldsymbol{r};t\right){\rho }_{\alpha }\left({\boldsymbol{r}}^{\prime };t\right)\right\rangle$ is the inter-species pair-correlation function. The functional ${\mathcal{F}}_{s}\left[\rho \right]$ satisfies the relation:

$$\begin{equation}{\rho }_{s}\left(\boldsymbol{r};t\right){\nabla }_{\boldsymbol{r}}\frac{\delta {\mathcal{F}}_{s}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}={k}_{\mathrm{B}}T{\nabla }_{\boldsymbol{r}}{\rho }_{s}\left(\boldsymbol{r},t\right)+\sum _{\alpha =1}^{K}\int \mathrm{d}{\boldsymbol{r}}^{\prime }{\nabla }_{\boldsymbol{r}}{V}_{\alpha }\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right){\rho }_{s,\alpha }^{\left(2\right)}\left(\boldsymbol{r},{\boldsymbol{r}}^{\prime };t\right)\end{equation} \tag{ 27 }$$

which represents the first equation of the Yvon–Born–Green hierarchy if the averages are taken at equilibrium. For this reason, the functional ${\mathcal{F}}_{s}\left[\rho \right]$ can be interpreted as the classical Helmholtz free-energy functional (a more extensive discussion is available in reference [29] and references therein). Notice that in local equilibrium, the term containing the Helmholtz free energy is associated with pressure. In fact, from the Gibbs–Duhem relation (dp = ∑_s ρ_s dμ_s ), and the Euler–Lagrange equation of equilibrium density functional theory (${\mu }_{s}=\delta {\mathcal{F}}_{s}\left[{\rho }_{s}\right]/\delta {\rho }_{s}$) we get:

$$\begin{equation}{\nabla }_{\boldsymbol{r}}p=\sum _{s}{\rho }_{s}\left(\boldsymbol{r};t\right){\nabla }_{\boldsymbol{r}}\frac{\delta {\mathcal{F}}_{s}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)},\end{equation} \tag{ 28 }$$

which is consistent with the well-known expressions for one-component systems [29, 59].

The terms w_ρ,s and w _{j ,s} represent the average density and momentum reaction rates. The density reaction rate w_ρ,s depends on the average number of particle changing species, and thus it can be modeled or directly measured from experiments [see for instance equation (76)]. The momentum reaction rate w _{j ,s} accounts for the average momentum transfer rate between particles during the reaction events. If reactions occur only between resolved particles, i.e. for systems with non-reactive thermal baths, the conservation of mass and momentum of the total system is satisfied, i.e. ${\sum }_{\alpha =1}^{K}{w}_{\rho ,\alpha }\left(\boldsymbol{r};t\right)=0$ and ${\sum }_{\alpha =1}^{K}{\boldsymbol{w}}_{\boldsymbol{j},\alpha }\left(\boldsymbol{r};t\right)=\boldsymbol{0}$. Moreover, for isotropic reactions, i.e. when the reaction depends only on the random collisions between the reacting particles, the momentum transfer rate w _{j ,s} has zero components. However, when reactions are favored along some specific directions, for instance in systems made of aligned polar particles, then the term w _{j ,s} will actively contribute to the evolution of the species' momentum.

In many cases the memory kernel decays to zero on a timescale much smaller than the characteristic timescale of the density dynamics. Under such circumstances, the density field relaxes much slower than the momentum field, and it is possible to assume that the material derivative of the momentum vanishes, D_t j _s ∼ 0. This approximation allows us to reduce the equations of FH to a single evolution equation for the density. Our system, however, contains a non-Markovian memory term which complicates taking the overdamped limit. We proceed by defining the following deconvolution operator:

$$\begin{equation}\mathcal{I}\left[{\int }_{{t}_{0}}^{t}\enspace \mathrm{d}{t}^{\prime }\phi \left(t-{t}^{\prime }\right)\boldsymbol{\pi }\left(\boldsymbol{r};t\right),\phi \right]=\boldsymbol{\pi }\left(\boldsymbol{r};t\right),\end{equation} \tag{ 29 }$$

which is applied to both sides of the momentum equation (24). With the above assumption D_t j _s ∼ 0, we obtain the non-Markovian time evolution of the density field of the species s in a multi-component system:

$$\begin{align}\hfill {\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)=& {\nabla }_{\boldsymbol{r}}\cdot \left\{{m}_{s}^{-1}\mathcal{I}\left[{\rho }_{s}\left(\boldsymbol{r};t\right)\left({\nabla }_{\boldsymbol{r}}\frac{\delta {\mathcal{F}}_{s}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}+{\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{U}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}\right)\right.\right.\hfill \\ \hfill & \left.+\left.\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)+{\boldsymbol{w}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right)\right]\right\}+{w}_{\rho ,s}\left(\boldsymbol{r};t\right).\hfill \end{align} \tag{ 30 }$$

In general, the operator $\mathcal{I}$ may be arbitrarily defined, with the only restriction that it needs to satisfy equation (29). In this work, we take advantage of the properties of the Laplace transform and define the operator $\mathcal{I}$ as:

$$\begin{align}\hfill \mathcal{I}\left[{\int }_{{t}_{0}}^{t}\enspace \mathrm{d}{t}^{\prime }\phi \left(t-{t}^{\prime }\right)\boldsymbol{\pi }\left(\boldsymbol{r};t\right),\phi \right]& ={\mathcal{L}}^{-1}\left\{\frac{\mathcal{L}\left[{\int }_{{t}_{0}}^{t}\enspace \mathrm{d}{t}^{\prime }\phi \left(t-{t}^{\prime }\right)\boldsymbol{\pi }\left(\boldsymbol{r};t\right)\right]}{\mathcal{L}\left[\phi \left(t\right)\right]}\right\}\hfill \\ \hfill & ={\mathcal{L}}^{-1}\left\{\frac{\mathcal{L}\left[\phi \left(t\right)\right]\mathcal{L}\left[\boldsymbol{\pi }\left(\boldsymbol{r};t\right)\right]}{\mathcal{L}\left[\phi \left(t\right)\right]}\right\}=\boldsymbol{\pi }\left(\boldsymbol{r};t\right)\hfill \end{align} \tag{ 31 }$$

where $\mathcal{L}\left[\boldsymbol{\pi }\left(\boldsymbol{r};t\right)\right]{=\int }_{0}^{\infty }\boldsymbol{\pi }\left(\boldsymbol{r};t\right){\mathrm{e}}^{-t/\lambda }\enspace \mathrm{d}t$ indicates the Laplace transform of a time-dependent field.

Applying the operator $\mathcal{I}$ in equation (29) for a multicomponent system of ideal gases, we obtain the following non-Markovian time evolution of the density field of the species s:

$$\begin{equation}\begin{aligned}\hfill {\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)={\nabla }_{\boldsymbol{r}}\cdot \left\{{m}_{s}^{-1}\mathcal{I}\left[{\nabla }_{\boldsymbol{r}}{k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)+\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)+{\boldsymbol{w}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right)\right]\right\}+{w}_{\rho ,s}\left(\boldsymbol{r};t\right).\end{aligned}\end{equation} \tag{ 32 }$$

Neglecting the reaction contributions to the momentum, we obtain:

$$\begin{equation}{\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)={\nabla }_{\boldsymbol{r}}\cdot \left\{{m}_{s}^{-1}\mathcal{I}\left[{k}_{\text{B}}T{\nabla }_{\boldsymbol{r}}{\rho }_{s}\left(\boldsymbol{r};t\right)+\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)\right]\right\}+{w}_{\rho ,s}\left(\boldsymbol{r};t\right)\end{equation} \tag{ 33 }$$

where we employed the space invariance of the operator $\mathcal{I}$. By utilising equation (31), we obtain

$$\begin{equation}{\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)={\nabla }_{\boldsymbol{r}}\cdot \left\{{m}_{s}^{-1}{\mathcal{L}}^{-1}\left[\frac{\mathcal{L}\left[{k}_{\text{B}}T{\nabla }_{\boldsymbol{r}}{\rho }_{s}\left(\boldsymbol{r};t\right)+\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)\right]}{\mathcal{L}\left[\theta \left(t\right)\right]}\right]\right\}+{w}_{\rho ,s}\left(\boldsymbol{r};t\right).\end{equation} \tag{ 34 }$$

Defining now the function $\tilde {D}\left(t\right)={\mathcal{L}}^{-1}\left[{\left(\mathcal{L}\left[\theta \left(t\right)\right]\right)}^{-1}\right]$, and using the properties of the Laplace transform we get

$$\begin{equation}{\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)=\tilde {D}\left(t\right){\ast}{\nabla }_{\boldsymbol{r}}\cdot \left\{\left({k}_{\text{B}}T{m}_{s}^{-1}{\nabla }_{\boldsymbol{r}}{\rho }_{s}\left(\boldsymbol{r};t\right)+\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)\right)\right\}+{w}_{\rho ,s}\left(\boldsymbol{r};t\right),\end{equation} \tag{ 35 }$$

where ∗ denotes convolution in time.

For a memory kernel in the form ${\theta }_{s}\left(t\right)={\sum }_{k}^{d}{A}_{s,k}{\mathrm{e}}^{{B}_{s,k}t}$, the function $\tilde {D}\left(t\right)$ is given by

$$\begin{equation}\tilde {D}\left(t\right)={\mathcal{L}}^{-1}\left[{\left(\sum _{k}^{d}\frac{{A}_{s,k}\lambda }{1-{B}_{s,k}\lambda }\right)}^{-1}\right].\end{equation} \tag{ 36 }$$

The highest values of λ in the Laplace transform correspond to the lowest frequencies, which are usually dominant. This allows us to simplify the above even further:

$$\begin{equation}\tilde {D}\left(t\right)\sim {\mathcal{L}}^{-1}\left[\underset{\lambda \to \infty }{\mathrm{lim}}{\left(\sum _{k}^{d}\frac{{A}_{s,k}\lambda }{1-{B}_{s,k}\lambda }\right)}^{-1}\right]={\mathcal{L}}^{-1}\left[\sum _{k}^{d}\frac{{B}_{s,k}}{{A}_{s,k}}\right]=\sum _{k}^{d}\frac{{B}_{s,k}}{{A}_{s,k}}\delta \left(t\right).\end{equation} \tag{ 37 }$$

Finally, equation (35) reduces to

$$\begin{equation}{\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)={D}_{\text{eff}}\enspace {\nabla }_{\boldsymbol{r}}\cdot \left\{{k}_{\text{B}}T{m}_{s}^{-1}{\nabla }_{\boldsymbol{r}}{\rho }_{s}\left(\boldsymbol{r};t\right)+\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right)\right\}+{w}_{\rho ,s}\left(\boldsymbol{r};t\right),\end{equation} \tag{ 38 }$$

where we have introduced the effective diffusion coefficient, ${D}_{\text{eff}}={\sum }_{k}^{d}\frac{{B}_{s,k}}{{A}_{s,k}}$.

An alternative derivation of equation (38) would be to first take the Markovian limit of an ideal-gas model followed by the overdamped limit assuming that the momentum reaction rate w _{j ,s} is negligible. This way equation (30) would be reduced to the usual reaction–diffusion equation for the density [57]:

$$\begin{equation}\begin{aligned}\hfill {\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)={\nabla }_{\boldsymbol{r}}\cdot \left[{k}_{\text{B}}T{\left({m}_{s}{\boldsymbol{\theta }}_{0,s}\right)}^{-1}\left({\nabla }_{\boldsymbol{r}}{\rho }_{s}\right)\right]+{\nabla }_{\boldsymbol{r}}\cdot \sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{0,s}\left(\boldsymbol{r};t\right)+{w}_{\rho ,s}\left(\boldsymbol{r};t\right).\end{aligned}\end{equation} \tag{ 39 }$$

In this section we have presented some of the main theoretical results of this work. Using a properly defined deconvolution operator and taking the overdamped limit of the previously derived FH, we derived a novel non-Markovian overdamped FDDFT equation for the density field. Moreover, for memory kernels in the form of exponential series, we simplified the convolutional structure of the above equation, defining an effective diffusion coefficient including low-frequencies contributions to memory effects.

Numerically solving the non-Markovian equations (23) and (24) is a highly non-trivial task. First, because it is not obvious how to sample colored (time-correlated) noise. Additionally, the direct evaluation of the convolution integral requires keeping track of the entire history of the momentum field, which drives up the computational cost. Within the so-called extended dynamics approach, additional variables are introduced to capture the memory and noise effects on the observables [56, 62, 63]. In this section we derive the extended dynamics for the non-Markovian FH, generalizing the work by Kawai [63]. For simplicity we assume that the global effects of the bath are isotropic [56] and the memory kernel can be expressed as the exponential series [62, 63] already introduced in section 2 and repeated here for convenience:

$$\begin{equation}\begin{aligned}\hfill {\boldsymbol{\theta }}_{s}\left(t\right)={\theta }_{s}\left(t\right)\mathbf{I}=\sum _{k}^{d}{A}_{s,k}{\mathrm{e}}^{{B}_{s,k}t}\mathbf{I},\quad {B}_{s,k}{\leqslant}0.\end{aligned}\end{equation} \tag{ 40 }$$

5.1. Convolution decomposition

Defining the extended field variable associated with the kth mode of the memory term of species s:

$$\begin{equation}{\boldsymbol{Z}}_{s,k}\left(\boldsymbol{r};t\right)={\int }_{{t}_{0}}^{t}\enspace \mathrm{d}{t}^{\prime }{A}_{s,k}{\mathrm{e}}^{{B}_{s,k}\left(t-{t}^{\prime }\right)}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};{t}^{\prime }\right),\end{equation} \tag{ 41 }$$

we can express the convolution term as ${\int }_{{t}_{0}}^{t}\enspace \mathrm{d}{t}^{\prime }{\theta }_{s}\left(t-{t}^{\prime }\right){\boldsymbol{j}}_{s}\left(\boldsymbol{r};{t}^{\prime }\right)={\sum }_{k}^{d}{\boldsymbol{Z}}_{s,k}\left(\boldsymbol{r};t\right)$. By differentiating Z_s,k with respect to t yields the time evolution of the field Z_s,k in the form of a stochastic partial differential equation (SPDE):

$$\begin{equation}{\partial }_{t}{\boldsymbol{Z}}_{s,k}\left(\boldsymbol{r};t\right)={B}_{s,k}\enspace {\boldsymbol{Z}}_{s,k}\left(\boldsymbol{r};t\right)+{A}_{s,k}\enspace {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right).\end{equation} \tag{ 42 }$$

5.2. Noise decomposition

The noise term is given by

$$\begin{equation}{\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right)=\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right),\end{equation} \tag{ 43 }$$

where R _s denotes the spatial-temporal noise with correlations ⟨ R _s ( r ; t) R _s ( r '; t')⟩ = k_B Tθ_s (t − t'δ( r − r ') I. Notice that because of the symmetry between t and t' in the fluctuation-dissipation theorem, θ_s (t) must be even: θ_s (t) = θ_s (−t). Additionally, its Fourier transform ${\tilde {\theta }}_{s}\left(\omega \right){=\int }_{-\infty }^{-\infty }{\theta }_{s}\left(t\right){\mathrm{e}}^{-\mathrm{i}\omega \tau }\mathrm{d}t$ must be both real and even when ω is real. It follows that in the ω-plane, the roots and singularities of ${\tilde {\theta }}_{s}\left(\omega \right)$ must be symmetric with respect to both the real and imaginary axes. This allows us to introduce the function $\tilde {\chi }\left(\omega \right)$ such that

$$\begin{equation}{\tilde {\theta }}_{s}\left(\omega \right)=2\tilde {\chi }\left(\omega \right){\tilde {\chi }}^{T}\left(-\omega \right),\end{equation} \tag{ 44 }$$

where

$$\begin{equation}\tilde {\chi }\left(\omega \right)=\sum _{k}-i{\left(\omega +i{B}_{s,k}\right)}^{-1}{b}_{s,k}\end{equation} \tag{ 45 }$$

with real coefficients b_s,k and B_s,k . The singular points of ${\tilde {\chi }}^{-1}\left(\omega \right)$ lie in the lower half of the complex ω-plane. Moreover, we define two functions in the Fourier space:

$$\begin{equation}\tilde {\zeta }\left(\omega \right)={\tilde {\chi }}^{-1}\left(\omega \right)=\sum _{k}\left(\omega +i{B}_{s,k}\right){\left(-i{b}_{s,k}\right)}^{-1},\end{equation} \tag{ 46 }$$

and

$$\begin{equation}{\tilde {k}}_{k}\left(\omega \right)=-i{\left(\omega +i{B}_{s,k}\right)}^{-1}{b}_{s,k}\tilde {\zeta }\left(\omega \right),\end{equation} \tag{ 47 }$$

and we denote their Fourier inverse transform with h(t) and k_k (t), respectively. By combining equations (45)–(47) it follows that:

$$\begin{equation}\sum _{k}{\tilde {k}}_{k}\left(\omega \right)=1\end{equation} \tag{ 48 }$$

or, equivalently,

$$\begin{equation}\sum _{k}{k}_{k}\left(t\right)=\delta \left(t\right).\end{equation} \tag{ 49 }$$

Moreover, equation (47) can be rewritten as $\left(i\omega -{B}_{s,k}\right){\tilde {k}}_{k}\left(\omega \right)={b}_{s,k}\tilde {\zeta }\left(\omega \right)$, which in the time domain gives:

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{k}_{k}\left(t\right)-{B}_{s,k}{k}_{k}\left(t\right)={b}_{s,k}\zeta \left(t\right).\end{equation} \tag{ 50 }$$

We introduce two more vector fields, corresponding to a delta correlated stochastic field (as will be proved later) and an auxiliary stochastic field, respectively:

$$\begin{equation}{\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};t\right)={\int }_{0}^{+\infty }\zeta \left(t-{t}^{\prime }\right){\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};{t}^{\prime }\right)\mathrm{d}{t}^{\prime }\end{equation} \tag{ 51 }$$

and

$$\begin{equation}{\boldsymbol{\eta }}_{s,k}\left(\boldsymbol{r};t\right)={\int }_{0}^{+\infty }{k}_{k}\left(t-{t}^{\prime }\right){\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};{t}^{\prime }\right)\mathrm{d}{t}^{\prime }.\end{equation} \tag{ 52 }$$

From equations (52) and (49) it follows that:

$$\begin{equation}\sum _{k}^{d}{\boldsymbol{\eta }}_{s,k}\left(\boldsymbol{r};t\right)={\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right),\end{equation} \tag{ 53 }$$

while, combining equations (52) and (50), we obtain

$$\begin{equation}{\partial }_{t}{\boldsymbol{\eta }}_{s,k}\left(\boldsymbol{r};t\right)={B}_{s,k}{\boldsymbol{\eta }}_{s,k}\left(\boldsymbol{r};t\right)+{b}_{s,k}{\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};t\right).\end{equation} \tag{ 54 }$$

Equations (53) and (54) allow us to express the original correlated noise as a function of the white noise ξ (t). A discussion of the properties of the stochastic process ξ ( r ; t) can be found in appendix B.

5.3. Non-Markovian FH equations

Defining an extended field S _s,k ( r ; t) = − Z _s,k ( r ; t) + η _s,k ( r ; t), the FH equations are expressed in the following form:

$$\begin{equation}{\partial }_{t}{\rho }_{s}\left(\boldsymbol{r};t\right)=-{\nabla }_{\boldsymbol{r}}\cdot \left({m}_{s}^{-1}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\right)+{w}_{\rho ,s}\left(\boldsymbol{r};t\right)\end{equation} \tag{ 55 }$$

$$\begin{align}\hfill {\partial }_{t}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)=& -{\nabla }_{\boldsymbol{r}}\cdot \left(\frac{\left({m}_{s}^{-1}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\right)\otimes {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)}{{\rho }_{s}\left(\boldsymbol{r};t\right)}\right)-{\rho }_{s}\left(\boldsymbol{r};t\right)\left[{\nabla }_{\boldsymbol{r}}\frac{\delta \mathcal{F}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}+{\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{U}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}\right]\hfill \\ \hfill & \quad +\sum _{k}^{d}{\boldsymbol{S}}_{s,k}\left(\boldsymbol{r};t\right)+{\boldsymbol{w}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right)\hfill \end{align} \tag{ 56 }$$

$$\begin{equation}{\partial }_{t}{\boldsymbol{S}}_{s,k}\left(\boldsymbol{r};t\right)={B}_{s,k}\enspace {\boldsymbol{S}}_{s,k}\left(\boldsymbol{r};t\right)-{A}_{s,k}\enspace {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)+{b}_{k}{\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};t\right)\end{equation} \tag{ 57 }$$

where ξ ( r ₁; t₁) is a stochastic process with zero mean and correlations $\left\langle {\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};{t}_{1}\right){\boldsymbol{\xi }}_{s}\left({\boldsymbol{r}}^{\prime };{t}_{2}\right)\right\rangle =2{k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\delta \left({t}_{1}-{t}_{2}\right)\mathrm{m}\mathrm{a}\mathrm{t}\mathrm{h}\mathrm{b}\mathrm{f}\mathrm{I}$.

The system in equation (57) can be also rewritten in vector form:

$$\begin{equation}{\partial }_{t}{\mathcal{U}}_{s}\left(\boldsymbol{r};t\right)=-{\nabla }_{\boldsymbol{r}}\cdot {\mathcal{M}}_{s}\left(\boldsymbol{r};t\right)+{\mathcal{N}}_{s}\left(\boldsymbol{r};t\right)+{\mathcal{R}}_{s}\left(\boldsymbol{r};t\right)+{\mathcal{Q}}_{s}.\end{equation} \tag{ 58 }$$

In the above equation ${\mathcal{U}}_{s}$ is the unknown, ${\mathcal{M}}_{s}$ is the conservative flux, and ${\mathcal{Q}}_{s}$ is the stochastic term for species s:

$$\begin{equation}{\mathcal{U}}_{s}=\left(\begin{matrix}\hfill {\rho }_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\boldsymbol{S}}_{s,1}\left(\boldsymbol{r};t\right)\hfill \\ \hfill \dots \hfill \\ \hfill {\boldsymbol{S}}_{s,k}\left(\boldsymbol{r};t\right)\hfill \\ \hfill \dots \hfill \end{matrix}\right),\quad {\mathcal{M}}_{s}=\left(\begin{matrix}\hfill \left({m}_{s}^{-1}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\right)\hfill \\ \hfill \frac{\left({m}_{s}^{-1}{\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\right)\otimes {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)}{{\rho }_{s}\left(\boldsymbol{r};t\right)}\hfill \\ \hfill 0\hfill \\ \hfill \dots \hfill \\ \hfill 0\hfill \\ \hfill \dots \hfill \end{matrix}\right),\quad {\mathcal{Q}}_{s}=\left(\begin{matrix}\hfill 0\hfill \\ \hfill \mathbf{0}\hfill \\ \hfill {b}_{1}\enspace {\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill \dots \hfill \\ \hfill {b}_{k}\enspace {\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill \dots \hfill \end{matrix}\right).\end{equation} \tag{ 59 }$$

The source terms are

$$\begin{equation}{\mathcal{N}}_{s}=\left(\begin{matrix}\hfill 0\hfill \\ \hfill -{\rho }_{s}\left(\boldsymbol{r};t\right)\left[{\nabla }_{\boldsymbol{r}}\frac{\delta \mathcal{F}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}+{\nabla }_{\boldsymbol{r}}\frac{\delta {H}_{U}^{\mathrm{e}\mathrm{x}\mathrm{t}}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}\right]+\sum _{k}^{d}{\boldsymbol{S}}_{s,k}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {B}_{s,1}\enspace {\boldsymbol{S}}_{s,1}\left(\boldsymbol{r};t\right)-{A}_{s,1}\enspace {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill \dots \hfill \\ \hfill {B}_{s,k}\enspace {\boldsymbol{S}}_{s,k}\left(\boldsymbol{r};t\right)-{A}_{s,k}\enspace {\boldsymbol{j}}_{s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill \dots \hfill \end{matrix}\right),\quad {\mathcal{R}}_{s}=\left(\begin{matrix}\hfill {w}_{\rho ,s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\boldsymbol{w}}_{\boldsymbol{j},s}\left(\boldsymbol{r};t\right)\hfill \\ \hfill 0\hfill \\ \hfill \dots \hfill \\ \hfill 0\hfill \\ \hfill \dots \hfill \end{matrix}\right).\end{equation} \tag{ 60 }$$

The FH system of equations and the corresponding deterministic and overdamped counterparts have been solved numerically using different techniques, including pseudo-spectral methods [64], and finite volumes [57, 65–67]. Here we employ the finite-volume method to discretize equation (58) because of its enhanced capability in dealing with non-regular fields and complex geometries. For simplicity, we present the discretization for a one dimensional system; the generalization to higher dimensions is then straightforward.

A domain of length L is partitioned into n identical cells of width Δx = L/n. The value of the discretized variable vector for the species s at the cell 1 ⩽ j ⩽ n is then defined as the following average over the cell volume V_j

$$\begin{equation}{\mathcal{U}}_{s,j}\left(t\right)=\frac{1}{{\Delta}x}{\int }_{{V}_{j}}{\mathcal{U}}_{s}\left(x,t\right)\enspace \mathrm{d}x.\end{equation} \tag{ 61 }$$

Similarly, time integration is performed by discretizing the total simulation time T in steps of magnitude Δt, such that we can assume ${\mathcal{U}}_{s,j}^{m}\equiv {\mathcal{U}}_{s,j}\left(m{\Delta}t\right)$, with 0 ⩽ mΔt ⩽ T. The multiplicative white noise ξ _s ( r ; t) is discretized, using the following spatio–temporal average [57, 65]

$$\begin{equation}{\boldsymbol{\xi }}_{s,j}^{m}\left(t\right)\equiv \frac{1}{{\Delta}x{\Delta}t}{\int }_{m{\Delta}t}^{\left(m+1\right){\Delta}t}{\int }_{{V}_{j}}\boldsymbol{\xi }\left(x,t\right)\enspace \mathrm{d}x\sim \sqrt{\frac{{\rho }_{s,j}\left(t\right)}{{\Delta}x{\Delta}t}}\enspace \boldsymbol{G}\end{equation} \tag{ 62 }$$

where G represents a vector with components sampled independently from a standard Gaussian distribution $\mathcal{N}\left(0,1\right)$. In order to implement time stepping, we define the forward and backward conservative fluxes as:

$$\begin{equation}\nabla \cdot {\mathcal{M}}_{s}\left(x;t\right)=\frac{{\mathcal{M}}_{s}\left(x;t\right)}{\mathrm{d}x}\sim \frac{{\mathcal{M}}_{s,j{\pm}1}\left(t\right)\mp {\mathcal{M}}_{s,j}\left(t\right)}{{\Delta}x}\end{equation} \tag{ 63 }$$

where ${\mathcal{M}}_{s,j}\left(t\right)$ denotes the average value of the field ${\mathcal{M}}_{s}\left(x;t\right)$ at the cell j, i.e.

$$\begin{equation}{\mathcal{M}}_{s}\left({\mathcal{U}}_{s,j}\left(x;t\right)\right)={\mathcal{M}}_{s,j}\left(t\right).\end{equation} \tag{ 64 }$$

Finally, both source terms are discretized simply:

$$\begin{equation}\begin{aligned}\hfill {\mathcal{N}}_{s}\left({\mathcal{U}}_{s,j}\left(x;t\right)\right)\equiv {\mathcal{N}}_{s,j}\left(t\right)\\ \hfill {\mathcal{R}}_{s}\left({\mathcal{U}}_{s,j}\left(x;t\right)\right)\equiv {\mathcal{R}}_{s,j}\left(t\right).\end{aligned}\end{equation} \tag{ 65 }$$

The source term ${\mathcal{N}}_{s}\left(x;t\right).$ contains the gradient of the variation of the free-energy functional. Within mean-field treatment, the latter has the following form:

$$\begin{equation}\frac{\delta \mathcal{F}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}\sim \frac{\delta {\mathcal{F}}^{id}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}+\sum _{\alpha =1}^{K}\int \mathrm{d}{\boldsymbol{r}}^{\prime }{V}_{\alpha }\left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right){\rho }_{\alpha }\left({\boldsymbol{r}}^{\prime };t\right).\end{equation} \tag{ 66 }$$

The above equation can be discretized as follows [66]:

$$\begin{equation}{\left(\frac{\delta \mathcal{F}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}\right)}_{i}\sim {\left.\frac{\delta {\mathcal{F}}^{\mathrm{i}\mathrm{d}}\left[{\rho }_{s}\right]}{\delta {\rho }_{s}\left(\boldsymbol{r},t\right)}\right\vert }_{{\rho }_{s}={\rho }_{s,i}}+{\Delta}x\sum _{\alpha =1}^{K}\sum _{i}{V}_{\alpha ,j-i}\left(t\right){\rho }_{\alpha ,i}\left(t\right).\end{equation} \tag{ 67 }$$

To improve the accuracy of time integration for nonlinear systems, each time step is usually further subdivided. In the present work, we adopt a two-stage MacCormack predictor–corrector scheme [68]:

$$\begin{equation}\begin{aligned}\hfill {\tilde {\mathcal{U}}}_{s,j}^{m}& ={\mathcal{U}}_{s,j}^{m}-\frac{{\Delta}t}{{\Delta}x}\left({\mathcal{M}}_{s,j}^{m}-{\mathcal{M}}_{s,j-1}^{m}\right)+{\Delta}t\left({\mathcal{N}}_{s,j}^{m}+{\mathcal{Q}}_{s,j}^{m}+{\mathcal{R}}_{s,j}^{m}\right)\hfill \\ \hfill {\mathcal{U}}_{s,j}^{m+1}& =\frac{1}{2}\left[{\mathcal{U}}_{s,j}^{m}+{\tilde {\mathcal{U}}}_{s,j}^{m}-\frac{{\Delta}t}{{\Delta}x}\left({\tilde {\mathcal{M}}}_{s,j+1}^{m}-{\tilde {\mathcal{M}}}_{s,j}^{m}\right)+{\Delta}t\left({\mathcal{N}}_{s,j}^{m}+{\mathcal{Q}}_{s,j}^{m}+{\mathcal{R}}_{s,j}^{m}\right)\right].\hfill \end{aligned}\end{equation} \tag{ 68 }$$

7.1. Equilibrium mono-component system

In order to test the numerical stochastic integrator, we consider a non-Markovian FDDFT approach to model a system of N non-reacting ideal gas particles with single exponential memory kernel. The particle momentum time correlation function of the corresponding GLE can be obtained in a closed analytical form [62, 69]:

$$\begin{equation}{\hat{C}}_{t,i}\left(t\right)=\frac{\langle {\boldsymbol{p}}_{i}\left(t\right){\boldsymbol{p}}_{i}\left(0\right)\rangle }{\langle {\boldsymbol{p}}_{i}\left(0\right){\boldsymbol{p}}_{i}\left(0\right)\rangle }=\begin{cases}{\mathrm{e}}^{\frac{tB}{2}}\left(\mathrm{cos}\left({\Omega}t\right)-\frac{B}{2{\Omega}}\mathrm{sin}\left({\Omega}t\right)\right)\quad \hfill & {\Omega}\ne 0,\hfill \\ {\mathrm{e}}^{\frac{tB}{2}}\left(1-\frac{Bt}{2}\right)\quad \hfill & {\Omega}=0,\hfill \end{cases}\end{equation} \tag{ 69 }$$

where we introduced a complex-valued parameter ${\Omega}=\sqrt{A-{B}^{2}/4}$.

To obtain the time-correlation function for the discretized non-Markovian FDDFT model, we proceed as follows. First, we define the time-correlation function of the discretized momentum in a cell as:

$$\begin{equation}{C}_{t,j}\left(t\right)=\left\langle \frac{1}{{\Delta}{x}^{2}}{\int }_{{V}_{j}}\mathrm{d}\boldsymbol{r}\sum _{i=1}^{{N}_{s}\left(t\right)}{\boldsymbol{p}}_{i}\left(t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\left(t\right)\right){\int }_{{V}_{j}}\mathrm{d}{\boldsymbol{r}}^{\prime }\sum _{j=1}^{{N}_{s}\left(t\right)}{\boldsymbol{p}}_{j}\left(0\right)\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{j}\left(0\right)\right)\right\rangle .\end{equation} \tag{ 70 }$$

Next, we employ the independence between p _i (t) and r _i (t), and between p _i (t) and p _j (t), to write

$$\begin{align}\hfill {C}_{t,j}\left(t\right)& =\frac{1}{{\Delta}{x}^{2}}{\int }_{{V}_{j}}\mathrm{d}\boldsymbol{r}{\int }_{{V}_{j}}\mathrm{d}{\boldsymbol{r}}^{\prime }\sum _{i=1}^{{N}_{s}\left(t\right)}\sum _{j=1}^{{N}_{s}\left(t\right)}\left\langle {\boldsymbol{p}}_{i}\left(t\right){\boldsymbol{p}}_{j}\left(0\right)\right\rangle \left\langle \delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\left(t\right)\right)\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{j}\left(0\right)\right)\right\rangle \hfill \\ \hfill & =\frac{1}{{\Delta}{x}^{2}}{\int }_{{V}_{j}}\mathrm{d}\boldsymbol{r}{\int }_{{V}_{j}}\mathrm{d}{\boldsymbol{r}}^{\prime }\sum _{i=1}^{{N}_{s}\left(t\right)}\left\langle {\boldsymbol{p}}_{i}\left(t\right){\boldsymbol{p}}_{i}\left(0\right)\right\rangle \left\langle \delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\left(t\right)\right)\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{i}\left(0\right)\right)\right\rangle .\hfill \end{align} \tag{ 71 }$$

Moreover, since all particles have the same momentum correlation, we introduce Ĉ_t (t) = Ĉ_t,i (t). This gives:

$$\begin{align}\hfill {C}_{t,j}\left(t\right)& =\frac{{\hat{C}}_{t}\left(t\right)}{{\Delta}{x}^{2}}{\int }_{{V}_{j}}\mathrm{d}\boldsymbol{r}{\int }_{{V}_{j}}\mathrm{d}{\boldsymbol{r}}^{\prime }\sum _{i=1}^{{N}_{s}\left(t\right)}\left\langle \delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\left(t\right)\right)\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{i}\left(0\right)\right)\right\rangle \hfill \\ \hfill & =\frac{{\hat{C}}_{t}\left(t\right)}{{\Delta}{x}^{2}}{\int }_{{V}_{j}}\mathrm{d}\boldsymbol{r}\left\langle \sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\left(t\right)\right){\int }_{{V}_{j}}\mathrm{d}{\boldsymbol{r}}^{\prime }\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{i}\left(0\right)\right)\right\rangle .\hfill \end{align} \tag{ 72 }$$

Finally, assuming that both r _i (0) and r _i (t) belong to the same cell V_j , i.e. that the cell size is much larger than the particle mean path, we obtain the time-correlation of the momentum field:

$$\begin{equation}{C}_{t,j}\left(t\right)=\frac{{\hat{C}}_{t}\left(t\right)}{{\Delta}{x}^{2}}{\int }_{{V}_{j}}\mathrm{d}\boldsymbol{r}\left\langle \sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i}\left(t\right)\right)\right\rangle =\frac{{\hat{C}}_{t}\left(t\right)}{{\Delta}x}\left\langle \frac{1}{{\Delta}x}{\int }_{{V}_{j}}\mathrm{d}\boldsymbol{r}\enspace \hat{\rho }\left(\boldsymbol{r},t\right)\right\rangle ={\hat{C}}_{t}\left(t\right)\frac{\left\langle {\rho }_{j}\right\rangle }{{\Delta}x}.\end{equation} \tag{ 73 }$$

We computed the momentum field time-correlation of a uniform ideal-gas system with non-Markovian FDDFT, atomistic stochastic simulations (GLE) and theory [equation (73)]. The ideal-gas system is coupled with a thermal bath at fixed temperature T = 1 through a single exponential memory kernel. We analyzed systems characterized by three different memory kernels, corresponding to the under damped (A = 1 and B = −1), the damped (A = 1 and B = −2) and the over damped (A = 1 and B = −4) cases. The domain is discretized in 40 cells of Δx = 100 (in reduced units).

The GLE simulations are run with the gld integrator, which is part of the LAMMPS package [62, 70]. The data are gathered over 500 independent trajectories, each running for 10⁴ steps with a time step Δt = 0.001.

The non-Markovian FDDFT is integrated with the finite-volume scheme discussed above. The data are averaged over 1000 independent trajectories, each running for 10² steps with a time step Δt = 0.1.

Figures 2(a)–(c) show that the momentum time-correlation functions obtained using atomistic, continuum and theoretical approaches converge, within statistical uncertainty, to the same result. It is worth noticing that only exponentially decaying momentum time-correlation functions can be obtained in the Markovian approximation, i.e. the time-correlation function in figure 2(a) can be recovered only by means of a non-Markovian framework. We also report the space correlations between every cell and the middle cell for the overdamped case in figure 2(d). Evidently, our numerical scheme does not introduce any spurious space correlations into the ideal gas.

Zoom In Zoom Out Reset image size

7.2. Non-equilibrium space transition

The main aim here is to demonstrate the significant role of the memory kernel in non-equilibrium processes and to validate the FH algorithm using GLE simulations. We study an ideal-gas system at temperature T = 1 inside the double-well potential:

$$\begin{equation}{U}^{\mathrm{e}\mathrm{x}\mathrm{t}}=C\left[{\left(\frac{x}{{x}_{0}}\right)}^{4}-{\left(\frac{x}{{x}_{0}}\right)}^{2}\right]\end{equation} \tag{ 74 }$$

with C = 5 and x₀ = 200. The time evolution of the system is computed in the underdamped limit, with A₁ = 1 and B₁ = −1, and in the overdamped limit, with A₁ = 1 and B₁ = −4. To force mass transport between the two wells, the system has non-uniform initial conditions with ${\overline{\rho }}^{+}=10$, ${\overline{\rho }}^{-}=5$, where ${\overline{\rho }}^{+}$ refers to the density in the left well and ${\overline{\rho }}^{-}$ to the density in the right well. As in the previous example, the GLE simulations are run with the gld integrator. The data are gathered over 10 independent trajectories, each running for 10⁹ steps with a time step Δt = 0.001.

The non-Markovian FDDFT is integrated with the finite-volume scheme discussed above. The data are averaged over 100 independent trajectories, where each was run for 10⁶ steps with a time step Δt = 0.1.

Figures 3(a) and (b) depict a comparison of the time-evolution of the mean density profiles obtained from GLE and FH in the underdamped and overdamped cases. Despite the presence of a sharp gradient, the numerical scheme is able to accurately reproduce the diffusion of the density step function and the mass transfer between the two wells. The evolution in time of the difference between the total mass contained in the left well m⁺(t) and the mass in the right well m⁻(t) is reported in figure 3(c). The noticeable difference between the overdamped and underdamped cases highlights the importance of the memory effects when modeling out-of-equilibrium systems. This prototypical example suggests that, during phase transitions, but also biological state transitions and chemical reactions, the widely adopted Markovian assumption may not be satisfied.

Zoom In Zoom Out Reset image size

7.3. Memory-driven Turing patterns in binary reacting system

Here we apply our framework to study Turing patterns [43] observed in reaction–diffusion systems. Specifically, our aim is to show how Turing patterns can arise in systems composed of two reacting chemicals characterized by dissimilar memory kernels. To observe Turing patterns in two-component systems one typically needs to use a third substance, which is fixed in space and reversibly bound to only one of the two components. As a result, the effective diffusion coefficient of the binding component is significantly smaller than that of the other components. Here we show that a similar scenario can be realized by using a thermal bath with dissimilar memory kernels for the two reacting components.

Consider the non-linear kinematic chemical-kinetics model proposed by Schnakenberg [44] involving three reactions with chemicals X₁, A, B and X₂:

$$\begin{equation}{X}_{1}\underset{{k}_{-1}}{\overset{{k}_{1}}{{\rightleftharpoons}}}A,\quad B{k}_{2}{\to }{X}_{2},\quad 2{X}_{1}+{X}_{2}{k}_{3}{\to }3{X}_{1},\end{equation} \tag{ 75 }$$

with k₋₁, k₁, k₂ and k₃ representing the reaction-rate constants. This reaction involves an autocatalyst X₁ (it is both product and catalyst) for the same elementary reaction, a reactant X₂ and two auxiliary species A and B, which are later assumed in excess. The prototypical Schnakenberg's reaction model exhibits interesting dynamical behaviors (i.e. steady-state solutions with spatially inhomogeneous patterns and bifurcations [71]). It also has many direct applications in the study of morphogenetic processes involving skin patterns [72]. Let us define the densities of X₁, A, B and X₂ as ρ₁, α, β and ρ₂, respectively. The law of mass action states that the rate of a chemical reaction is proportional to the concentration of the reacting species. This law, when applied to equation (75), gives the time-evolution of the components, namely:

$$\begin{equation}\begin{aligned}\hfill & {w}_{{\rho }_{1}}\left(\boldsymbol{r};t\right)={k}_{-1}\alpha -{k}_{1}{\rho }_{1}\left(\boldsymbol{r};t\right)+{k}_{3}{\rho }_{1}^{2}\left(\boldsymbol{r};t\right){\rho }_{2}\left(\boldsymbol{r};t\right),\hfill \\ \hfill & {w}_{{\rho }_{2}}\left(\boldsymbol{r};t\right)={k}_{2}\beta -{k}_{3}{\rho }_{1}^{2}\left(\boldsymbol{r};t\right){\rho }_{2}\left(\boldsymbol{r};t\right),\hfill \\ \hfill & {w}_{\alpha }\left(\boldsymbol{r};t\right)={k}_{1}{\rho }_{1}\left(\boldsymbol{r};t\right)-{k}_{-1}\alpha ,\hfill \\ \hfill & {w}_{\beta }\left(\boldsymbol{r};t\right)=-{k}_{2}\beta .\hfill \end{aligned}\end{equation} \tag{ 76 }$$

If A and B are in excess, then α and β can be assumed to be constant. It follows that w_α ( r ; t) ∼ w_β ( r ; t) ∼ 0, the so-called 'quasi-steady approximation' in chemical kinetics, and the following system of continuity equations governs the time evolution of the system:

$$\begin{equation}\begin{aligned}\hfill {\partial }_{t}{\rho }_{1}\left(\boldsymbol{r};t\right)=& -{\nabla }_{\boldsymbol{r}}\cdot \left({m}_{s}^{-1}{\boldsymbol{j}}_{1}\left(\boldsymbol{r};t\right)\right)+{w}_{{\rho }_{1}}\left(\boldsymbol{r};t\right)\hfill \\ \hfill {\partial }_{t}{\rho }_{2}\left(\boldsymbol{r};t\right)=& -{\nabla }_{\boldsymbol{r}}\cdot \left({m}_{s}^{-1}{\boldsymbol{j}}_{2}\left(\boldsymbol{r};t\right)\right)+{w}_{{\rho }_{2}}\left(\boldsymbol{r};t\right).\hfill \end{aligned}\end{equation} \tag{ 77 }$$

The following parameters are adopted in our numerical simulations: k₋₁ α = 1.0, k₁ = 0.02, k₂ β = 2.0, k₃ = 10⁻⁶.

Turing patterns are observed when a uniform dynamical system is linearly stable in the absence of diffusion, but becomes linearly unstable in the presence of non-uniform diffusion. In our analysis, we use the one-dimensional overdamped density equation, equation (38), in the weak-noise limit:

$$\begin{equation}\begin{aligned}\hfill {\partial }_{t}{\rho }_{1}\left(x;t\right)& ={d}_{1}{\partial }_{x}^{2}{\rho }_{1}\left(x;t\right)+{w}_{{\rho }_{1}}\left(x;t\right)\hfill \\ \hfill {\partial }_{t}{\rho }_{2}\left(x;t\right)& ={d}_{2}{\partial }_{x}^{2}{\rho }_{1}\left(x;t\right)+{w}_{{\rho }_{2}}\left(x;t\right)\hfill \end{aligned}\end{equation} \tag{ 78 }$$

where we defined ${d}_{1}={k}_{\text{B}}T{D}_{\text{eff},1}{m}_{1}^{-1}$ and ${d}_{2}={m}_{2}^{-1}{k}_{\text{B}}T{D}_{\text{eff},2}$. First, we compute the stationary states of the Schnakenberg model by setting:

$$\begin{equation}\begin{cases}_{{\rho }_{1}}\left(\boldsymbol{r};t\right)=0\quad \hfill \\ {w}_{{\rho }_{2}}\left(\boldsymbol{r};t\right)=0.\quad \hfill \end{cases}\end{equation} \tag{ 79 }$$

A unique uniform stationary state then exists given by:

$$\begin{equation}\begin{aligned}\hfill & {\rho }_{1}^{b}=\frac{{k}_{2}\beta +{k}_{-1}\alpha }{{k}_{1}}\sim 150\hfill \\ \hfill & {\rho }_{2}^{b}=\frac{{k}_{2}\beta {k}_{1}^{2}}{{k}_{3}{\left({k}_{2}\beta +{k}_{-1}\alpha \right)}^{2}}\sim 89.\hfill \end{aligned}\end{equation} \tag{ 80 }$$

The Jacobian matrix for the purely reactive system evaluated at the stationary state, i.e. with components ${J}_{i,j}={\left.\frac{\partial {w}_{{\rho }_{i}}}{\partial {\rho }_{j}}\right\vert }_{{\rho }_{i}^{b}}$, is given by:

$$\begin{equation}\boldsymbol{J}=\left(\begin{matrix}\hfill {k}_{1}\frac{{k}_{2}\beta -{k}_{-1}\alpha }{{k}_{2}\beta +{k}_{-1}\alpha }\hfill & \hfill \frac{{k}_{3}}{{k}_{1}^{2}}{\left({k}_{2}\beta +{k}_{-1}\alpha \right)}^{2}\hfill \\ \hfill \frac{2{k}_{2}\beta {K}_{1}}{{k}_{2}\beta +{k}_{-1}\alpha }\hfill & \hfill -\frac{{k}_{3}}{{k}_{1}^{2}}{\left({k}_{2}\beta +{k}_{-1}\alpha \right)}^{2}\hfill \end{matrix}\right).\end{equation} \tag{ 81 }$$

At the same time, the Jacobian matrix of the reaction–diffusion system is given by J _k = J − [d₁, d₂] I . The reaction parameters chosen in this study satisfy the following conditions

$$\begin{equation}\mathrm{Tr}\enspace J={J}_{11}+{J}_{22}\sim -0.0158{< }0\end{equation} \tag{ 82 }$$

$$\begin{equation}\text{Det}\enspace J={J}_{11}{J}_{22}-{J}_{12}{J}_{21}\sim 0.0004{ >}0\end{equation} \tag{ 83 }$$

to ensure the stability of the purely reactive system, as required for the appearance of Turing patterns (see appendix C for more details).

The effective diffusion parameters d₁ and d₂ are functions of the memory kernels, as shown in section 4. We consider two cases, both with single exponential memory kernels, with parameters A_1,1 = 1.0, A_2,1 = 1.0 and B_2,1 = −20.0. The two cases differ due to B_1,1, which is set to = −0.4 in case (a) and to −0.2 in case (b). The respective values of d₁ are 0.4 and 0.2. Notice that in both cases the Turing instability conditions are satisfied because De tJ_km < 0 (see appendix C for more details). The minimum wavenumber for which Turing instability occurs is given by:

$$\begin{equation}{k}_{m}=\sqrt{\frac{{d}_{1}{J}_{22}+{d}_{2}{J}_{11}}{2{d}_{1}{d}_{2}}}.\end{equation} \tag{ 84 }$$

Such value corresponds to a wave-length L = 2π/k_m .

Figure 4 shows the time-evolution of the density for case (a) obtained from both deterministic and stochastic formulations. The characteristic wavelength for such a system, L ∼ 70, is theoretically predicted by equation (84) and is recovered in both deterministic and stochastic solutions. As expected, for the stochastic evolution the instability is enhanced by the noise, thus the patterns arise sooner. Interestingly, the structure factor reported in figure 5 and discussed in appendix D, exhibits a clear peak at about the predicted wave number. However, in the stochastic case, higher frequency contributions due to the correlated noise are also visible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 6 we report the deterministic and stochastic density evolution for case (b). In this case we also observe the wavelength predicted by the theoretical analysis, L ∼ 50. However, the evolution of this system shows the presence of secondary peaks with low intensity, as also suggested by the corresponding structure factor (figure 7).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Let us conclude this section by emphasizing that the prototypical Schnakenberg model was used to exemplify our FFDFT framework. But the framework is quite general and can be used to model any reaction scheme. Specifically, obtaining the memory kernel for each species (via e.g. molecular dynamics simulations) opens the road to analyzing the influence of interparticle interactions and correlated fluctuations (accounting for thermal bath effects) on pattern formation.

We have derived a theoretical–computational framework for reactive multispecies colloidal systems. Starting from the non-Markovian GLEs describing the motion of the reactive particles, and by appropriately defining the local number and momentum density fields, we obtained the corresponding non-Markovian Navier–Stokes-like system of microscopic FH. Introducing the Helmholtz free-energy functional from classical density-functional theory yields a non-Markovian FDDFT for reacting multispecies fluid systems, the main result of our theoretical work. Introducing an appropriate deconvolution operator enables us to obtain the overdamped limit. This resulted in the non-Markovian and overdamped equation for the density field, a generalisation of FDDFT to include chemical reactions and memory effects. Moreover, in order to make the equations for the non-Markovian FDDFT computationally tractable, we introduced additional field variables whose purpose is twofold: they include the memory term contributions without actually computing any convolution explicitly, and they account for the effects of the time-correlated random term. This additional fields allowed us to develop an extended dynamics framework equivalent to the non-Markovian FDDFT. Additionally, we presented a finite-volume scheme to discretize the non-Markovian FDDFT, and we provided several numerical experiments to exemplify our framework. These include ideal gas in-and-out-of equilibrium, and a two-species reacting system exhibiting memory-driven Turing patterns.

The framework can be adopted to investigate a wide spectrum of physicochemical phenomena in colloidal systems. We have demonstrated that memory effects, through inducing time correlations and variations of the fluid mobility that alter the transition path, can play a crucial role in the dynamical evolution of systems undergoing state transitions. As a result, the theoretical–computational treatment of memory effects requires careful consideration.

We believe that our results can motivate further theoretical–computational studies on reactive colloidal systems and may even be useful in the experimental scrutiny of pattern formation including Turing patterns. There are also a number of interesting questions related to the work presented here. For example, an interesting extension would be to apply the methods proposed here to the study of phase transitions in biological fluids where memory effects should play an important role. Such systems are often characterized by space- and time-scales, which are not accessible with atomistic simulations [7, 8]. Another possible line of future work is to include the effects of geometric confinement [73–75]. We shall examine these and related issues in future studies.

We gratefully acknowledge financial support from the Imperial College (IC) Department of Chemical Engineering PhD Scholarship scheme, European Research Council through Advanced Grant No. 247031 and the Engineering and Physical Sciences Research Council of the UK via Grants No. EP/L027186 and EP/L020564. The computations where performed at the High Performance Computing center of IC. Finally, we are grateful to the anonymous referees for insightful comments and suggestions.

The noise term in equation (10) is ${\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right)={\sum }_{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right){\boldsymbol{R}}_{i,s}\left(t\right)$, where R _i,s (t) has zero mean and correlation function ⟨ R _i,s (t) R _i,s (t')⟩ = k_B T θ _s (t − t'). It follows that the spatiotemporal correlations of η _s can be expressed as:

$$\begin{equation}\langle {\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right){\boldsymbol{\eta }}_{s}\left({\boldsymbol{r}}^{\prime };{t}^{\prime }\right)\rangle =\left\langle \sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)\sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{i,s}\right){\boldsymbol{R}}_{i,s}\left(t\right){\boldsymbol{R}}_{i,s}\left({t}^{\prime }\right)\right\rangle .\end{equation} \tag{ A.1 }$$

Using the property of the Dirac delta function

$$\begin{equation}\sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)\sum _{i=1}^{{N}_{s}\left({t}^{\prime }\right)}\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{i,s}\right)=\sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)=\sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left({\boldsymbol{r}}^{\prime }-{\boldsymbol{r}}_{i,s}\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right),\end{equation} \tag{ A.2 }$$

the spatio–temporal correlations of η_s can be rewritten as:

$$\begin{align}\hfill \langle {\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right){\boldsymbol{\eta }}_{s}\left({\boldsymbol{r}}^{\prime };{t}^{\prime }\right)\rangle & =\left\langle \sum _{i=1}^{{N}_{s}\left(t\right)}\delta \left(\boldsymbol{r}-{\boldsymbol{r}}_{i,s}\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right){\boldsymbol{R}}_{i,s}\left(t\right){\boldsymbol{R}}_{i,s}\left({t}^{\prime }\right)\right\rangle \hfill \\ \hfill & =\left\langle {\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right){k}_{\text{B}}T{\boldsymbol{\theta }}_{s}\left(t-{t}^{\prime }\right)\right\rangle .\hfill \end{align} \tag{ A.3 }$$

Thus, introducing a spatiotemporal noise satisfying ⟨ R _s ( r ; t) R _s ( r '; t')⟩ = k_B T θ _s (t − t')δ( r − r '), the local fluctuation term can be rewritten as

$$\begin{equation}{\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};t\right)=\sqrt{{\rho }_{s}\left(\boldsymbol{r};t\right)}{\boldsymbol{R}}_{s}\left(\boldsymbol{r};t\right).\end{equation} \tag{ A.4 }$$

In this section, we discuss the properties of the stochastic process ξ ( r ; t). First, since all the singularities of $\tilde {\zeta }\left(\omega \right)={\tilde {\chi }}^{-1}\left(\omega \right)$ lie in the lower-half complex ω-plane, then for τ > 0:

$$\begin{align}\hfill \zeta \left(\tau \right)& =\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\omega \enspace \tilde {\zeta }\left(\omega \right){\mathrm{e}}^{\mathrm{i}\omega \tau }=\underset{a\to \infty }{\mathrm{lim}}\frac{1}{2\pi }{\int }_{-a}^{a}\mathrm{d}\omega \enspace \tilde {\zeta }\left(\omega \right){\mathrm{e}}^{\mathrm{i}\omega \tau }\hfill \\ \hfill & ={\frac{1}{2\pi }\oint }_{{C}^{+}}\mathrm{d}\omega \enspace \tilde {\zeta }\left(\omega \right){\mathrm{e}}^{\mathrm{i}\omega \tau }-\underset{a\to \infty }{\mathrm{lim}}\frac{1}{2\pi }{\int }_{\text{arc}\left(a\to -a\right)}\mathrm{d}\omega \enspace \tilde {\zeta }\left(\omega \right){\mathrm{e}}^{\mathrm{i}\omega \tau }=0\hfill \end{align} \tag{ B.1 }$$

where ${\oint }_{{C}^{+}}\mathrm{d}\omega$ indicates the integral over a closed contour C⁺ that goes along the real line from −a to a and then along a semicircle centered at 0 from a to a, while ∫_{arc(a→−a)}dω is the integral along an arc centered at 0 from a to −a. Hence, for t > 0 we can write

$$\begin{equation}{\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};t\right)={\int }_{0}^{+\infty }\zeta \left(t-{t}^{\prime }\right){\boldsymbol{\eta }}_{s}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }={\int }_{-\infty }^{+\infty }\zeta \left(t-{t}^{\prime }\right){\boldsymbol{\eta }}_{s}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }.\end{equation} \tag{ B.2 }$$

Thus, the spatio–temporal correlation function of ξ _s ( r ; t) is given by:

$$\begin{align}\hfill \left\langle {\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};{t}_{1}\right){\boldsymbol{\xi }}_{s}\left({\boldsymbol{r}}^{\prime };{t}_{2}\right)\right\rangle & =\left\langle {\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{1}^{\prime }{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{2}^{\prime }\zeta \left({t}_{1}-{t}_{1}^{\prime }\right){\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};{t}_{1}^{\prime }\right){\left[\zeta \left({t}_{2}-{t}_{2}^{\prime }\right){\boldsymbol{\eta }}_{s}\left({\boldsymbol{r}}^{\prime };{t}_{2}^{\prime }\right)\right]}^{T}\right\rangle \hfill \\ \hfill & ={\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{1}^{\prime }{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{2}^{\prime }\zeta \left({t}_{1}-{t}_{1}^{\prime }\right)\left\langle {\boldsymbol{\eta }}_{s}\left(\boldsymbol{r};{t}_{1}^{\prime }\right){\boldsymbol{\eta }}_{s}^{T}\left({\boldsymbol{r}}^{\prime };{t}_{2}^{\prime }\right)\right\rangle \zeta \left({t}_{2}-{t}_{2}^{\prime }\right)\hfill \end{align} \tag{ B.3 }$$

where we used the fluctuation-dissipation theorem.

Employing the definition of Fourier transform of θ_s and ζ(t), together with equations (44)–(46), we obtain:

$$\begin{align}\hfill \left\langle {\boldsymbol{\xi }}_{s}\left(\boldsymbol{r};{t}_{1}\right){\boldsymbol{\xi }}_{s}\left({\boldsymbol{r}}^{\prime };{t}_{2}\right)\right\rangle & ={k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\enspace \mathbf{I}{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{1}^{\prime }{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{2}^{\prime }\enspace \zeta \left({t}_{1}-{t}_{1}^{\prime }\right)\zeta \left({t}_{2}-{t}_{2}^{\prime }\right){\theta }_{s}\left({t}_{2}^{\prime }-{t}_{1}^{\prime }\right)\hfill \\ \hfill & ={k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\enspace \mathbf{I}{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{1}^{\prime }{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{2}^{\prime }\enspace \zeta \left({t}_{1}-{t}_{1}^{\prime }\right)\frac{1}{2\pi }\hfill \\ \hfill & \quad {\times}{\int }_{-\infty }^{+\infty }\mathrm{d}\omega \enspace {\tilde {\theta }}_{s}\left(\omega \right){\mathrm{e}}^{\mathrm{i}\omega \left({t}_{1}^{\prime }-{t}_{2}^{\prime }\right)}\zeta \left({t}_{2}-{t}_{2}^{\prime }\right)\hfill \\ \hfill & ={k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\enspace \mathbf{I}{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{1}^{\prime }{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{2}^{\prime }\enspace \zeta \left({t}_{1}-{t}_{1}^{\prime }\right)\frac{1}{2\pi }\hfill \\ \hfill & \quad {\times}{\int }_{-\infty }^{+\infty }\mathrm{d}\omega \enspace {\tilde {\theta }}_{s}\left(\omega \right){\mathrm{e}}^{\mathrm{i}\omega \left({t}_{1}^{\prime }-{t}_{2}^{\prime }\right)}\zeta \left({t}_{2}-{t}_{2}^{\prime }\right){\mathrm{e}}^{-\mathrm{i}\omega \left({t}_{1}-{t}_{2}\right)}{\mathrm{e}}^{\mathrm{i}\omega \left({t}_{1}-{t}_{2}\right)}\hfill \\ \hfill & ={k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\enspace \mathbf{I}\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\mathrm{d}\omega {\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{1}^{\prime }\zeta \left({t}_{1}-{t}_{1}^{\prime }\right){\mathrm{e}}^{-\mathrm{i}\omega \left({t}_{1}-{t}_{1}^{\prime }\right)}{\tilde {\theta }}_{s}\left(\omega \right)\hfill \\ \hfill & \quad {\times}{\int }_{-\infty }^{+\infty }\enspace \mathrm{d}{t}_{2}^{\prime }\enspace \zeta \left({t}_{2}-{t}_{2}^{\prime }\right){\mathrm{e}}^{\mathrm{i}\omega \left({t}_{2}^{\prime }-{t}_{2}^{\prime }\right)}{\mathrm{e}}^{\mathrm{i}\omega \left({t}_{1}-{t}_{2}\right)}\hfill \\ \hfill & ={k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\enspace \mathbf{I}\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\mathrm{d}\omega \enspace \tilde {\zeta }\left(\omega \right){\tilde {\theta }}_{s}\left(\omega \right)\tilde {\zeta }\left(-\omega \right){\mathrm{e}}^{\mathrm{i}\omega \left({t}_{1}-{t}_{2}\right)}\hfill \\ \hfill & ={k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\enspace \mathbf{I}\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\mathrm{d}\omega \enspace \tilde {\zeta }\left(\omega \right)2\tilde {\chi }\left(\omega \right)\tilde {\chi }\left(-\omega \right)\tilde {\zeta }\left(-\omega \right){\mathrm{e}}^{\mathrm{i}\omega \left({t}_{1}-{t}_{2}\right)}\hfill \\ \hfill & ={k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\enspace \mathbf{I}\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }\mathrm{d}\omega \enspace 2{\mathrm{e}}^{\mathrm{i}\omega \left({t}_{1}-{t}_{2}\right)}\hfill \\ \hfill & =2{k}_{\text{B}}T{\rho }_{s}\left(\boldsymbol{r};t\right)\delta \left(\boldsymbol{r}-{\boldsymbol{r}}^{\prime }\right)\delta \left({t}_{1}-{t}_{2}\right)\enspace \mathbf{I}.\hfill \end{align} \tag{ B.4 }$$

It follows that ξ _s ( r ; t) is delta-correlated in space and time, and consequently can be easily generated, i.e. as in reference [65].

Turing patterns are observed when a uniform system is stable in the absence of diffusion, but become unstable to perturbation with diffusion. For the analysis, we consider the one-dimensional density equation:

$$\begin{equation}\begin{aligned}\hfill {\partial }_{t}{\rho }_{1}\left(x;t\right)& ={d}_{1}{\partial }_{x}^{2}{\rho }_{1}\left(x;t\right)+{w}_{{\rho }_{1}}\left(x;t\right)\hfill \\ \hfill {\partial }_{t}{\rho }_{2}\left(x;t\right)& ={d}_{2}{\partial }_{x}^{2}{\rho }_{1}\left(x;t\right)+{w}_{{\rho }_{2}}\left(x;t\right)\hfill \end{aligned}\end{equation} \tag{ C.1 }$$

with stationary states given by:

$$\begin{equation}\begin{aligned}\hfill & {w}_{{\rho }_{1}}\left(\boldsymbol{r};t\right)=0\hfill \\ \hfill & {w}_{{\rho }_{2}}\left(\boldsymbol{r};t\right)=0.\hfill \end{aligned}\end{equation} \tag{ C.2 }$$

Hence, we study the linear stability of the system around the stationary state for both the corresponding reacting-only system and for the full system. The time evolution of the infinitesimal perturbations δρ₁ and δρ₂ is given by:

$$\begin{equation}\begin{aligned}\hfill {\partial }_{t}\delta {\rho }_{1}& ={d}_{1}{\partial }_{x}^{2}\delta {\rho }_{1}+{J}_{11}\delta {\rho }_{1}+{J}_{12}\delta {\rho }_{2}\hfill \\ \hfill {\partial }_{t}\delta {\rho }_{2}& ={d}_{2}{\partial }_{x}^{2}\delta {\rho }_{1}+{J}_{21}\delta {\rho }_{1}+{J}_{22}\delta {\rho }_{2}\hfill \end{aligned}\end{equation} \tag{ C.3 }$$

where the coefficients of the Jacobian matrix evaluated at the stationary state ${J}_{i,j}={\left.\frac{\partial {w}_{{\rho }_{i}}}{\partial {\rho }_{j}}\right\vert }_{{\rho }_{i}^{b}}$ are introduced. A particular solution is given by

$$\begin{equation}\left(\begin{matrix}\hfill \delta {\rho }_{1}\hfill \\ \hfill \delta {\rho }_{2}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill \delta {\rho }_{1}^{0}\hfill \\ \hfill \delta {\rho }_{2}^{0}\hfill \end{matrix}\right){\mathrm{e}}^{{\lambda }_{k}t}{\mathrm{e}}^{\mathrm{i}kx}.\end{equation} \tag{ C.4 }$$

Combining equations (C.3) and (C.4), we obtain the eigenvalue problem:

$$\begin{equation}\left({\boldsymbol{J}}_{k}-k\boldsymbol{I}\right)\left(\begin{matrix}\hfill \delta {\rho }_{1}\hfill \\ \hfill \delta {\rho }_{2}\hfill \end{matrix}\right)=\mathbf{0}\end{equation} \tag{ C.5 }$$

where

$$\begin{equation}{\boldsymbol{J}}_{k}=\left(\begin{matrix}\hfill {J}_{1,1}-{d}_{1}{k}^{2}\hfill & \hfill {J}_{1,2j}\hfill \\ \hfill {J}_{2,1}\hfill & \hfill {J}_{2,2}-{d}_{2}{k}^{2}\hfill \end{matrix}\right).\end{equation} \tag{ C.6 }$$

The uniform stationary solution is stable if both eigenvalues λ_k,i at a wave number k have negative real part. Such a condition is verified if:

$$\begin{equation}\mathrm{Tr}\enspace {J}_{k}={J}_{11}+{J}_{22}-\left({d}_{1}+{d}_{2}\right){k}^{2}{< }0\end{equation} \tag{ C.7 }$$

$$\begin{equation}\text{Det}\enspace {J}_{k}=\left({J}_{11}-{d}_{1}{k}^{2}\right)\left({J}_{22}-{d}_{2}{k}^{2}\right)-{J}_{12}{J}_{21}{ >}0.\end{equation} \tag{ C.8 }$$

Turing patterns can arise when a reaction-only system, corresponding to d₁ = d₂ = 0, is stable, but the corresponding reaction–diffusion system is unstable. The reaction-only system is stable when the following two conditions are satisfied:

$$\begin{equation}\mathrm{Tr}\enspace J={J}_{11}+{J}_{22}{< }0\end{equation} \tag{ C.9 }$$

$$\begin{equation}\text{Det}\enspace J={J}_{11}{J}_{22}-{J}_{12}{J}_{21}{ >}0.\end{equation} \tag{ C.10 }$$

The reaction–diffusion system is unstable when at least one of equations (C.7) and (C.8) is not verified. Because of equation (C.9), Tr  J_k < Tr  J < 0, thus equation (C.7) is always verified. Instability then occurs if and only if Det  J_k < 0 for at least a value of k. The minimum value of k for which such a condition occurs is given by:

$$\begin{equation}{k}_{m}=\sqrt{\frac{{d}_{1}{J}_{22}+{d}_{2}{J}_{11}}{2{d}_{1}{d}_{2}}}\end{equation} \tag{ C.11 }$$

which corresponds to a wave number L_k = 2π/k. It directly follows that the corresponding Det  J_km < 0 when:

$$\begin{equation}\text{Det}\enspace {J}_{km}={J}_{11}{J}_{22}-{J}_{12}{J}_{21}-\frac{{d}_{1}{J}_{22}+{d}_{2}{J}_{11}}{4{d}_{1}{d}_{2}}{< }0.\end{equation} \tag{ C.12 }$$

As shown in previous works [57, 65], the structure factor is useful not only to study the stability of the numerical integrator, but also in multi-phase systems to obtain the frequency components characterizing the system. If we consider a periodic domain of volume V, the spatial Fourier transform of the number density for the species s is given by:

$$\begin{equation}{\hat{\rho }}_{s,\mathbf{k}}\left(t\right)=\frac{1}{V}{\int }_{V}{\rho }_{s}\left(\mathbf{r},t\right){\mathrm{e}}^{-\mathrm{i}\mathbf{k}\cdot \mathbf{r}}\mathrm{d}\mathbf{r}.\end{equation} \tag{ D.1 }$$

Thus, the structure factor is defined as the variance of the Fourier transform of the density fluctuations of species s,

$$\begin{equation}S\left(\mathbf{k}\right)=V\langle \delta {\hat{\rho }}_{s,\mathbf{k}}\delta {\hat{\rho }}_{s,\mathbf{k}}^{{\ast}}\rangle \end{equation} \tag{ D.2 }$$

with the Fourier transform of the density fluctuations being $\delta {\hat{\rho }}_{s,\mathbf{k}}={\hat{\rho }}_{s,\mathbf{k}}-\langle {\hat{\rho }}_{s,\mathbf{k}}\rangle$. In case of uniform ideal gases, the structure factor is computed as [57]:

$$\begin{equation}S\left(\mathbf{k}\right)={\overline{\rho }}_{s}/\beta .\end{equation} \tag{ D.3 }$$