Universality class of the motility-induced critical point in large scale off-lattice simulations of active particles

Abstract

We perform large-scale computer simulations of an off-lattice two-dimensional model of active particles undergoing a motility-induced phase separation (MIPS) to investigate the system's critical behaviour close to the critical point of the MIPS curve. By sampling steady-state configurations for large system sizes and performing finite size scaling analysis we provide exhaustive evidence that the critical behaviour of this active system belongs to the Ising universality class. In addition to the scaling observables that are also typical of passive systems, we study the critical behaviour of the kinetic temperature difference between the two active phases. This quantity, which is always zero in equilibrium, displays instead a critical behavior in the active system which is well described by the same exponent of the order parameter in agreement with mean-field theory.

---

Introduction

One of the pillars of statistical physics is the concept of universality in critical phenomena. In equilibrium systems, close to a second-order phase transition, universality can be ultimately attributed to the divergence of the correlation length of the order parameter. The behavior of this growing length-scale is found to be independent of the microscopic details of the systems but is determined only by few specific features, i.e. the spatial dimensionality and the symmetries of the order parameter as firstly hypothesized by Kadanoff.¹ Depending on these parameters it is possible to trace back the critical behaviour of disparate systems within few groups, called universality classes.

One of the biggest challenges of recent years is to transfer the vast knowledge acquired on the universal behaviour of equilibrium systems into active matter physics. Active matter represents a peculiar class of non-equilibrium systems where the elementary units or agents are self-propelled objects capable of converting energy into systematic movement.^2,3 The interacting agents are often complex biological objects that exhibit self-organized behavior on large scales giving rise to many diverse living materials.⁴ In particular, self-propulsion can trigger a feedback between motility and local density causing an effective attractive interaction between active particles.⁵ This attractive force can cause a phase separation in active fluids, remarkably similar to the gas–liquid coexistence in equilibrium systems, which is called Motility-Induced Phase Separation (MIPS).⁶ Since MIPS is a very general feature of active dynamics it might play a role also in biological systems. For instance, it has been recently observed that multicellular aggregates of Myxococcus xanthus might take advantage of density-motility feedback for developing large scale collective behaviors that have been well described by MIPS.⁷

In equilibrium physics the standard gas–liquid coexistence ends with a critical point that belongs to the Ising universality class.^8,9 A natural question is whether or not the MIPS curve ends in a critical point and whether there is a region close to phase separation in which the active critical behaviour can be traced back to a specific universality class. Effective equilibrium approaches,^10,11 and field-theoretic computations¹² have previously pointed to the Ising universality class. Concerning numerical simulations, although numerous works have addressed the properties of phase-separated MIPS states,^13–20 the study of the critical region and the determination of the critical properties still remains challenging and controversial. In particular it has been shown that active Brownian particles in two-dimensions (2d) display some critical exponents deviating considerably from the Ising ones.²¹ Differently on-lattice simulations of an active model have shown results in agreement with the Ising universality class¹² and suggested that off-lattice simulations have been performed far from the scaling regime due to the small size of the systems. However, more recent simulations of the same lattice model have shown a deviation of the order parameter exponent from the Ising value.²² Since it is not a priori clear if these models lie in the same universality class, one must rely on large-scale computer simulations which are often a necessary tool to understand possible discrepancies between off-lattice and on-lattice critical exponents. This has been the case, for example, for Heisenberg fluids at equilibrium.^23,24

In this work, we report results of large-scale off-lattice simulations of Active Ornstein–Uhlenbeck Particles (AOUPs) at criticality. Performing Finite-Size-Scaling (FSS) analysis, we show that the system's critical exponents agree with the Ising universality class. We also show that the kinetic temperature difference between the two phases, emerging at criticality, displays a critical behavior which is well described by the exponent of the order parameter in agreement with mean-field theory combined with a small-τ expansion of the AOUPs model.

Model and methods

We consider a system composed of N self-propelled AOUP particles in 2d.^25,26 This model is perhaps the simplest active particle model (due to the linearity of the process producing the “active noise”) which has led to numerous novel theoretical developments.^11,27–31 It has been shown²⁷ that AOUPs exhibit MIPS for large values of the persistence time of the activity, as also displayed in Fig. 1(a–d). The equations of motion of AOUPs read

ṙ_i = μ(ψ_i + F_i) (1)

where r_i indicates the i-th particle position, ψ_i is the self-propulsion force, and is the conservative force acting on the particle. We consider two-body interactions, i.e., f_ij = −∇_riϕ(r_ij), with r_ij = |r_i − r_j| and we use a simple inverse power-law potential ϕ(r) = (r/σ)⁻¹²/12 with a cut-off at r = 2.5σ. Here σ represents the diameter of the particle and is set to 1. In eqn (2)τ is the persistence time of the active force and η is a standard white noise source, i.e. 〈η^α_i(t)〉 = 0 and , where the greek indices indicate Cartesian components. Here D is the diffusivity of the non-interacting particles and μ is the particles mobility (set to one in simulations). In the absence of deterministic forces eqn (1) and (2) yield^32,33 〈ṙ²〉 = 2v² = 2D/τ, allowing to express the diffusivity as D = v²τ, i.e. as function of v which is the (axial) free root-mean-squared velocity. In all simulations we fix v = 1 and increase τ from small to large values to observe the transition as shown by the schematic phase diagram in Fig. 1(e). Eqn (1) and (2) have been integrated numerically using the Euler scheme with a time step Δt = 10⁻³ up to N_t = 10⁹ time steps for the largest system size which is adequate to observe full relaxation of the density auto-correlation function as shown in the ESI.† In addition we start each simulation from a random initial configuration and we perform up to N_e = 10⁸ preparation steps which guarantee that all runs reach the steady-state, even close to the critical point. We perform averages over several independent runs and errors reported represent twice the standard error of the mean (see the ESI† for details on the error estimation).

Fig. 1 (a–d) Smoothed density field in a rectangular geometry for four active systems of different sizes at τ = 16.5 (where phase separation becomes appreciable). The color map encodes the local density value from yellow (high density) to blue (low density). The configurations have been shifted so that the dense and dilute phases are centered onto the four sub-boxes (panel (d)) considered in the analysis. The average density in the right sub-boxes is denoted by ρ_h (the density of the high-density phase), while the average density of the left sub-boxes by ρ_l (low-density). (e) Coexistence curve constructed with the densities of the dense and dilute phases (filled data-points), the estimated critical point is shown as an open circle.

The active particles move in a rectangular box of size L_x × L_y with a 1 : 3 ratio (L_x = 3L_y) and periodic boundary conditions. The simulated system sizes are N = (7.5, 15, 30, 60) × 10³. To avoid spurious effects due to the presence of an interface^21,34 we compute the quantities of interest only in four sub-boxes of size L = L_y/2 centered on the dense and dilute phases. These boxes are located at x = L_x/2 ± L_x/4 and at y = L_y/2 ± L_y/4. To ensure that these positions coincide with the locations of the dense and dilute phases, as in ref. 21, for each configuration we first find the center of mass of the system along x (with periodic boundaries³⁵), and we shift all particles so that the x-coordinate of the center of mass coincides with x = 3L_x/4 as shown in Fig. 1(d). We stress that the method also works at the critical point if the two phases are distinguishable enough. This is typically the case when the critical order parameter distribution shows two peaks as in many Ising³⁶ and Lennard–Jones systems.³⁷ This technique has been successfully applied to the 2d Ising model²¹ and to the lattice active model.¹²

All simulations are performed at a fixed density ρ = 0.95 by varying accordingly the box size L_x = 3L_y = (3N/ρ)^1/2. This value approximately corresponds to the critical density estimated for the smallest investigated system size, which is found to be ρ_c = 0.953(0.037) (errors from fit are given in brackets, see the ESI† for details on the estimation of ρ_c).

Results

Although finite systems cannot develop any diverging correlation length, the finite-size scaling hypothesis allows us to systematically study the critical properties away from the thermodynamic limit.³⁸ Using the finite-size scaling ansatz, we assume that a generic observable [scr O, script letter O] near the critical point behaves as , where ζ_{[scr O, script letter O]} is the critical exponent associated with the observable [scr O, script letter O], F_{[scr O, script letter O]} is a universal finite-size scaling function and ω is the power of the (subleading) correction-to-scaling exponent.³⁸ Here ν is the exponent associated with the divergence of the correlation length ξ as the control parameter is varied across the transition. In our active particle system the relaxation time of the noise τ is the control parameter, therefore we assume ξ ∼ (τ − τ_c)^−ν. Using this and ignoring sub-leading corrections we get [scr O, script letter O] = L^{ζ_{[scr O, script letter O]}/ν}G_{[scr O, script letter O]}(L^1/ν(τ − τ_c)) (where G_{[scr O, script letter O]} is a universal scaling function). This implies that, if the correct τ_c, ν and ζ_{[scr O, script letter O]} are known, all values of [scr O, script letter O] measured for different sizes should collapse onto each other when L^{−ζ_{[scr O, script letter O]}/ν}[scr O, script letter O] is plotted as a function of L^1/ν(τ − τ_c).

A particularly interesting observable is the fourth order cumulant of density fluctuations 〈Δρ²〉²/〈Δρ⁴〉 (the Binder parameter^39–41), where brackets indicate averages over configurations and over sub-boxes. The density fluctuations are computed in the four L × L sub-boxes described above, specifically 〈Δρ²〉 = 〈(N_b/L² − 〈N_b/L²〉)²〉 where N_b is the number of particles found in one single sub-box. For the Binder parameter we expect ζ_{[scr O, script letter O]} = 0 and thus it should be size-independent at τ = τ_c.

Exploiting this property we locate τ_c = 16.361(0.058) and [scr B, script letter B] = [〈Δρ²〉²/〈Δρ⁴〉]_τ=τc = 0.781(0.017) as the intersection of the data for N = 15 × 10³ and N = 60 × 10³ (Fig. 2(a)). We chose to use these two sizes because N = 60 × 10³ is the largest simulated size and N = 15 × 10³ is two times smaller in linear size. The estimated value of [scr B, script letter B] = 0.781(0.017) is lower than the corresponding value found in the triangular lattice gas ([scr B, script letter B] = 0.8321(0.0023), see the ESI† for discussion on the triangular lattice), but it is close to that found in the active lattice model¹² ([scr B, script letter B] ≈ 0.75). Note that τ_c = 16.36 is approximately the value at which cumulants of all sizes cross as shown in the inset of Fig. 2(a) where we report a magnification of the main panel in a small τ-interval around τ_c (see the ESI† for a systematic study of the crossing points). Fig. 2(b) shows a good data collapse of the cumulant data-points with the Ising exponent ν = 1. A direct way²¹ to determine ν is to consider the size dependence of the slope of the cumulants at τ = τ_c, this method yields ν = 1.03(0.10) as shown in the ESI.†

Fig. 2 FSS analysis. (a) Binder parameter for different system sizes. The intersection of the curves allows locating τ_c (vertical line). The inset shows a magnification of a small region close to τ_c (i.e. |τ − τ_c|/τ_c < 0.1), where the straight black lines indicate the crossing point of N = 15 × 10³ and N = 60 × 10³ determining the critical point; (b) collapse of data in (a) as a function of the scaling variable (τ − τ_c)L^1/ν with ν = 1; (c) susceptibility χ as a function of τ for different system sizes; (d) collapse of data in (c) onto a universal scaling function with the exponents γ = 7/4 and ν = 1; (e) order parameter for different system sizes; (f) collapse of data in (e) with the exponents β = 1/8 and ν = 1. The color code is the same for all panels (see legend in (a)).

Next we test the scaling of the susceptibility χ = 〈(N_b − 〈N_b〉)²〉/〈N_b〉 (shown in Fig. 2(c)). Fig. 2(d) shows that the scaling is very good using the Ising critical exponent γ = 7/4. In the ESI† we also show that if we fit directly the size-dependent values of χ at τ_c, we obtain γ = 1.84(0.20) which is compatible with the Ising γ. Note also that the χ in Fig. 2(c) does not show the typical peak as in the Ising model. This is due to the fact that the χ is obtained here by averaging the values of N_b in both the dense and dilute phases. In the ESI† we show that the χ (computed in the same way) for the lattice gas displays a similar s-shaped curve as a function of the inverse temperature and scales with γ = 7/4. Furthermore, in Fig. 2(e) we consider the density difference between the boxes centered in the high and low-density phases (ρ_h − ρ_l) (see Fig. 1(d)), which corresponds to the order parameter of the system. Note that this quantity is on average greater than zero since, even in the homogeneous phase, since the density around the center of mass for a given configuration is expected to be larger than the density in the low-density sub-boxes. We find that this quantity displays a good scaling with the Ising exponent β = 1/8 = 0.125 (Fig. 2(f)). It is worth stressing that in the thermodynamic limit, the order parameter would be different from zero only for τ > τ_c. It is, however, expected that for finite systems, a smooth variation of the order parameter could also be found below τ_c and that this should scale with the appropriate exponent. We have also checked that (ρ_h − ρ_l), computed as in the active system, scales with β = 0.125 in the case of the 2d equilibrium lattice gas for temperatures above the critical temperature (see the ESI† for details). A direct fit of the size-dependent critical (ρ_h − ρ_l) gives β = 0.113(0.055), which is compatible with the Ising β. To improve the accuracy of the β estimate we apply a finer technique finding the value of the exponent which optimizes the data collapse (see the ESI† for details on β estimation). This yields β = 0.133(0.022).

Next we test if our results are consistent with the Ising exponent η = 1/4 which controls the decay of the static structure factor S(k) near the critical point, i.e. S(k) ∼ k^−2+η. To this aim we compute where ρ_k is the Fourier transform of the density fluctuations and the normalization factor A is chosen so that S(0) = 〈ΔN_b²〉/〈N_b〉. To avoid the interface and focus only on the bulk phase we compute the S(k) for particles in the dilute phase considering only those particles having |x − L_x/4| < L/2, i.e. all particles in the left sub-boxes in Fig. 1(d). As shown in Fig. 3 the resulting S(k) close to criticality becomes fairly linear in the double log scale at low k. The data at low k are well fitted by the power law k^−2+η with η = 1/4 (full line) which is appreciably different from the mean-field decay S(k) ∼ k⁻² (dashed line). A direct power-law fit of these points gives 2 − η = 1.709(0.090) and η = 0.290(0.090) which are compatible with the Ising values 2 − η = 1.75 and η = 0.25.

Fig. 3 Static structure factor S(k) computed in the dilute phase for the largest system (N = 60 × 10³). Different colors refer to different values of τ (see the legend). Approaching the critical point the structure factor is well fitted by a power law S(k) ∼ k^−2+η at low k with η = 1/4 (full line). The best fit with S(k) ∼ k⁻² (mean-field) is also shown as a dashed line for comparison.

Up to this point, we have discussed quantities which display a critical behavior also in equilibrium fluids. We now show an observable that is zero in equilibrium while it exhibits a singular behavior in the active case. Since in active systems the instantaneous velocities are coupled to positions,²⁷ whenever MIPS occurs, dense regions of slow particles coexist with dilute regions of fast ones.²⁰ We then consider the average squared speed of particles in the dense and dilute sub-boxes that we indicate, respectively, with 〈|ṙ|²〉_h and 〈|ṙ|²〉_l. The quantity can be seen as the (effective) kinetic temperature difference between the two phases and its behavior is shown in Fig. 4(a). As expected ΔT̃ decreases, as the two phases progressively mix upon lowering τ. By combining mean-field theory with a small-τ approximation of the AOUP model we have derived the scaling ΔT̃ ∼ (τ − τ_c)^κ (see Appendix A). Moreover these approximations lead to the identification of the exponent κ of ΔT̃ with the exponent β of the order parameter. Indeed we find a good data collapse for ΔT̃ if we use the exponents κ = β = 1/8 and ν = 1 (Fig. 4(b)). A direct estimate of the exponent gives κ = 0.122(0.022) satisfying κ = β within the errors (see the ESI† for details).

Fig. 4 (a) Kinetic temperature difference between the dilute and dense phases plotted as a function of τ (different colors indicate different system sizes, legends are the same as Fig. 2(a)). (b) Data of (a) scaled with the exponents ν = 1 and κ = β = 1/8.

Conclusions

In this article, we have studied the critical properties of an active system undergoing MIPS in two spatial dimensions. Performing large-scale numerical simulations on GPU we have demonstrated that the critical behavior of the system agrees well with the Ising universality class. The importance of simulating large system sizes is worth stressing since previous studies on ABPs have reported different values of the critical exponents.²¹ Although it has been speculated¹² that the limited sizes employed did not allow the observation of the scaling regime, it is true that similar sizes have been exploited for the study of critical passive attractive liquids, providing numerical results compatible with the Ising universality class.³⁴ We instead suspect that for those sizes another correlation length, different from the critical one, may interfere with the scaling behaviour of the active system. A recent work⁴² has shown that the dense phase formed by active particles undergoing MIPS is made of a mosaic of hexatic micro-domains. We find that, at the critical point, the hexatic correlation length is comparable with the size of the sub-boxes employed for the FSS analysis when the system is small (N = 3750), justifying the choice of larger system sizes (see the ESI† for discussion). Indeed for a size as small as N = 3750 we find that the crossing point of the Binder cumulant is observed at quite low values, although a reasonable scaling is found. Although we expect that the microscopic dynamical details of the model (ABPs or AOUPs) should not affect the universality class of the system, large-scale simulations of critical ABPs are necessary to further clarify this issue.

Our large-scale simulation results are also consistent with recent works taking into account non-integrable active terms in a field-theoretical framework.⁴³ These models, in some parameter range, may show micro-phase separation instead of a full MIPS and must be included in a different (non-Ising) universality class. Differently when parameters allow for a full MIPS the extra terms in these active field-theories are irrelevant (in a renormalization group sense) and the system belongs to the Ising universality class. On the other hand, far from criticality, these non-equilibrium contributions could produce significant differences with respect to an equilibrium gas–liquid phase separation.^20,44–47 Within this context it would be interesting to understand how one could make the active critical point unstable⁴³ by altering the microscopic interactions and/or the dynamics.

Conflicts of interest

There are no conflicts of interest to declare.

Appendix A

Derivation of the exponent identity κ = β

To derive the exponent identity κ = β we considered an AOU particle in d = 1 and subjected it to an external potential Φ(x). It is known^33,48 that for small τ the velocity distribution of the particle is a zero-centered Gaussian with varianceto first order in τ. Here v² = D/τ is the free particle mean square velocity (which is assumed to remain constant as in our simulations) and is the potential curvature. We now assume that the total potential curvature Φ experienced by the probe particle in x is generated by the interactions with other particles: , where ϕ′′(x − x_i) is the second derivative of the pair interaction potential. This can be rewritten as where the integral extends over all space, and we have introduced the density field . By ignoring density fluctuations (mean-field approximation) we set [small rho, Greek, circumflex](x) = ρ = constant and we obtain , where is the mean potential curvature which is assumed to be positive. By using this in eqn (1) and neglecting higher order corrections we get 〈ẋ²〉 = v²(1 − τ[small phi, Greek, macron]₂ρ). We now consider the difference between the averaged squared speed in the low and high density phases, i.e.. If we now assume that (ρ_h − ρ_l) ∼ (τ − τ_c)^β near the critical point we have:

ΔT̃ ∼ (τ − τ_c)^κ (A.2)

with κ = β, which is the relation verified by the simulation data within errors. The derivation of this identity can be easily generalized to higher dimensions leading to the same result.

Acknowledgements

MP acknowledges financial support from the H2020 program and from the Secretary of Universities and Research of the Government of Catalonia through Beatriu de Pinós program Grant No. 2018 BP 00088.

References

1. L. Kadanoff, Critical behavior, universality and scaling in critical phenomena, 1971.
2. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
3. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
4. A. Klopper, Nat. Phys., 2018, 14, 645.
5. J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2008, 100, 218103.
6. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
7. G. Liu, A. Patch, F. Bahar, D. Yllanes, R. D. Welch, M. C. Marchetti, S. Thutupalli and J. W. Shaevitz, Phys. Rev. Lett., 2019, 122, 248102.
8. H. B. Callen, Thermodynamics and an Introduction to Thermostatistics, 1998.
9. C. Domb, Phase transitions and critical phenomena, Elsevier, 2000.
10. M. Paoluzzi, C. Maggi and A. Crisanti, Phys. Rev. Res., 2020, 2, 023207.
11. M. Paoluzzi, C. Maggi, U. Marini Bettolo Marconi and N. Gnan, Phys. Rev. E, 2016, 94, 052602.
12. B. Partridge and C. F. Lee, Phys. Rev. Lett., 2019, 123, 068002.
13. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489–1499.
14. A. Patch, D. Yllanes and M. C. Marchetti, Phys. Rev. E, 2017, 95, 012601.
15. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
16. D. Levis, J. Codina and I. Pagonabarraga, Soft Matter, 2017, 13, 8113–8119.
17. A. Patch, D. M. Sussman, D. Yllanes and M. C. Marchetti, Soft Matter, 2018, 14, 7435–7445.
18. S. Hermann, P. Krinninger, D. de las Heras and M. Schmidt, Phys. Rev. E, 2019, 100, 052604.
19. P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella and I. Pagonabarraga, Phys. Rev. Lett., 2018, 121, 098003.
20. S. Mandal, B. Liebchen and H. Löwen, Phys. Rev. Lett., 2019, 123, 228001.
21. J. T. Siebert, F. Dittrich, F. Schmid, K. Binder, T. Speck and P. Virnau, Phys. Rev. E, 2018, 98, 030601.
22. F. Dittrich, T. Speck and P. Virnau, 2020, arXiv preprint arXiv:2010.08387.
23. M. Nijmeijer and J. Weis, Phys. Rev. Lett., 1995, 75, 2887.
24. I. Mryglod, I. Omelyan and R. Folk, Phys. Rev. Lett., 2001, 86, 3156.
25. C. Maggi, U. M. B. Marconi, N. Gnan and R. Di Leonardo, Sci. Rep., 2015, 5, 10742.
26. G. Szamel, E. Flenner and L. Berthier, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 062304.
27. E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco and F. van Wijland, Phys. Rev. Lett., 2016, 117, 038103.
28. S. Dal Cengio, D. Levis and I. Pagonabarraga, Phys. Rev. Lett., 2019, 123, 238003.
29. L. L. Bonilla, Phys. Rev. E, 2019, 100, 022601.
30. R. Wittmann, U. M. B. Marconi, C. Maggi and J. M. Brader, J. Stat. Mech.: Theory Exp., 2017, 2017, 113208.
31. U. M. B. Marconi, M. Paoluzzi and C. Maggi, Mol. Phys., 2016, 114, 2400–2410.
32. G. Szamel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 012111.
33. U. M. B. Marconi, N. Gnan, M. Paoluzzi, C. Maggi and R. Di Leonardo, Sci. Rep., 2016, 6, 23297.
34. M. Rovere, P. Nielaba and K. Binder, Z. Phys. B: Condens. Matter, 1993, 90, 215–228.
35. L. Bai and D. Breen, J. Graph. Tools, 2008, 13, 53–60.
36. J. A. Plascak and P. Martins, Comput. Phys. Commun., 2013, 184, 259–269.
37. J. J. Potoff and A. Z. Panagiotopoulos, J. Chem. Phys., 1998, 109, 10914–10920.
38. D. J. Amit and V. Martin-Mayor, Field Theory, the Renormalization Group, and Critical Phenomena: Graphs to Computers Third Edition, World Scientific Publishing Company, 2005.
39. K. Binder, Z. Phys. B: Condens. Matter, 1981, 43, 119–140.
40. M. Rovere, D. Hermann and K. Binder, Europhys. Lett., 1988, 6, 585.
41. M. Rovere, D. W. Heermann and K. Binder, J. Phys.: Condens. Matter, 1990, 2, 7009.
42. C. B. Caporusso, P. Digregorio, D. Levis, L. F. Cugliandolo and G. Gonnella, 2020, arXiv preprint arXiv:2005.06893.
43. F. Caballero, C. Nardini and M. E. Cates, J. Stat. Mech.: Theory Exp., 2018, 2018, 123208.
44. R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo and M. E. Cates, Nat. Commun., 2014, 5, 4351.
45. C. Nardini, É. Fodor, E. Tjhung, F. Van Wijland, J. Tailleur and M. E. Cates, Phys. Rev. X, 2017, 7, 021007.
46. R. Singh and M. E. Cates, Phys. Rev. Lett., 2019, 123, 148005.
47. E. Tjhung, C. Nardini and M. E. Cates, Phys. Rev. X, 2018, 8, 031080.
48. D. Martin, J. O'Byrne, M. E. Cates, É. Fodor, C. Nardini, J. Tailleur and F. van Wijland, 2020, arXiv preprint arXiv:2008.12972.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d0sm02162h

---