Long-time diffusion and energy transfer in polydisperse mixtures of particles with different temperatures

Evidence suggests that the transport rate of a passive particle at long timescales is enhanced due to interactions with the surrounding active ones in a size- and composition-dependent manner. Using a system of particles with different temperatures, we probe these effects in dilute solutions and derive long-time friction and self-diffusion coefficients as functions of volume fractions, sizes and temperatures of particles in $d=2$ and 3 dimensions. Thus, we model excluded-volume interactions for nonequilibrium systems but also extend the scope to short-range soft potentials and compare our results to Brownian-dynamics simulations. Remarkably, we show that both viscosity and energy flux display a nonlinear dependence on size. The simplicity of our formalism allows to discover various interesting scenarios that can be relevant for biological systems and active colloids.


I. INTRODUCTION
A particle in a solvent is constantly bombarded by the solvent molecules (receiving random kicks at a rate ∼ τ −1 s ) which push the particle, while the particle dissipates the excess energy through dynamical friction also exerted by the solvent at longer timescales. This leads to the well-known Brownian motion [1]. The transport is measured by the mean-square displacement (MSD) of the particle, which is determined by its diffusion constant D = T /ζ at time periods t τ s . Here, T is the temperature of the solvent, i.e., the "bath" temperature, and ζ is the friction coefficient of the particle, which depends on the viscosity of the solvent and the size of the particle. For finite concentrations of interacting particles, at much longer time spans t τ c the diffusion constant of a single tagged particle differs from the bare diffusion constant D, because the particle experiences numerous collisions with the surrounding particles at a rate ∼ τ −1 c . In equilibrium fluids, this effect is a consequence of the effective friction due to the excluded-volume interactions [2] and can only reduce the value of long-time value diffusion constant D s of the tagged particle, i.e., its self-diffusion coefficient.
By contrast, in an active fluid, a passive tracer particle can gain extra energy through interactions with the surrounding active particles [3][4][5]. The self-diffusion coefficient is then determined by the interplay between the effective friction discussed above and this energy transfer. In biological systems, crowding, energy transmission and composition are key factors on the transport of material [6][7][8][9][10], and out-of-equilibrium aspects enrich the dynamics. Notably in the cytoplasm, not only the longtime friction scales nonlinearly [11,12], but also the energy flux varies with the size of the probe particle [13]. Disentangling these elements is difficult in general polar * ilker@pks.mpg.de active models. Yet, we can attempt to establish simpler models which keep the relevant features of these systems. In particular, activity brings extra persistence with velocity v to the translational motion of active particles while the direction is randomized at a rotational timescale τ r . At intermediate timecales τ c > t τ r , this motion appears diffusive with D ∼ v 2 τ r , and can be described by an effective temperature, higher than the thermal one. Similar characteristics can be achieved by phoretic motility in artifical self-propelled particles [14]. While τ r is very short for biomolecules, it gains importance for self-propelled particles. At timescales t < τ c , the tagged particle diffusion is only perturbed by hydrodynamic interactions [15,16]. A recent study explores how the hydrodynamic field of the microswimmer leads to a ballistic motion of tracer particles at short to intermediate timescales [17]. On the other hand, at long timescales t τ c , the dynamics converges to diffusive motion where direct collisions dominate and the effect of the hydrodynamic field is decreased in the long-time self-diffusion coefficient D s [18,19]. This allows for a consistent separation of timescales to implement a multipletemperature model for studying the long-time transport properties in out-of-equilibrium systems.
In this work, we provide a theoretical and numerical analysis of the long-time behavior of transport properties and energy transfer as functions of sizes, temperatures and concentrations in mixtures of particles with different temperatures. We start from the Langevin dynamics with excluded-volume interactions while neglecting hydrodynamic interactions, and derive the long-time transport coefficients expanded up to first order in concentrations. In Section II, we describe the model and highlight the central results of this paper which we compare with the Brownian dynamics simulations. In Section III, we discuss how modeling of excluded-volume interactions in nonequilibrium systems differs from the equilibrium systems. The detailed derivations are included in Supplemental Material [41] which we hope to be useful for further studies in nonequilibrium systems as we extend the scope of microrheology approaches developed in equilibrium systems and simplify the framework.

II. MODEL AND RESULTS
We consider the motion of colloidal particles in d = 2, 3 dimensions following overdamped Langevin dynamics The total potential energy includes both the interactions between particles and an external potential which confines our system in a volume V . The position of particle i is denoted by r i , its friction coefficient by ζ i , its temperature by T i and ξ i (t) is a d−dimensional standard Gaussian white noise vector with zero mean and unit variance, i.e., ξ i (t)·ξ j (t ) = d δ(t−t )δ ij . We consider multiple species of particles. Each particle of species α is in contact with a thermostat at temperature T α and friction coefficient ζ α . In general, the observables of the multi-particle system can be expressed as a series expansion in the particle densities. At first order in concentrations, the correction to isolated particle properties results from the sum of the contributions of all two-particle clusters.
For two particles of species α and β, respectively at positions r 1 and r 2 and interacting through a pairwise potential which depends only on the distance between the particles, u αβ ≡ u αβ (|r 2 − r 1 |), we choose a coordinate system where r = r 2 − r 1 is the separation vector and R = wr 1 + (1 − w)r 2 the center of motion, with w ≡ D β /(D α + D β ) = ζ α T β /(ζ α T β + ζ β T α ). With this choice, the diffusive motions along the directions r and R are statistically independent [20]. As a result, we obtain the two-particle Langevin dynamics in the form: as ξ k (t) · ξ l (t ) = d δ(t − t )δ kl and ξ k,l (t) = 0 where D αβ r = T αβ /ζ r and ζ αβ ζαT β +ζ β Tα , and ζ αβ R = ζ β Tα+ζαT β Tα−T β . The advective term inṘ leads to energy transfer, it only exist for nonequilibrium systems, and it vanishes as ζ −1 We may now proceed to derive long-time dynamical coefficients for a tagged particle of species α that has temperature T α , friction coefficient ζ α and diameter σ α which plays a role in excluded-volume interactions. Throughout this paper, we define each secondary particle as type-β, hence {β} represents a set of different particle types surrounding the tagged particle. Accordingly, for an α-type tagged particle, there are N β pairs for αβ interactions (we take N α − 1 ≈ N α ), and we leverage the statistical equivalence of pair interactions described in Eqs. (2) and (3).

A. Long-time friction coefficient
At long times, the collisions with the surrounding particles alters the effective friction on the tagged particle. In order to determine the density dependent long-time friction coefficient, we apply a small constant test force F on a tagged particle of species α at position r 1 along the lines of Refs. [42,43]. The test force F breaks the symmetry of the steady state and creates a non-homogeneous distribution of the other particles around the tagged one. In return, this alters the force balance and induces an effective mean force F in on the tagged particle through the excluded-volume interactions with the surrounding particles. The average velocity follows a linear relation with the total applied forces u α = (F + F in ) /ζ α from which we can infer the long-time friction coefficient ζ L α as we look for u α = F/ζ L α . To account for all pair interactions with the tagged particle, we define the concentrations for each species {β} with N β particles of diameter σ β in a volume V as c β = N β /V , and the volume frac- where Ω d is the solid angle in d dimensions. We obtain [41] The summation is performed over all particle species, and hence our theory is valid for arbitrary polydisperse compositions with many particle species at linear order in concentrations. Including the size effects on solvent-based friction constants [44], this suggests a weak compositional dependence in d = 2 in accordance with experimental results [19] and a strong dependence in d = 3 as we illustrate with simulations below. We should note that the long-time friction coefficients are temperature independent and identical for equilibrium systems  for two examples; cold particle in hot bath (blue), hot particle in cold bath (orange). In (b)-(f), the color code depicts size ratios σA/σB = 0.5, 1, 2 in gray, purple and green, respectively, and dots are obtained from simulations, solid lines represent theoretical results for friction coefficient, self-diffusion coefficient, and effective temperature shift respectively from Eqs. (4),(6), (7) with scaled contact distances σ αβ → σ αβ . In (b), the long-time friction coefficient increases with crowder volume fraction φB and displays a behavior beyond Stokes' law as ζ L show how both self-diffusion (shown in log scale) and temperature vary with sizes, temperatures, and volume fraction of crowders. A slight discrepancy between the theory, Eq. (6), and simulations for self-diffusion arises for larger tracer particles for which Eq. (10) (dashed lines) show better success. In (f), we show three illustrative cases where the energy transfer balances frictional forces for temperature ratios satisfying Eq. (11).
of particles in contact with a single thermostat and for non-equilibrium systems with the same types of particles having multiple temperatures. This calculation is sufficient to determine the self-diffusion coefficient at linear order in equilibrium systems where D s α = D α ζ α /ζ L α . However, the non-equilibrium counterpart for the selfdiffusion coefficient D s α requires an investigation of the dynamics of the tagged-particle density.

B. Long-time diffusion constant
The free diffusion constant of an α-type particle is given by D α = T α /ζ α and reflects the behavior at short timescales in the absence of collisions. The long-time diffusion constant is obtained from the asymptotic behavior of the MSD (r α (t) − r α (0)) 2 of tagged particles as t → ∞ (averaged over all the α particles from many realizations). We may definē α is linked to the taggedparticle scattering function F (k, t) = e −ik·(r1(t)−r1(0)) which is the Fourier transform of the tagged-particle density auto-correlation function. Here, each incident wave is weighted with the steady-state probability distribution of finding the tagged particle at r 1 (0) at t = 0 and the conditional probability of finding the particle at r 1 (t) at time t given that it started at r 1 (0) (see Supplemental Material [41]). For small values of the Laplace variable s conjugate to t, and small wave vector k, the Laplace transform of the tagged-particle scattering function can be expanded as . To first order in concentrations, we obtain in d = 2, 3 dimensions [45]. The second term in the sum reflects the non-equilibrium contributions to the selfdiffusion constant due to the different temperatures.
To illustrate certain aspects of Eq. (6), we consider a binary system of particles with one single particle of type-A (tracer) surrounded by many type-B particles (crowders), i.e., φ A → 0 and φ B > 0 while the temperatures differ by an amount ∆T = T B − T A . When T B T A and σ A σ B , the fractional change in root-MSD scales linearly with the size of the tracer, i.e., remarkably similar to the observations on tracer diffusion in bacterial cytoplasm where the metabolic activity is tuned by energy depletion procedures [13]. For equal-sized A and B particles with particle size σ A = σ B and ζ A = ζ B = ζ, the self- For equilibrium systems, ∆T = 0, and the self-diffusion constant recovers the known value [18,46]. For non-equilibrium systems, as ∆T = 0, we investigate two opposite limits: (i) increasing the slowing down of the hot particle by the cold crowders. This effect is also observed with active Brownian particles for which the propulsion speed is more damped [47] when the surrounding bath is composed of less active particles [48].
such that the diffusion of cold particles is mainly driven by hot crowders.

C. Modified Einstein relation and effective temperatures
-We can define an effective temperature T eff α for the tagged particle by imposing an Einstein relation D s , and we find Thus, this measure corresponds only to the translational motion and the difference between the effective temperature and the real thermostat temperature reflects the energy exchange between hot and cold particles. For binary mixtures of equal-sized hard spheres, this simplifies to T eff Since the interaction potential is conservative the total energy conservation imposes that [22,24] which naturally follows from Eq. (7).

D. Brownian dynamics simulations and results in 3 dimensions
To test our predictions, we perform Brownian dynamics simulations of Eq. (1) with a soft repulsive pairwise potential, v αβ (r) = k 2 (σ αβ − r) 2 for r < σ αβ and v αβ (r) = 0 otherwise. For a short-range potential with sufficiently narrow cut-off, an effective hard-core interaction diameter can be defined in d = 2, 3, and it reads where [20,28] which recovers σ αβ = (σ α + σ β )/2 for hard spheres. This transformation allows for an equivalent hard-sphere description with scaled contact distances [49][50][51][52] using which we can compare the simulation results with our analytical formulas from Eqs.(4),(6),(7) by replacing σ αβ with σ αβ (see Supplemental Material for details [41]). In Fig.1(a), we showD α (t) (Eq. (5)) obtained from simulations and the fitting curve from which we determine the D s α values for two examples a tagged hot particle in a cold particle bath (orange) and a tagged cold particle in a hot particle bath (blue). We extract D s α by extrapolating Eq. (5) with a stretched exponential, i.e.,D α (t) = D α + (D s α − D α ) 1 − e −(t/τ ) γ and typically γ ≈ 1/2 as previously employed in hard-sphere simulations [53]. The details on the simulations and their analysis are given in Sec.V of the Supplemental Material [41]. The size dependence enters in Eqs. (4),(6), (7), both through the contact distance σ αβ and the friction coefficient which follow Stokes' law ζ α ∼ σ α . We first consider the behavior of a tagged particle A in a bath of crowders B. In all examples, we set σ B = σ, total volume V = 5 × 10 4 πσ 3 /6, κ ≡ kσ 2 /T max = 10 2 , where T max = max(T A , T B ) . In Fig. 1(b), we observe the size and concentration dependence of the long-time friction coefficient. The increase with radius is stronger than the Stokes' law as the concentration of crowders increases. Equation (4) This nonlinear scaling of the friction constant with increasing tracer size is observed in cytoplasmic transport and polymeric systems [11]. Although the interactions in those systems are not limited to steric interactions but also include hydrophobic and electrostatic interactions, one can use our formalism extended to narrow-range soft potentials, which defines an effective contact distance σ αβ . In Fig. 1(c)-(e), we show the self-diffusion constant and energy transfer respectively for T B = 10T A , T B = 100T A , T B = 0.1T A . Our theory fits particularly well the simulations for the effective temperatures, whereas we observe slight discrepancy in self-diffusion constant as the size ratio σ A /σ B increases. To propose a correction, we attempted to write using T eff α from Eq. (7), and ζ L α from Eq. (4). The outcome is shown by dashed lines and and slightly improves the agreement with the simulation results. It obviously recovers the rigorously-derived Eq. (6) at linear order, while it is an uncontrolled expansion at higher orders.
An interesting limit is that where energy transfer balances the friction from the same type of surrounding particles. This happens when then D s A = D A becomes insensitive to φ B , see Fig. 1(f). Thus, the friction of the B bath vanishes and the crowder particles appear completely transparent to the tracer particle. The motion remains diffusive (only solvent-based) at all timescales (but strictly t τ r for systems with active swimmers) as the only process is collisions with non-viscous crowders.
Next, we study binary mixtures of equal-sized particles at finite concentrations, e.g., φ A = φ B = 0.025 and T B /T A > 1. We also investigate the dependence on potential stiffness and the self-diffusion coefficients converge to the hard-sphere values with increasing κ as lim κ→∞ σ αβ = σ αβ . Upon increasing the temperature ratio, the fractional increment in diffusion coefficient of cold particles is significant and follows a linear trend. On the other hand, the fractional loss for hot particles is less sensitive to T B /T A and converges to D s T A (and κ → ∞) predicted by our theory. We also verified the energy conservation (8) for the same system ( Fig. 2(b)). This holds irrespective of the potential stiffness and the energy transfer in this example more strictly suggests T eff

III. NATURE OF EXCLUDED-VOLUME INTERACTIONS
On a final note, we stress a crucial distinction on how to model excluded-volume interactions in non-equilibrium systems. For real-world systems, the excluded-volume refers to repulsive steric interactions. In equilibrium systems, these can be treated as a boundary condition for the moving particle in which interparticle distances r < σ αβ are not allowed. This would produce the same statistics for more realistic potentials with nearly hard-sphere descriptions (e.g., Lennard-Jones potential at high temperatures) while both simplifying analytical approaches and allowing for Monte-Carlo methods. However, in such description, the diffusion of a moving particle is always constrained by the existence of the surrounding particles. Thus, even for the colder particle, we would expect no increase in diffusivity, such that D s α ≤ D α for all temperature ratios. On the other hand, this would no longer be an obligation for repulsive interactions, and the cold particle can gain energy through collisions. Conventional Monte-Carlo methods would fail to mimic this effect as jump rates are only constrained -and not enhanced-by the existence of another particle. Yet, a proper description can be recovered in lattice models by adding a simple collision rule as demonstrated in Ref. [54].
Here, our derivation conserves the analytical properties of pair distribution functions, and enables us to approach the hard-sphere limit through a smooth function (and vice versa towards soft potentials). We are therefore able to capture non-equilibrium aspects while keeping the simplicity of a hard-sphere description.

IV. CONCLUSION
In this work, we have identified the effect of energy flux and friction on the long-time self-diffusion constant of colloidal particles. We intended to study non-equilibrium mixtures, but Eq. (6) also provides an explicit formula that can be used in equilibrium systems of polydisperse solutions at T α = T β = T (more generally with σ αβ → σ αβ ). We left out the effect of hydrodynamic interactions in our current analysis, where we focused on the long-time dynamics. However, hydrodynamic interactions could be included in our framework. Yet this would be a more ambitious undertaking, because also there is no consensus even for equilibrium systems on the significance of the contribution of hydrodynamic interactions to the long-time value of the self-diffusion coefficient D s α ,compared to the value only due to direct collisions [16,18,[55][56][57]. It would be interesting to investigate these effects on tagged particle dynamics at different time windows in a study complemented by experimental results.
We hope that our theory will also shed light on experimental studies of active colloids, microswimmers and single-molecule dynamics in biological mixtures such as the cytoplasm. Notably, the theoretical predictions on the size scalings of diffusion enhancement and long-time friction at dilute concentrations of crowders show close resemblance with earlier experimental results on the cytoplasmic transport. Even though we set our theory in the mixing regime, it can be adapted to complex fluids that have spatial heterogeneities by introducing local values of the volume fractions. This would require further study and we leave this to a future work. The Langevin dynamics can be reformulated as a Fokker-Planck equation for a multi-particle probability distribution P ({r}) where {r} is a vector whose components are the positions of all the particles. The time evolution obeys ∂P/∂t = L N P where L N is the N -particle Liouville operator is the total potential energy function. As we shall see later on, the general properties of this operator are important to show the statistical equivalence of certain manipulations that we will use. On the other hand, our main mathematical formulation relies mostly on two-particle interactions. Accordingly, next we consider the two-particle system.

B. Two-particle system
The time evolution of the two-particle system P (r 1 , r 2 ) is governed by the Liouvilian operator in (S1) for N = 2.
In the main text, we have written the two-particle Langevin dynamics transformed into relative and center of motion coordinates. We follow the same analysis and distinguish particle types α and β that have temperature T α and T β , friction coefficient ζ α and ζ β , diameter σ α and σ β . We use the coordinates r = r 2 − r 1 , the separation vector and R = wr 1 + (1 − w)r 2 the center of motion with w ≡ D β /(D α + D β ) = ζ α T β /(ζ α T β + ζ β T α ). The corresponding Fokker-Planck equation reads: where the Liouvillian operators are given by: Tα−T β . The steady-state solution of the Fokker Planck equation obtained by separation of variables, P ss αβ ≡ G ss αβ (R)g ss αβ (r)/V 2 reads: and G ss αβ (R) is uniform. Note that this sets the flux along r to zero while a non-zero flux would exist along R direction reflecting the non-equilibrium nature of the problem.

II. NECESSARY IDENTITIES AND RELATIONS
We now introduce the mathematical identities, which are required to treat the pairwise hardcore interaction terms in our derivation. The steady-state pair distribution function for the two-body problem as given above reads where the second equality holds for a hard-sphere potential only for which it can be expressed in terms of a Heaviside step function where σ αβ = (σ α + σ β )/2 is the contact distance between the corresponding two particles. The spatial derivative of the distribution function is then given by a Dirac delta function: These identities will allow us to keep track of pairwise quantities and to simplify the solution, because the Dirac delta function can be treated as a boundary condition as we will demonstrate below in our calculations.

III. LONG-TIME DYNAMIC BEHAVIOR OF A TAGGED PARTICLE
Here, we detail the steps to calculate the long-time dynamical constants of a tagged particle of type α in the presence of surrounding {β} particles. As mentioned in the main text, the first-order corrections in particle densities emerge from summing over all two-particle clusters. Thus, once we solve the two-body problem, we can simply sum over all clusters that the tagged particle forms with the surrounding {β} particles.

A. Long-time viscosity
At long times, the long-time friction constant on the tagged particle is a function of the volume fractions of the surrounding particles, i.e., ζ L α ≡ ζ L α ({φ β }). In order to determine the density-dependent friction constant, we apply a small constant test force F on a single particle of type α at position r 1 (the tagged particle) along the lines of Refs. [42,43]. In the presence of the external test force, the modified Fokker-Planck equation for the pair distribution function is The steady-state solution of this equation can be obtained by introducing a perturbation to the solution of Eq.(S5). This suggests that up to linear order in F, we have g ss αβ (r, F) = e −u αβ (r)/T αβ 1 + where q(r) is a function of the distance r only. As a result, the equation for the steady-state pair distribution function reduces to: Then, in d = 2, 3, the radial-dependent part q(r) satisfies the relation where we used Eq. (S7). We first look for a solution in the range where r > σ αβ and then use the part that is coupled with the Dirac delta function as a boundary condition at r = σ αβ such that we satisfy Eq.(S11) everywhere. Thus, we look for a solution of the following equation: with boundary conditions The boundary condition (S13b) is necessary to obtain a valid solution in a large volume V . The solution of Eq. (S12) is of the form q(r) = A 1 r 1−d + A 2 r, and the condition q(∞) = 0 implies that A 2 = 0. By applying the other boundary condition (S13a), we obtain: which can then be plugged in to obtain the full form of the modified steady-state pair distribution function Eq.(S9). For the N -particle system with N = β N β , the modified steady-state pair distribution function (Eq.(S9)) induces a relaxation force on the tagged particle exerted by N − 1 particles. By summing over N − 1 pairs, and accounting for different particle species, we obtain: which is a summation over all secondary particle types (including α) wherer = r/r, c β = N β /V are the particle concentrations [58]. The integral is in d dimensions with a volume element dr = r d−1 drdΩ d where Ω d is the solid angle in d dimensions. Next, we use spherical symmetry rdΩ d = 0 and (F ·r)rdΩ d = FΩ d /d and plug in the solution (S9) along with Eq.(S14). As a result, we get where the volume fraction is φ β = Ω d 2 d d σ β d c β . The long-time friction constant is defined by the relation between average velocity and force, u α = (F + F in ) /ζ α = F/ζ L α , leading to which is Eq.(4) in the main text.

B. Self-diffusion coefficient
As discussed in the main text, the self-diffusion coefficient can be obtained from the tagged-particle scattering function. Using the time-evolution operator, the tagged-particle scattering function can be written as: where the brackets ss denotes right averaging, i.e., f ss = f P ss N dX, dX ≡ dr 1 dr 2 ...dr N , P ss N is the multi-particle steady-state distribution function, and any operator inside f acts also on P ss N (see Appendix A) . For small s and k the Laplace transform of the tagged-particle scattering function becomes We next use F (k, s) = e −ik·r1 (s − L N ) −1 e ik·r1 ss , and isolate the self-diffusion coefficient: By rearranging the terms, we have the following relation [18]: The first term in Eq. (S21) recovers the bare diffusion constant in the absence of interactions between the particles.
Applying the operator L N on the right of the resolvent operator (s − L N ) −1 and integrating by parts from the left side we obtain, in the limit k → 0, In a system of N particles, the effect of the interactions must be summed over all the particles. At lowest order in particle densities, we may consider only the contribution from all two-body clusters. Then the two-particle ensemble average of a function is f ss = 1 V 2 f g ss αβ (r)G ss αβ (R)drdR and we consider all αβ pairs identical. Using r 1 = R − (1 − w)r for the tagged particle and its differential form, we obtain: where x = x/x is a unit vector parallel to x, and we performed the summation over N −1 pairs composed of a tagged particle of type α and the surrounding {β} particles. The pairwise Liouvillian operator is given by L αβ = L r +L R using definitions from Eqs.(S2)-(S4). Since G ss αβ is uniform, we may simply integrate over R, and write f ss = 1 V f g ss αβ (r)dr. Then, by using g αβ (r) ∂u ∂rr ·k = −T αβk · ∇ r g αβ (r), we get where we replaced the Liouvillian operator as L αβ → L r since there is no longer any R dependence. We also defined a function χ(r, t = 0) ≡k · ∇ r g(r). We can then use the definition of Laplace transforms to write (s − L r ) −1 χ(r, t = 0) = e Lrt χ(r, t = 0)e −st dt and identify the time-evolution operator as e Lrt to write χ(r, t) = e Lrt χ(r, t = 0). Now, we take s → 0 and get: lim s→0 (s − L r ) −1 χ(r, t = 0) =χ(r, s = 0) (S27) whereχ(r, s = 0) is the Laplace transform of χ(r, t) at s = 0. As L r is independent of time, the time evolution imposes that the function χ is a solution of ∂χ(r, t)/∂t = L r χ(r, t). Taking the Laplace transform of this equation gives: L rχ (r, s) = sχ(r, s) −k · ∇ r g αβ (r).
This equation is equivalent to (S10) with q(r) → X(r)/ζ α and F →k. Thus, adapting the solution (S14) for X(r), we obtain:χ p(X t |X 0 ) = δ(X − X t )e L N t δ(X − X 0 )dX, where L N is the N -particle operator, X ≡ {r 1 , r 2 , ..., r N } and δ(X − X t ) = N i=1 δ(r i − r i (t)). We can then perform the integral in F (k, t) over X t to obtain: F (k, t) = P ss (X 0 )e −ik·(r1−r1(0)) e Lt δ(X − X 0 )dXdX 0 . (S3) Now, we may introduce the adjoint form in order to carry the integral over dX 0 , and hence we first rewrite: where L † is the adjoint Liouvillian operator. Then, we perform the integral over X 0 using δ(X − X 0 ) = N i=1 δ(r i − r i (0)) and obtain: F (k, t) = P ss (X)e ik·r1 e L † N t e −ik·r1 dX.