Collectivity in diffusion of colloidal particles: from effective interactions to spatially correlated noise

The classical picture of single-particle diffusion, dating back to Einstein and Smoluchowski, describes the observed particle (a tracer) as constantly bombarded by the much smaller particles of its molecular environment [1, 2]. This results in stochastic driving, which is accompanied by Stokesian, hydrodynamic friction that dissipates the energy of the tracer. In the following decades, the major modification to this picture came from Mori and Zwanzig, who provided a formal link between the microscopic interactions in the system and the Langevin dynamics [3]. They have shown that the stochastic force is time-correlated and must be accompanied by friction with a memory kernel. This establishes the formalism of the generalized Langevin equations [3, 4]. Nowadays, time-correlated noise and its consequences, most notably the sub-diffusive processes [5–7], are thoroughly researched phenomena with multiple applications (e.g. [8, 9]).

While the concept of correlations in noise can be extended to the spatial domain, this aspect of stochastic dynamics remains widely unrecognized. In section 2 we provide a summary of phenomena related to what we call spatially correlated noise (SCN). Let us just mention now that SCN manifests itself in such problems as self-assembly and its non-equilibrium dynamics [10], and these problems are increasingly important in many branches of physics.

Few attempts have been aimed at understanding the role of SCN in molecular systems, and even fewer have tried to explain its origin. In our recent article ([10]), we employed the binary mixture theory to show that spatial correlations are inherently present in particle-thermal bath interactions. We have also introduced a new type of thermodynamically consistent SCN-driven Langevin dynamics [10]. This new formalism has been developed for two particles, and in this article we extend it to the multi-particle case. We also introduce the idea of collectivity in diffusion, i.e. the notion that the friction and stochastic force affecting different particles are not independent. Finally, we introduce another approach to the spatially correlated behavior, which is the multi-particle Mori–Zwanzig model. We will show that this model proves complementary to the SCN-based theory and can help us resolve some of its ambiguities.

We will now give a short overview of [10], which sets the context for this current work. Let us assume that, microscopically, the system consists of N tracer particles and the $\tilde{N}$ particles of the thermal bath. The total microscopic Hamiltonian of this system reads:

$$\begin{eqnarray}&&\mathcal{H}={{\mathcal{H}}_{\text{tt}}}+{{\mathcal{H}}_{\text{tb}}}+{{\mathcal{H}}_{\text{bb}}}\end{eqnarray} \tag{ 1 }$$

where ${{\mathcal{H}}_{\text{tt}}}$ contains the microscopic tracer–tracer interactions, ${{\mathcal{H}}_{\text{tb}}}$ is the tracer-bath interaction and ${{\mathcal{H}}_{\text{bb}}}$ is the Hamiltonian of the thermal bath, including its internal interactions. More specifically:

$$\begin{eqnarray}&&{{\mathcal{H}}_{\text{tt}}}=\underset{i}{\overset{N}{\sum}}\,\frac{P_{i}^{2}}{2M}+\underset{i>j}{\overset{N}{\sum}}\,{{U}_{0}}\left({{x}_{i}}-{{x}_{j}}\right)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{\mathcal{H}}_{\text{tb}}}=\underset{i}{\overset{N}{\sum}}\,\underset{j}{\overset{{\tilde{N}}}{\sum}}\,V\left({{x}_{i}}-{{q}_{j}}\right)\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\mathcal{H}}_{\text{bb}}}=\underset{i}{\overset{{\tilde{N}}}{\sum}}\,\frac{p_{i}^{2}}{2m}+\underset{i>j}{\overset{{\tilde{N}}}{\sum}}\,v\left({{q}_{i}}-{{q}_{j}}\right)\end{eqnarray} \tag{ 4 }$$

where $\left({{x}_{i}},{{P}_{i}}\right)$ and $\left({{q}_{i}},\,{{p}_{i}}\right)$ are the phase-space coordinates of the tracers and environment particles, respectively. In this setting, we can determine the effective interactions, i.e. the interactions that in equilibrium are equivalent to the influence of the thermal bath on the tracers [11, 12]. These interactions read:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {} & {{\mathcal{H}}_{\text{eff}}}=\underset{i>j}{\overset{N}{\sum}}\,{{U}_{\text{eff}}}\left({{x}_{i}}-{{x}_{j}}\right) \\ {} & =-\frac{1}{\beta}\ln \left({\int}^{}\left\{\text{d}p\text{d}q\right\}\exp \left(-\beta \left({{\mathcal{H}}_{Tb}}+{{\mathcal{H}}_{\text{bb}}}\right)\right)\right) \end{array}\end{eqnarray} \tag{ 5 }$$

where $\left\{\text{d}p\text{d}g\right\}={\prod}_{i}^{{\tilde{N}}}\text{d}{{p}_{i}}\text{d}{{q}_{i}}$ and $\beta ={{\left({{k}_{\text{B}}}T\right)}^{-1}}$ is the temperature factor. We can now define the deterministic tracer-bath coupling force as:

$$\begin{eqnarray}&&\xi \left({{x}_{i}}\right)=-\underset{j}{\overset{{\tilde{N}}}{\sum}}\,{{\partial}_{{{x}_{i}}}}V\left({{x}_{i}}-{{q}_{j}}\right).\end{eqnarray} \tag{ 6 }$$

We have shown that the covariance of this coupling force reads [10]:

$$\begin{eqnarray}&&\langle \xi (0)\xi (r)\rangle =\langle \underset{\begin{array}{c} j\ne i, \\ k\ne i \end{array}}{\overset{N}{\sum}}\,{{F}_{\text{eff}}}\left({{x}_{i}}-{{x}_{j}}\right){{F}_{\text{eff}}}\left({{x}_{i}}-{{x}_{k}}+r\right)\rangle \end{eqnarray} \tag{ 7 }$$

The above average is taken with respect to the Boltzmann distribution:

$$\begin{eqnarray}&&{{P}_{\text{B}}}\left({{x}_{1}},...,{{x}_{N}}\right)={{\mathcal{N}}^{-1}}\exp \left(-\beta \underset{i>j}{\sum}\,U\left({{x}_{i}}-{{x}_{j}}\right)\right)\end{eqnarray} \tag{ 8 }$$

where $U(r)={{U}_{0}}(r)+{{U}_{\text{eff}}}(r)$ and $\mathcal{N}$ is the normalization constant. Then, the spatial correlation function of the tracer-bath coupling force reads:

$$\begin{eqnarray}&&h(r)=\langle \xi (0)\xi (r)\rangle /\langle {{\xi}^{2}}\left({{x}_{i}}\right)\rangle.\end{eqnarray} \tag{ 9 }$$

Thus, these correlations prove to be directly related to the correlations of the effective forces. This result provides the microscopic foundation for SCN, which relates it to the wide class of soft matter phenomena. Furthermore, we introduce the Langevin dynamics of two tracers driven by the Gaussian SCN, in which the noise correlation function is given by h(r). We have found that the thermodynamic consistency (i.e. recovering the Boltzmann distribution (8) in the steady state) requires that the friction coefficient for the relative variable $r={{x}_{2}}-{{x}_{1}}$ is spatially variant. More precisely, we have derived the following over-damped equation of motion:

$$\begin{eqnarray}&&K(r)\overset{\centerdot}{{r}}\,=2F(r)+\sqrt{2}\sigma {{g}_{-}}(r)\eta \end{eqnarray} \tag{ 10 }$$

where $F(r)={{F}_{0}}(r)+{{F}_{\text{eff}}}(r)$, ${{g}_{-}}(r)=\sqrt{1-h(r)}$, η is the uncorrelated Gaussian noise, h(r) is the correlation function given by (9), σ is the noise amplitude and

$$\begin{eqnarray}&&K(r)=\frac{{{g}_{-}}(r){{\text{e}}^{-\beta U(r)}}}{\frac{1}{\gamma}-\frac{2}{{{\sigma}^{2}}}{\int}_{r}^{+\infty}\text{d}{{r}^{\prime}}\frac{F\left({{r}^{\prime}}\right)}{g\left({{r}^{\prime}}\right)}{{\text{e}}^{-\beta U\left({{r}^{\prime}}\right)}}}\end{eqnarray} \tag{ 11 }$$

is the spatially variant friction coefficient (SVFC). The application of this SCN-driven Langevin dynamics to the system of two charged spheres and counter-ions leads to several non-equilibrium effects, e.g. the emergence of the friction-less regime and the transient attraction effect [10].

These results settle the technical core of our theory, and in this article we want to propose its generalizations as well as discuss some of its additional aspects. First, in section 2 we review numerous contexts in which SCN is encountered, speculating on the applicability of our theory. In section 3 we explain the notion of collectivity in diffusion, which is the major consequence of [10], and in section 4 we extend our formalism to systems with multiple (N  >  2) tracers. Furthermore, we will confront our theory with another approach, which is the multi-particle Mori–Zwanzig model. In [10], we have shown that the SVFC tends to zero as the particles get close to each other. This behavior arises in the over-damped theory, which indicates its breakdown. Thus, an insight from the inertial dynamics is desired, and we propose the multi-tracer Mori–Zwanzig model as a solution to this problem. The Mori–Zwanzig model is also the simplest, microscopically consistent approach to collectivity in diffusion, and we show that it provides a complementary perspective for our SCN-driven Langevin dynamics. However, we will also show that the Mori–Zwanzig approach is severely limited in its predictions and, in fact, that it can only serve as an auxiliary short-distance approximation. The Mori–Zwanzig model is solved in section 5 and the comparison to SCN-driven dynamics is given in section 6.

SCN is a random disturbance statistically characterized by a certain correlation length-scale, i.e. the stochastic forces 'felt' by two nearby particles are similar in terms of their direction and amplitude, though the pattern of forces at a greater length-scale is completely random. Such a pattern can also randomly evolve in time and if its memory is short, this system might be thought of as purely SCN-driven.

One prominent class of phenomenon affected by spatially correlated behavior is dense molecular systems, especially in the glassy state [13]. It has been shown via both simulations and experiments that the velocities of particles in the glassy regime are spatially correlated [14, 15], the rearrangement of molecules is highly cooperative [16] and that clusters of particles moving together are also observed [17]. These features strongly resemble SCN. Yet another discipline in which SCN occurs is in intracellular biophysics, e.g. it has been shown that the flows of cytoplasm inside a cell are spatially correlated [18]. Recent research has proven that most of the diffusive transport in cytoplasm results from the constant stirring by molecular motors [19] and other active components [20]. Since the random flows caused by stirring are spatially extensive, the intracellular diffusion should be considered as being SCN-driven. This is partially related to the field of active matter in which the self-propelling particles are also known to produce SCN, e.g. in the interaction between the polymer chain and active swimmers [21] or when active particles form clusters [22]. Spatially correlated behavior is also encountered in fluids, e.g. in sedimentation experiments (as a result of hydrodynamic interactions, e.g. [23]) or in turbulent flows [24].

While the systems we mention span a wide range of physical phenomena, they share a common approach in modeling, which is the Langevin dynamics. However, while the spatial correlations seem to play a crucial role in these systems, the ordinary Langevin dynamics neglects these correlations completely. Thus, an understanding of these systems might be improved with the perspective of SCN. However, SCN-driven dynamics, which we discuss here, is based on the equilibrium correlation function, so it is only applicable in the close to equilibrium regime. This means that far-from-equilibrium systems—especially the flows and active particles—lie beyond its current scope. However, the general observation that the SCN requires SVFC should a specific stationary distribution be recovered, might remain valid for these systems. The adaptation of our approach would require the adequate non-equilibrium distributions to be identified. Glassy systems are a more promising field, since they are the limiting case of the equilibrium colloidal systems, for which our theory is applicable.

Two-component (binary) colloids are a class of systems for which our SCN-driven Langevin dynamics has been derived. In these systems, self-organization is induced by effective interactions [11, 12] as given by (5). The relation (7) shows that the presence of effective interactions entails the presence of spatial correlations in tracer-thermal bath interactions. This suggests that self-assembly might be dynamically perceived as SCN-driven. The simplest example of this phenomenon is, e.g. the clustering of colloidal spheres due to the presence of smaller particles [25]. In this example, described on the statistical level by the renowned Asakura–Oosawa model [26], the bigger spheres do not interact in the long-range fashion, but they are kept together by interaction with the environment. Although for the hard spheres one might explain this effect by a simple imbalance in the number of molecular collisions [12], there are more subtle aspects of this phenomenon, e.g. the multiple energy minima of this interaction in a dense environment [25] or the long-range 'attraction-through-repulsion'/'repulsion-through-attraction' effects for screened-charged spheres [11, 27], which suggests cooperativeness in the influence of the environment. In fact, in light of our SCN-driven Langevin dynamics, any type of self-organization driven by effective interactions (e.g. the separation of polymer blends [28], like-charge attraction [29] etc) can be interpreted as a manifestation of SCN. We should also mention that many of the biological self-assembly mechanisms (e.g. [30–32]) belong to this category, since molecular crowding and excluded volume effects are important factors in the intracellular environment. However, we should mention that we are not aware of any direct experimental observations of SCN in colloidal systems. The reason is that most of the experimental techniques can directly measure the positions of particles, but deducing random forces usually requires one to pre-assume some microscopic model. Thus, the spatial correlations are usually associated with deterministic rather than stochastic forces.

Phenomena qualitatively similar to colloidal self-assembly are also encountered in plasma physics. The electro-magnetic field is per se a spatially extensive and fluctuating entity, and a few attempts have been made to incorporate this fact into the theoretical description of the diffusion in plasma (e.g. [33, 34]). However, these considerations address solely the diffusion of a single particle. On the other hand, the two-component plasma can behave in a similar way to the binary colloids, and analogous separation phenomena emerge in both simulations and experiments [35, 36]. Some theoretical attempts (e.g. [37]) describe these effects from the perspective of effective interactions.

The given examples show that SCN emerges in many different branches of physics. Despite this fact, the literature regarding the influence of SCN on molecular systems is rather scarce. One example is the research on SCN-driven single file diffusion, i.e. the effectively one-dimensional diffusion of impenetrable particles in narrow pores, which leads to sub-diffusive dynamics [40]. Another piece of research utilizes SCN as a stimulus for a neural network [41], concluding that it is a crucial factor in obtaining the synchronized output. Yet further research has been carried out by the current authors, who have analyzed the influence of SCN on a model polymeric chain. We have found that SCN can cause the spontaneous unfolding effect and synchronization of monomer motion [38, 39], possibly leading to non-Gaussian chain statistics [28]. This last idea finds qualitative confirmation in the research on the influence of active swimmers on a polymer chain by Shin et al [21].

Finally, from a purely theoretical perspective, a distinct class of stochastic differential equations with multiplicative noise can be related to the SCN-driven system. This comes from the properties of the correlated Gaussian noise. Let $\xi \left({{x}_{i}}\right)$ be the Gaussian random variable that satisfies:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle \xi \left({{x}_{i}}\right)\rangle =0 \\ \langle {{\xi}^{2}}\left({{x}_{i}}\right)\rangle ={{\sigma}^{2}} \\ \langle \xi (r)\xi (0)\rangle /\langle {{\xi}^{2}}\left({{x}_{i}}\right)\rangle =h(r) \end{array}\end{eqnarray} \tag{ 12 }$$

where h(r) is now some spatial correlation function. Knowing that the sum of the Gaussian variables is also a Gaussian variable, one can write [10]:

$$\begin{eqnarray}&&\xi \left({{x}_{1}}\right)\pm \xi \left({{x}_{1}}+r\right)=\sqrt{2}\sigma \sqrt{1\pm h(r)}\eta =\sqrt{2}\sigma {{g}_{\pm }}(r)\eta \end{eqnarray} \tag{ 13 }$$

where η is now the uncorrelated Gaussian noise. This shows that the linear combination of additive, spatially correlated noise terms can be translated into one, uncorrelated, but multiplicative noise term. In what follows, any stochastic differential equation of the form

$$\begin{eqnarray}&&\gamma \overset{\centerdot}{{r}}\,=2F(r)+\sqrt{2}\sigma {{g}_{-}}(r)\eta \end{eqnarray}$$

where F(r) is an odd function, can be turned into the over-damped dynamics of two particles driven by SCN. Namely, given that $r={{x}_{2}}-{{x}_{1}}$, their equations of motion read:

$$\begin{eqnarray}&&\gamma {{\overset{\centerdot}{{x}}\,}_{i}}=F\left({{x}_{i}}-{{x}_{j}}\right)+\xi \left({{x}_{i}}\right)\end{eqnarray}$$

where $\langle \xi \left({{x}_{i}}\right)\xi \left({{x}_{j}}\right)\rangle /{{\sigma}^{2}}=1-g_{-}^{2}(r)$. From this perspective references [42, 43], for example, can be considered as the analysis of the over-damped oscillators driven by SCN with $\langle \xi \left({{x}_{i}}\right)\xi \left({{x}_{j}}\right)\rangle /{{\sigma}^{2}}=1-{{r}^{2\alpha}}$. It is shown there that such an SCN induces the localization of particles, and the probability distribution for r is discontinuous in α.

Unfortunately, the disadvantage of the existing research on SCN-driven systems is that it makes use of the constant friction coefficient. The crucial observation from our theory is that one cannot achieve thermodynamic consistency without the SVFC given by (11) [10]. This means that none of these results can be related to the equilibrium conditions in a simple way. A similar observation regarding SCN appears independently in [44]. The notion of SVFC itself is also not entirely new and can be encountered in the research on diffusion in viscosity landscapes, e.g. in [45]. Interestingly, in this context the authors intuitively conclude that the system must be driven by multiplicative noise. This suggests a weak analogy between their approach and our SCN-driven dynamics, i.e. one could interpret our SVFC as diffusion in some effective viscosity landscape. In this context, it is particularly intriguing whether the sub-diffusive behavior predicted in [45] could emerge in our approach. Unfortunately, as our equations of motion are highly nonlinear and the dependence between K(r), F(r) and g(r) is strictly fixed (unlike in [45]), finding an exactly solvable toy-model for our theory is still an open problem.

Finally, we should also mention that systems with SCN are often subjected to some significant temporal correlations. Unfortunately, the problem of spatio-temporal dynamics is extremely challenging and reaches beyond the scope of this paper. In this context, our SCN-based approach should be perceived as a step towards the ultimate spatio-temporal solution.

Let us now discuss a major consequence of relation (7), which is collectivity in diffusion. When we describe the diffusion of the group of tracers driven by the common thermal bath, this resembles the situation of the binary mixture. In other words, we divide the entire system into two species of particles (tracers and the particles of a thermal bath) and try to replace one of them with some effective impact on the other. The difference is that in soft matter we are usually interested in the stationary behavior of the system, and in diffusion we ask for the dynamics. Nevertheless, the traced-out species of particles is always the source of effective forces and this fact should be reflected by the properties of the thermal bath.

First, let us discuss the most basic approach. The mean value of any observable O({x}) dependent on the positions of tracers (but not thermal bath particles) can be calculated according to the expression:

$$\begin{eqnarray}&&\langle O\left(\{x\}\right)\rangle ={{\mathcal{N}}^{-1}}{\int}^{}\left\{\text{d}P\text{d}x\right\}O\left(\{x\}\right){{\text{e}}^{-\beta \left({{H}_{\text{tt}}}+{{H}_{\text{eff}}}\right).}}\end{eqnarray} \tag{ 14 }$$

The effective interactions enter this formula just like another regular potential. Thus, the simplest and very common approach is to describe a dynamical system in terms of the ordinary over-damped Langevin equations, with the effective interactions included as an additional force, namely:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \gamma {{\overset{\centerdot}{{x}}\,}_{i}}=\underset{j}{\overset{N}{\sum}}\,{{F}_{0}}\left({{x}_{i}}-{{x}_{j}}\right)+\underset{j}{\overset{N}{\sum}}\,{{F}_{\text{eff}}}\left({{x}_{i}}-{{x}_{j}}\right)+{{\eta}_{i}}(t) \\ \langle {{\eta}_{i}}(t){{\eta}_{j}}\left({{t}^{\prime}}\right)\rangle ={{\sigma}^{2}}{{\delta}_{ij}}\delta \left(t-{{t}^{\prime}}\right). \end{array}\end{eqnarray} \tag{ 15 }$$

Here, the friction coefficient is constant and the noise is uncorrelated. Writing the complementary Fokker–Planck equation for this system and solving it in the stationary state, one obtains the Boltzmann distribution (8). This means that one can use (15) to obtain the averaged-out observables $\langle O\left(\{x\}\right)\rangle$ in the equilibrium. However, this does not guarantee that the trajectories generated by (15) are close to the actual physical trajectories—it only ensures the statistical agreement. Another problem is that in this model the average values of the observables related to the thermal bath are directly affected by the choice of uncorrelated Gaussian noise in (15).

The equation (7) suggests that the role of effective interactions in the dynamics, even in equilibrium, is more complicated than being simply another potential. In light of (7), effective interactions might be understood as dynamically induced by the correlated behavior of the thermal bath. Intuitively, one expects then that the friction and noise (which replace the deterministic tracer-bath coupling) should also reproduce these spatial correlations. Satisfying this intuition alongside the requirement of thermodynamic consistency leads to the inclusion of SCN and SVFC. However, SCN and SVFC make the diffusion a collective phenomenon, i.e. the thermal noise affecting different particles is not independent, and SVFC takes into account that one particle disturbs the environment perceived by the other ones. In [10] this manifests itself via the additional response forces that arise when the absolute positions of particles are considered. Obviously, the thermodynamically consistent Langevin equations with SCN [10] are not perfectly equivalent to the microscopic, deterministic theory. One reason for this is, e.g. the choice of Gaussian noise, and currently we cannot assess how different it is from the actual distribution of the coupling forces. However, via SCN and SVFC we are able to transfer more properties from the microscopic level into the stochastic dynamics, e.g. some microscopic self-assembly mechanisms. This is because both SCN and SVFC depend explicitly on the spatial correlation function h(r), which is sensitive to the thermodynamical state of the system [46]. This opens the possibility of formulating a single stochastic model adequate for a few different regimes.

Another aspect of this is the non-equilibrium regime. The agreement between the microscopic theory (given by the Boltzmann distribution) and the stochastic process approach can be ensured only in equilibrium conditions. However, the Langevin dynamics also provides an insight into the non-equilibrium regime, if it is started from non-equilibrium initial conditions. In this case, one expects that the Langevin equations of motion resemble the microscopic equations of motion closely enough, so the inference about the system evolution is allowed. From this perspective, by incorporating the spatial correlations in the noise, we construct a model which is closer to the actual physics than the non-correlated one. Thus, using SCN is possibly more accurate in the non-equilibrium regime than neglecting the correlations entirely. Obviously, since our model makes use of the equilibrium correlation function, it is still limited to the close-to-equilibrium regime. However, even under this constraint, our formalism can reveal new collective, non-equilibrium effects such as, e.g. the transient attraction [10].

We will now propose the extension of our formalism to one-dimensional, multi-particle systems (N  >  2). Let us introduce the random vector:

$$\begin{eqnarray}&&{{\vec{\xi}}^{T}}=\left({{\xi}_{1}},\ldots,{{\xi}_{N}}\right)\end{eqnarray} \tag{ 16 }$$

whose components satisfy:

$$\begin{eqnarray}&&\langle {{\xi}_{i}}(t)\rangle =0\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\langle {{\xi}_{i}}(t){{\xi}_{j}}\left({{t}^{\prime}}\right)\rangle ={{\sigma}^{2}}h\left({{x}_{i}}-{{x}_{j}}\right)\delta \left(t-{{t}^{\prime}}\right)\end{eqnarray} \tag{ 18 }$$

so $\vec{\xi}$ contains the correlated Gaussian variables. We shall omit the argument t in the notation if it is not necessary. We can introduce the following correlation matrix:

$$\begin{eqnarray}&&\langle \vec{\xi}{{\vec{\xi}}^{T}}\rangle =H\end{eqnarray} \tag{ 19 }$$

where the entries are ${{H}_{ij}}={{\sigma}^{2}}h\left({{x}_{i}}-{{x}_{j}}\right)$. H is a symmetric, real and positive definite matrix, so its eigen-decomposition can be written in the following form:

$$\begin{eqnarray}&&H=Q{{ \Lambda }^{2}}{{Q}^{T}}=Q \Lambda {{Q}^{T}}Q \Lambda {{Q}^{T}}={{G}^{2}}\end{eqnarray} \tag{ 20 }$$

where $\Lambda$ is the diagonal matrix such that ${{ \Lambda }^{2}}$ contains the eigenvalues of H and the matrix Q is orthonormal. One can observe that the matrix $G=Q \Lambda {{Q}^{T}}={{G}^{T}}$ is also symmetric, so, in fact:

$$\begin{eqnarray}&&H={{G}^{2}}=G{{G}^{T}}.\end{eqnarray} \tag{ 21 }$$

The important fact is that we can generate the vector of correlated variables $\vec{\xi}$ from the linear combination of the uncorrelated Gaussian variables with the aid of G, namely:

$$\begin{eqnarray}&&\vec{\xi}=G\vec{\eta}\end{eqnarray} \tag{ 22 }$$

where:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{{\vec{\eta}}}^{T}}=\left({{\eta}_{1}},\ldots,{{\eta}_{N}}\right) \\ \langle {{\eta}_{i}}\rangle =0 \\ \langle {{\eta}_{i}}(t){{\eta}_{j}}\left({{t}^{\prime}}\right)\rangle ={{\delta}_{ij}}\delta \left(t-{{t}^{\prime}}\right). \end{array}\end{eqnarray} \tag{ 23 }$$

One can check that:

$$\begin{eqnarray}&&\langle \left(G\vec{\eta}\right){{\left(G\vec{\eta}\right)}^{T}}\rangle =G\langle \vec{\eta}{{\vec{\eta}}^{T}}\rangle {{G}^{T}}=H.\end{eqnarray} \tag{ 24 }$$

From the two-particle case we know that the use of SCN requires SVFC. We also know that in absolute variables, the SVFC is translated into the friction coefficient and the response forces [10]. Therefore, we postulate the following generalization from the two-tracer model to the case of the N-tracer system:

$$\begin{eqnarray}&&\underset{i}{\overset{N}{\sum}}\,{{K}_{ij}}{{\overset{\centerdot}{{x}}\,}_{j}}=\underset{j}{\overset{N}{\sum}}\,F\left({{x}_{i}}-{{x}_{j}}\right)+{{\xi}_{j}}\end{eqnarray} \tag{ 25 }$$

where K_ii is the SVFC for the ith tracer and ${{K}_{ij\ne i}}$ are the response forces. We assume that K_ij are the functions of ${{x}_{1}},\ldots,{{x}_{2}}$, but we will omit the arguments for a more compact notation. Using (22), we can rewrite the equations of motion in the matrix form:

$$\begin{eqnarray}&&K\overset{\centerdot}{{{\vec{x}}}}\,=\vec{F}+G\vec{\eta}\end{eqnarray} \tag{ 26 }$$

where the ith component of the vector $\vec{F}$ reads:

$$\begin{eqnarray}&&{{F}_{i}}=\underset{j}{\overset{N}{\sum}}\,F\left({{x}_{i}}-{{x}_{j}}\right).\end{eqnarray} \tag{ 27 }$$

Our goal is to determine the entries of the matrix K. We begin with putting (26) into the following form:

$$\begin{eqnarray}&&\overset{\centerdot}{{{\vec{x}}}}\,={{K}^{-1}}G{{G}^{-1}}\vec{F}+{{K}^{-1}}G\vec{\eta}\end{eqnarray} \tag{ 28 }$$

where we assume that K and G are invertible. Let us now introduce the auxiliary matrix S:

$$\begin{eqnarray}&&S={{K}^{-1}}G\end{eqnarray} \tag{ 29 }$$

so our problem reads:

$$\begin{eqnarray}&&\overset{\centerdot}{{{\vec{x}}}}\,=S{{G}^{-1}}\vec{F}+S\vec{\eta}.\end{eqnarray} \tag{ 30 }$$

Applying the Stratonovich interpretation for this system, its stationary Fokker–Planck equation reads:

$$\begin{eqnarray}&&0=\underset{i}{\overset{N}{\sum}}\,{{\partial}_{{{x}_{i}}}}\left[\underset{j}{\overset{N}{\sum}}\,{{S}_{ij}}\left(\underset{k}{\overset{N}{\sum}}\,{{\tilde{g}}_{jk}}{{F}_{k}}P+\frac{1}{2}\underset{k}{\overset{N}{\sum}}\,{{\partial}_{{{x}_{k}}}}\left({{S}_{kj}}P\right)\right)\right]\end{eqnarray} \tag{ 31 }$$

where ${{\tilde{g}}_{ij}}={{\left({{G}^{-1}}\right)}_{ij}}$ and ${{P}_{s}}=P\left({{x}_{1}},\ldots,{{x}_{N}}\right)$ is the stationary probability distribution describing the system. Let us now demand that $P={{P}_{\text{B}}}$, i.e. that our system tends to the Boltzmann distribution in the stationary state. Employing the definition (8) of P_B in (31) we obtain the equation:

$$\begin{eqnarray}&&0=\underset{i}{\overset{N}{\sum}}\,{{\partial}_{{{x}_{i}}}}\left[{{P}_{\text{B}}}\underset{j}{\overset{N}{\sum}}\,{{S}_{ij}}\underset{k}{\overset{N}{\sum}}\,\left({{\tilde{g}}_{jk}}{{F}_{k}}+\frac{\beta}{2}{{F}_{k}}{{S}_{kj}}+\frac{1}{2}{{\partial}_{{{x}_{k}}}}{{S}_{kj}}\right)\right].\end{eqnarray} \tag{ 32 }$$

At this point we should specify the boundary conditions. In general, one might observe that when particles are far from each other, the correlation matrix H becomes diagonal (${{H}_{ij}}\to {{\sigma}^{2}}{{\delta}_{ij}}$) and so does ${{G}_{ij}}\to \sigma {{\delta}_{ij}}$. In this situation, we also expect that ${{K}_{ij}}\to \gamma {{\delta}_{ij}}$, so, eventually, the ${{S}_{ij}}\to \frac{\sigma}{\gamma}{{\delta}_{ij}}$. However, even under these conditions we still have some significant freedom in the structure of matrix S. One particularly simple choice is to demand for S to be in fact diagonal, so ${{S}_{ij}}={{S}_{ii}}{{\delta}_{ij}}$. One can check now that:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{S}^{2}} & =\langle \left(S\vec{\eta}\right){{\left(S\vec{\eta}\right)}^{T}}\rangle =\left({{K}^{-1}}\right)G{{G}^{T}}{{\left({{K}^{-1}}\right)}^{T}} \\ {} & ={{K}^{-1}}H{{\left({{K}^{-1}}\right)}^{T}} \end{array}\end{eqnarray} \tag{ 33 }$$

and for the diagonal S this means that K⁻¹ is the transformation which diagonalizes H, i.e. K⁻¹  =  Q.

Under this choice, the Fokker–Planck equation simplifies into:

$$\begin{eqnarray}&&0=\underset{i}{\overset{N}{\sum}}\,{{\partial}_{{{x}_{i}}}}\left[{{P}_{\text{B}}}{{S}_{ii}}\left(\underset{k}{\overset{N}{\sum}}\,{{\tilde{g}}_{ik}}{{F}_{k}}+\frac{\beta}{2}{{F}_{i}}{{S}_{ii}}+\frac{1}{2}{{\partial}_{{{x}_{i}}}}{{S}_{ii}}\right)\right].\end{eqnarray} \tag{ 34 }$$

Taking into account the boundary condition, the solution for each S_ii reads:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{S}_{ii}}={{\text{e}}^{\beta {\sum}_{m\ne i}^{N}U\left({{x}_{i}}-{{x}_{m}}\right)}} \\ \,\,\,\,\,\,\,\,\,\,\,\,\times \,\left(\frac{\sigma}{\gamma}-{\int}_{{{x}_{i}}}^{+\infty}\text{d}{{y}_{i}}\underset{k}{\overset{N}{\sum}}\,{{\tilde{{g}}\,}_{ik}}{{F}_{k}}{{\text{e}}^{-\beta {\sum}_{n\ne i}^{N}U\left(\,{{y}_{i}}-{{x}_{n}}\right)}}\right). \end{array}\end{eqnarray} \tag{ 35 }$$

Note that in the product of functions ${{\tilde{g}}_{ik}}{{F}_{k}}$ the argument x_i is replaced by the integration variable y_i.

Having found the matrix S, one obtains the matrix of SFVC and the response forces via the transformation K  =  GS⁻¹. One can see that certain properties from the two-particle model might transfer into the multi-particle model. In particular, it might be possible that for a certain choice of $F\left({{x}_{i}}-{{x}_{j}}\right)$ and $h\left({{x}_{i}}-{{x}_{j}}\right)$ there exists such a set of positions $x_{1}^{0},\ldots,x_{N}^{0}$ that ${{S}_{ii}}\left(x_{1}^{0},\ldots,x_{N}^{0}\right)=0$. In this case, the SVFC for the ith particle would become singular and the effects analyzed in [10] might arise. However, this time we deal with the collective behavior in a sense that all particles have an influence on the value of S_ii. Practically, when both h(r) and F(r) have a finite range, only the interactions with a limited number of neighbors are significant. It is then possible to truncate the sum ${\sum}_{k}^{N}{{\tilde{g}}_{ik}}{{F}_{k}}$ in (35), e.g. to k  =  i  −  1,i,i  +  1. This means that S_ii  =  0 for a single particle might be ensured by the proper choice of $x_{i-1}^{0},x_{i}^{0},x_{i+1}^{0}$, and so it could be possible to find such a set of $\left\{x_{i}^{0}\right\}$ that satisfies S_ii  =  0 for all i simultaneously. This would indicate a global change in the characteristics of the system; however, obtaining this solution is numerically difficult and we shall not purse it in this paper.

Finally, we must comment that the solution (35) is not unique, and it is possible to solve (34) e.g. without demanding that S be diagonal. We should also denote that (35) leads to the non-symmetric matrix $K\ne {{K}^{T}}$. Whether the response forces should be reciprocal or not (${{K}_{ij}}={{K}_{ji}}$) should be a subject of further research. However, imposing the constraints that ensure the symmetry of K on (34) complicates the problem significantly and the solution in this case is also unknown.

The SCN-driven Langevin dynamics which we have discussed avoids direct microscopic considerations. For this reason it has a few ambiguities that cannot be addressed within its own framework. The first thing is that one cannot tell whether ${{F}_{\text{eff}}}(r)$ should be present in (10) and (25) explicitly, or whether it should manifest itself only indirectly, via the SCN and SVFC. In fact, it is equally possible to construct a dynamics similar to (10), but without the ${{F}_{\text{eff}}}(r)$ term, and our approach gives no criterion to distinguish between them. Another issue is that in equation (10), as $r\to 0$ we also have ${{g}_{-}}(r)\to 0$ and $K(r)\to 0$. This means that our over-damped theory leads to the essentially frictionless (hence inertial) dynamics for small r. It needs to be answered whether this is an artifact of the theory or a valid prediction. Finally, our dynamics incorporates SCN only, while completely abandoning the temporal correlations; one might ask then if there is a way to circumvent this limitation. We will answer these questions by undertaking a completely alternative approach, which is the multi-tracer variant of the renowned Mori–Zwanzig theory [3, 4].

In this approach, we assume that N tracers are coupled to the common heat-bath of $\tilde{N}$ oscillators. The Hamiltonian of our system is given by the equations (1)–(4), where we specify:

$$\begin{eqnarray}&&V\left({{x}_{i}}-{{q}_{n}}\right)=\frac{mw_{in}^{2}}{4}{{\left({{q}_{n}}-\frac{2{{c}_{in}}}{mw_{in}^{2}}{{x}_{i}}\right)}^{2}},\end{eqnarray} \tag{ 36 }$$

${{U}_{0}}\left({{x}_{i}}-{{x}_{j}}\right)$ is some arbitrary interaction and $v\left({{q}_{i}}-{{q}_{j}}\right)=0$. Adopting the most general approach we assume that c_in and w_in might have arbitrary values. Following the approach of Mori and Zwanzig [3, 4], we can write the equations of motion for this system:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} M{{\overset{\centerdot}{{x}}\,}_{i}}={{P}_{i}} \\ {{\overset{\centerdot}{{P}}\,}_{i}}=-\underset{j}{\overset{N}{\sum}}\,{{\partial}_{{{x}_{i}}}}{{U}_{0}}\left({{x}_{i}}-{{x}_{j}}\right)+\underset{n}{\overset{{\tilde{N}}}{\sum}}\,{{c}_{in}}\left({{q}_{n}}-\frac{2{{c}_{in}}}{mw_{in}^{2}}{{x}_{i}}\right) \\ m{{\overset{\centerdot}{{q}}\,}_{n}}={{p}_{n}} \\ {{\overset{\centerdot}{{p}}\,}_{n}}=-m{{q}_{n}}\underset{i}{\overset{N}{\sum}}\,\frac{w_{in}^{2}}{2}+\underset{i}{\overset{N}{\sum}}\,{{c}_{in}}{{x}_{i}}. \end{array}\end{eqnarray} \tag{ 37 }$$

Let us denote:

$$\begin{eqnarray}&&\omega _{n}^{2}=\underset{i}{\overset{N}{\sum}}\,\frac{w_{in}^{2}}{2}.\end{eqnarray} \tag{ 38 }$$

The solution for q_n(t) reads [4]:

$$\begin{eqnarray}\begin{array}{*{35}{l}} {{q}_{n}}(t)= & {{q}_{n}}\left({{t}_{0}}\right)\cos \left({{\omega}_{n}}\left(t-{{t}_{0}}\right)\right)+\frac{{{p}_{n}}\left({{t}_{0}}\right)}{m{{\omega}_{n}}}\sin \left({{\omega}_{n}}\left(t-{{t}_{0}}\right)\right) \\ {} & +\underset{i}{\overset{N}{\sum}}\,\frac{{{c}_{in}}}{m{{\omega}_{n}}}{\int}_{{{t}_{0}}}^{t}\text{d}s\sin \left({{\omega}_{n}}(t-s)\right){{x}_{i}}(s). \end{array}\end{eqnarray} \tag{ 39 }$$

The classical move now is to apply the integration by parts in the convoluted term:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{i}{\overset{N}{\sum}}\,\frac{{{c}_{in}}}{m{{\omega}_{n}}}{\int}_{{{t}_{0}}}^{t}\text{d}s\sin \left({{\omega}_{n}}\left(t-{{t}_{0}}\right)\right){{x}_{i}}(s)=\underset{i}{\overset{N}{\sum}}\,\frac{{{c}_{in}}}{m\omega _{n}^{2}}{{x}_{i}}(t)-\underset{i}{\overset{N}{\sum}}\,\frac{{{c}_{in}}}{m\omega _{n}^{2}}{{x}_{i}}\left({{t}_{0}}\right)\cos \left({{\omega}_{n}}\left(t-{{t}_{0}}\right)\right) \\ \qquad \qquad \qquad \qquad \qquad \qquad \ -\,\underset{i}{\overset{N}{\sum}}\,\frac{{{c}_{in}}}{m\omega _{n}^{2}}{\int}_{{{t}_{0}}}^{t}\text{d}s\cos \left({{\omega}_{n}}(t-s)\right){{\overset{\centerdot}{{x}}\,}_{i}}(s) \end{array}\end{eqnarray} \tag{ 40 }$$

and carefully substitute q_n(t) with the use of (40) into the equation of motion (37). Assigning ${{F}_{0}}\left({{x}_{i}}-{{x}_{j}}\right)=-{{\partial}_{{{x}_{i}}}}{{U}_{0}}\left({{x}_{i}}-{{x}_{j}}\right)$ we obtain the following equation:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} M{{\ddot{{x}}\,}_{i}}(t)=\underset{j}{\overset{N}{\sum}}\,{{F}_{0}}\left({{x}_{i}}-{{x}_{j}}\right)+\underset{j}{\overset{N}{\sum}}\,{{F}_{\text{eff}}}\left({{x}_{i}},{{x}_{j}}\right) \\ \qquad \quad \;-\underset{j}{\overset{N}{\sum}}\,{\int}_{{{t}_{0}}}^{t}\text{d}s{{\mathcal{K}}_{ij}}(t-s){{\overset{\centerdot}{{x}}\,}_{j}}(s)+{{\xi}_{i}}(t) \end{array}\end{eqnarray} \tag{ 41 }$$

where:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{F}_{\text{eff}}}\left({{x}_{i}},{{x}_{j}}\right)=-\underset{n}{\overset{{\tilde{N}}}{\sum}}\,\left(\frac{2c_{in}^{2}}{Nmw_{in}^{2}}{{x}_{i}}-\frac{{{c}_{in}}{{c}_{jn}}}{m\omega _{n}^{2}}{{x}_{j}}\right) \end{array}\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\mathcal{K}}_{ij}}(t-s)=\underset{n}{\overset{{\tilde{N}}}{\sum}}\,\frac{{{c}_{in}}{{c}_{jn}}}{m\omega _{n}^{2}}\cos \left({{\omega}_{n}}(t-s)\right) \end{array}\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} {{\xi}_{i}}(t)=\underset{n}{\overset{\tilde{{N}}\,}{\sum}}\,\frac{{{c}_{in}}{{p}_{n}}\left({{t}_{0}}\right)}{m{{\omega}_{n}}}\sin \left({{\omega}_{n}}\left(t-{{t}_{0}}\right)\right) \\ \qquad \\ +\underset{n}{\overset{\tilde{{N}}\,}{\sum}}\,{{c}_{in}}\cos \left({{\omega}_{n}}\left(t-{{t}_{0}}\right)\right)\left({{q}_{n}}\left({{t}_{0}}\right)-\underset{j}{\overset{N}{\sum}}\,{{x}_{j}}\left({{t}_{0}}\right)\frac{{{c}_{jn}}}{m\omega _{n}^{2}}\right). \end{array}\end{eqnarray} \tag{ 44 }$$

The above theory differs from the single-particle generalized Langevin equation [3, 4] in a few significant aspects. First, let us look at the friction term ${\sum}_{j}{\int}_{{{t}_{0}}}^{t}\text{d}s{{\mathcal{K}}_{ij}}(t-s){{\overset{\centerdot}{{x}}\,}_{j}}(s)$. While it has the usual form of the convoluted memory kernel, it also has the contributions from every tracer. This closely resembles the response forces appearing in [10], so the dissipation of energy by one tracer is affected by the others. Another issue is the stochastic force term ${{\xi}_{i}}(t)$, which also contains the contributions from all the tracers.

Let us now restrict ourselves to the N  =  2 system and assume that ${{c}_{in}}={{c}_{jn}}$ for every pair of i and j. When we switch to the relative distance $r={{x}_{2}}-{{x}_{1}}$ and the center of mass $X=\left({{x}_{1}}+{{x}_{2}}\right)/2$ variables, the equations of motion read:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} M\ddot{{r}}\,=2{{F}_{0}}(r)+2{{F}_{\text{eff}}}\left(X-\frac{r}{2},X+\frac{r}{2}\right) \end{array}\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}\begin{array}{*{35}{l}} \frac{M}{2}\ddot{{X}}\,= & \xi (t)-2{\int}_{{{t}_{0}}}^{t}\text{d}sK(t-s)\overset{\centerdot}{{X}}\,(s). \end{array}\end{eqnarray} \tag{ 46 }$$

One can instantly notice that the dynamics of r becomes purely deterministic, i.e. equation (45) is only affected by F_eff, and not by ${{\xi}_{i}}(t)$ or the memory kernel. All the 'stochasticity' is poured into the dynamics of the mass center; this is possible due to the symmetric choice of c_in. Yet another issue is the emergence of the additional inter-particle force (42). This result might be linked to the effective interactions theory. Inasmuch as ${{\mathcal{H}}_{\text{eff}}}$ given by (5) traces out q_n and p_n from the partition function, the Mori–Zwanzig approach renormalizes the equations of motion. It is then welcome that some additional force ${{F}_{\text{eff}}}(r)$ emerges in (45). However, since we do not assume the equilibrium condition in the Mori–Zwanzig approach, this additional force might be seen as a non-equilibrium effective interaction. However, the force (42) is generally non-local, i.e. it depends on the absolute positions of particles. This effect is not present provided that the symmetry in the coupling (for both ${{c}_{in}}={{c}_{jn}}$ and ${{w}_{in}}={{w}_{jn}}$, $i\ne j$) is assumed. In this case ${{F}_{\text{eff}}}(r)=-{\sum}_{n}^{{\tilde{N}}}\frac{c_{n}^{2}}{mw_{n}^{2}}r$, which is local.

We now get back to the general N-tracer case and examine what kind of spatial correlations in the noise are predicted by the Mori–Zwanzig model. Typically, we shall assume that the initial positions and momenta satisfy the equilibrium Boltzmann distribution [4]. From (44), one can see that ${{\xi}_{i}}(t)$ is the function of ${{x}_{i}}\left({{t}_{0}}\right)$, so we will denote ${{\xi}_{i}}(t)=\xi \left({{x}_{i}}\left({{t}_{0}}\right),t\right)$ for the ith particle. We want to calculate the covariance function $\langle \xi \left({{x}_{i}}\left({{t}_{0}}\right),t\right)\xi \left({{x}_{j}}\left({{t}_{0}}\right)+r,s\right)\rangle$, where the averages are taken over every degree of freedom in t  =  t₀. From the structure of (44) we conclude that $\xi \left({{x}_{j}}\left({{t}_{0}}\right)+r,t\right)$ simply reads:

$$\begin{eqnarray}&&\xi \left({{x}_{j}}\left({{t}_{0}}\right)+r,t\right)=\xi \left({{x}_{j}}\left({{t}_{0}}\right),t\right)-r\underset{n}{\overset{{\tilde{N}}}{\sum}}\,\frac{c_{jn}^{2}}{m{{\omega}_{n}}}\cos \left({{\omega}_{n}}\left(t-{{t}_{0}}\right)\right).\end{eqnarray} \tag{ 47 }$$

Thus, the correlation function reduces to:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \langle \xi \left({{x}_{i}}\left({{t}_{0}}\right),t\right)\xi \left({{x}_{j}}\left({{t}_{0}}\right)+r,s\right)\rangle =\langle {{\xi}_{i}}(t){{\xi}_{j}}(s)\rangle -r\langle {{\xi}_{j}}(s)\rangle \underset{n}{\overset{\tilde{{N}}\,}{\sum}}\,\frac{c_{jn}^{2}}{m\omega _{n}^{2}}\cos \left({{\omega}_{n}}\left(s-{{t}_{0}}\right)\right). \end{array}\end{eqnarray} \tag{ 48 }$$

In order to quickly calculate $<\xi \left({{x}_{j}}\left({{t}_{0}}\right),t\right)>$ it is feasible to rearrange the interaction terms that enter the Boltzmann distribution:

$$\begin{eqnarray}&&\begin{array}{*{35}{l}} \underset{i}{\overset{N}{\sum}}\,\underset{n}{\overset{\tilde{{N}}\,}{\sum}}\,V\left({{x}_{i}}-{{q}_{n}}\right)=\underset{n}{\overset{\tilde{{N}}\,}{\sum}}\,\frac{m\omega _{n}^{2}}{2}{{\left({{q}_{n}}-\underset{i}{\overset{N}{\sum}}\,\frac{{{c}_{in}}}{m\omega _{n}^{2}}{{x}_{i}}\right)}^{2}} \\ \qquad \qquad \qquad \\ \;-\underset{n}{\overset{N}{\sum}}\,\frac{1}{m\omega _{n}^{2}}{{\left(\underset{i}{\overset{N}{\sum}}\,{{c}_{in}}{{x}_{i}}\right)}^{2}}+\underset{i}{\overset{N}{\sum}}\,\underset{n}{\overset{\tilde{{N}}\,}{\sum}}\,\frac{{{c}_{in}}}{mw_{in}^{2}}x_{i}^{2}. \end{array}\end{eqnarray} \tag{ 49 }$$

This first sum instantly matches the structure of $<\xi \left({{x}_{j}}\left({{t}_{0}}\right),t\right)>$ given by (44), so the integrations over ${{p}_{n}}\left({{t}_{0}}\right)$ and ${{q}_{n}}\left({{t}_{0}}\right)$ in $<\xi \left({{x}_{j}}\left({{t}_{0}}\right),t\right)>$ lead to:

$$\begin{eqnarray}&&\langle \xi \left({{x}_{j}}\left({{t}_{0}}\right),t\right)\rangle =0.\end{eqnarray} \tag{ 50 }$$

In fact, the above result is expected for the stochastic force. We can also calculate $\langle {{\xi}_{i}}(t){{\xi}_{j}}(s)\rangle$ to conclude:

$$\begin{eqnarray}&&\langle \xi \left({{x}_{i}}\left({{t}_{0}}\right),t\right)\xi \left({{x}_{j}}\left({{t}_{0}}\right)+r,s\right)\rangle =\frac{1}{\beta}\underset{n}{\overset{{\tilde{N}}}{\sum}}\,\frac{{{c}_{in}}{{c}_{jn}}}{m\omega _{n}^{2}}\cos \left({{\omega}_{n}}(t-s)\right).\end{eqnarray} \tag{ 51 }$$

This recovers the celebrated relation between the temporal correlations in the noise and friction memory kernel [3–5]. In the purely spatial context, the stochastic forces prove extremely strongly spatially correlated, i.e. for t  =  s:

$$\begin{eqnarray}&&\langle \xi \left({{x}_{i}}\left({{t}_{0}}\right),t\right)\xi \left({{x}_{j}}\left({{t}_{0}}\right)+r,t\right)\rangle =\frac{1}{\beta}\underset{n}{\overset{{\tilde{N}}}{\sum}}\,\frac{{{c}_{in}}{{c}_{jn}}}{m\omega _{n}^{2}}\end{eqnarray} \tag{ 52 }$$

which is constant in the entire space. In fact, on average, all the tracers are moved by the thermal bath as one entity.

We will now summarize the properties of the Mori–Zwanzig approach and compare them to the SCN-based dynamics of the multi-tracer diffusion. The Mori–Zwanzig approach is fully microscopically consistent and, at least formally, it describes the spatio-temporal aspects of simultaneous diffusion. Yet, the Mori–Zwanzig model is limited to the linear tracer-bath coupling and neglects the bath–bath interactions, unlike our model, which is not restricted in these aspects. Both models agree that the dissipative term for a single tracer can be decomposed into the proper friction term and the sum of the response forces, which are the contributions from the other tracers. This supports the idea of collective diffusion.

One question regarding our SCN-driven dynamics is whether or not the effective forces should be explicitly included in the Langevin equation (10). The observation retrieved from the Mori–Zwanzig model (45) is that the effective force emerges in the renormalized equation of motion, thus justifying the former option. Although the effective interactions predicted by the Mori–Zwanzig model are limited only to the relatively simple form (42), it is also remarkable that they can be both equilibrium or non-equilibrium. This is also in contrast with [10], where we are limited to the equilibrium effective interactions.

Another problem is the $K(r)\to 0$ and ${{g}_{-}}(r)\to 0$ behavior as $r\to 0$ in (10), which is the limiting case for our over-damped dynamics. Our theory predicts that the influence of both friction and noise should vanish. The Mori–Zwanzig model shows that, indeed, the dynamics of r can become purely deterministic provided that there is a symmetry in the coefficients ${{c}_{in}}={{c}_{jn}}$ for $i\ne j$ in (36). In this case, the equation of motion for r reads (45), which has no stochastic contribution.

The problem of the symmetry ${{c}_{in}}={{c}_{jn}}$ for $i\ne j$ is related to the interpretation of the coupling potential (36). If we perceive it as a quadratic expansion of some otherwise non-linear function, there is the question of whether we expand this nonlinear $V\left({{x}_{i}}-{{q}_{n}}\right)$ about the same position for every tracer (e.g. $V\left({{x}_{i}}-{{q}_{n}}\right)\simeq V(0)+{{V}^{\prime}}(0)\left({{x}_{i}}-q\right)+\frac{1}{2}{{V}^{\prime \prime}}(0){{\left({{x}_{i}}-q\right)}^{2}}$) or not. In the former case, the symmetry is naturally ensured, but not in the latter. This is also partially related to the range of applicability of the Mori–Zwanzig model. On the one hand, it can be constructed wherever the linear approximation applies to the coupling forces; this usually requires a local energy minimum. On the other hand, the Mori–Zwanzig model predicts the constant spatial correlation function of the stochastic forces, but the realistic expectation (confirmed by the microscopic relation (7)) is a function which decays with distance. Thus, one can see the Mori–Zwanzig model as a limit of strong correlations and identify it with the short-distance limit of SCN-based theory. At the same time, this latter theory seems applicable as long as the equilibrium correlation function and the Gaussian SCN are adequate approximations.

Yet another important aspect of the multi-tracer Mori–Zwanzig model is the predicted spatio-temporal correlation function. This function shows that the entire group of tracers within a range can move in a highly coordinated manner. This resembles the behavior of glasses [13–17]. It is also well known that the temporal correlations in the motion of a single tracer can result in sub-diffusive dynamics, i.e. its mean square displacement (MSD) grows as $\propto {{t}^{\alpha}}$ (where t  =  time, $0<\alpha <1$). The single-tracer generalized Langevi equations (in fact: the Mori–Zwanzig theory) are commonly used to model this effect (e.g. [5, 6]). The multi-particle Mori–Zwanzig model extrapolates these predictions on the group of tracers, provided that the constants c_in are chosen so the memory kernel (43) also takes the form of the power-law [5, 6]. On the other hand, while the sub-diffusive dynamics is observed in dense colloids, it is claimed to manifest itself via the logarithmic growth of MSD ($\propto \ln t/\tau$) [47, 48]. This imposes an intriguing problem: on the one hand, the logarithmic MSD is usually modeled as a result of the non-Gaussian noise (e.g. [49]), not the memory effect; on the other hand, the logarithmic and fractional power-law functions can mimic each other even over several orders of magnitude, when properly adjusted. Supposing that the actual experimental dependence is, in fact, of the power-law type, the Mori–Zwanzig model could be readily applied to reproduce this data. However, this would require a change in their current interpretation.

Finally, let us comment on the limitations of the Mori–Zwanzig approach. The main restriction is the linearity in coupling, but, unfortunately it is also its essential ingredient. Linearity results in the somewhat over-simplified form of the effective interaction (42) and in the constant spatial correlation function. This last result is particularly unrealistic, especially for a huge distance r between the tracers. These impediments could be partially remedied by the coupling force of the quasi-non-linear form e.g. $\left({{q}_{n}}-{{x}_{i}}\right)G\left({{x}_{i}}\right)$, for which it is also possible to obtain some analytical results [4]. However, such a coupling is not invariant under the global translations (i.e. $\left(\left({{q}_{n}}+\epsilon \right)-(x+\epsilon )\right)G\left({{x}_{i}}+\epsilon \right)\ne \left({{q}_{n}}-{{x}_{i}}\right)G\left({{x}_{i}}\right)$), so this solution is also dissatisfying. All of these problems show that the Mori–Zwanzig approach can only serve as the auxiliary short-distance model and a more versatile workaround is necessary, such as, e.g. our SCN-driven Langevin dynamics. Obviously, our theory is restricted to SCN-dominated systems as it cannot handle the temporal correlations. For this reason, the Mori–Zwanzig model suggests that a similar approach, but covering the non-linear coupling, could be the ultimate solution to the problem of spatio-temporal correlations in dynamics.

The theory of diffusion in colloids has gone a long way from its initial formulations up to the present day. Its significance has only grown as it has been involved in the understanding of biological and self-organizing systems. In this paper we have attempted to review and summarize the growing field of research on SCN-related phenomena. We have also expounded the idea of collectivity in diffusion and provided two complementary models that reflect this idea. The collective perspective expands the importance of the diffusion theory to the field of molecular self-assembly and sheds new light on the stochastic dynamics of transient, non-equilibrium states. The theory of thermodynamically consistent SCN-driven Langevin equations [10] and the multi-particle Mori–Zwanzig model provided in this article amount to a coherent picture, establishing the framework for collective diffusion. However they must be considered as intermediate steps towards an even more comprehensive theory of spatio-temporal correlations, which could possibly deal with non-equilibrium states without resorting to equilibrium or linear approximations.

The authors gratefully acknowledge the National Science Center, Poland for the grant support (2014/13/B/ST2/02014).