Dimension reduction of noisy interacting systems

We consider a class of models describing an ensemble of identical interacting agents subject to multiplicative noise. In the thermodynamic limit, these systems exhibit continuous and discontinuous phase transitions in a, generally, nonequilibrium setting. We provide a systematic dimension reduction methodology for constructing low dimensional, reduced-order dynamics based on the cumulants of the probability distribution of the infinite system. We show that the low dimensional dynamics returns the correct diagnostic properties since it produces a quantitatively accurate representation of the stationary phase diagram of the system that we compare with exact analytical results and numerical simulations. Moreover, we prove that the reduced order dynamics yields also the prognostic, i.e., time dependent properties, as it provides the correct response of the system to external perturbations. On the one hand, this validates the use of our complexity reduction methodology since it retains information not only of the invariant measure of the system but also of the transition probabilities and time dependent correlation properties of the stochastic dynamics. On the other hand, the breakdown of linear response properties is a key signature of the occurrence of a phase transition. We show that the reduced response operators capture the correct diverging resonant behaviour by quantitatively assessing the singular nature of the susceptibility of the system and the appearance of a pole for real value of frequencies. Hence, this methodology can be interpreted as a low dimensional, reduced order approach to the investigation and detection of critical phenomena in high dimensional interacting systems in settings where order parameters are not known.

The investigation of dynamical phenomena in complex networks constructed according to different topologies is an extremely active research area [1][2][3].Interacting agent based models are commonly employed to model various phenomena in the natural sciences, social sciences and engineering [4,5], such as cooperation [6], synchronisation [7], systemic risk [8] and consensus formation [9].Several algorithms for sampling, optimization and the training of neural networks can be interpreted as interacting particle systems [10][11][12].In the thermodynamic limit, such models often exhibit phase transitions as a result of the complex interplay between the interacting dynamics and the noise.Clearly, singularities associated to phase transitions, such as the divergence of correlation properties [6] and the breakdown of linear response properties [13,14], can only be observed in the mean field (thermodynamic) limit.Consequently, their investigation involves the study of (nonlinear and nonlocal) mean field Fokker-Planck equation or a brute force approach, i.e., * niccolo.zagli@su.seextensive numerical simulations of very large ensemble of agents.Reduction of complexity can be achieved by defining collective variables (reaction coordinates) able to accurately describe the full dynamics in a low dimensional space.Nonetheless, while order parameters like magnetization can in many cases easily deduced for equilibrium systems using, e.g.symmetry arguments, the definition of reaction coordinates for nonequilibrium system is far more challenging [15][16][17].
The goal of this paper is to present a model reduction approach for the study of such infinite systems based on a systematic approximation of the full infinite dimensional dynamics in terms of a low number of ODEs.This methodology can be applied to any interacting systems model with mean field polynomial dynamics, with numerous applications including cooperation phenomena [6], synchronisation of nonlinear, possibly chaotic, oscillators [18,19] and emergent phenomena in neural networks and life sciences [20,21].The dimension reduction procedure we propose is based on a suitable closure method of the infinite hierarchy of equations for the moments or, equivalently, cumulants of the probability distribution of the infinite dimensional system.Such closure method results in a deterministic parametrization of the full dynamics in terms of a low number of cumulants.One could potentially improve on this by using the Mori-Zwanzig formalism [22,23] to construct a stochastic, possibly non-Markovian, parametrization [24][25][26].From a data-driven perspective, one could rely on empirical model reduction [27] techniques to obtain closures from partial observations of the system.The resulting closure structure is given in terms of multilayer stochastic systems whose relevance and robustness has also been highlighted from an alternative, theory-informed parametrization perspective [28].As validation case studies, we apply our dimension reduction methodology to investigate the nonequilibrium continuous phase transition in a model featuring noiseinduced stabilisation phenomena [29] and a model featuring an equilibrium discontinuous transition [30].

I. THE CLASS OF MODELS
We consider a system of exchangeable weakly interacting one-dimensional diffusions whose dynamics is governed by the following Stratonovich SDE with initial condition x i ∼ ρ(x) and i = 1, . . ., N .Each agent undergoes an internal dynamics given by the vector field F α (x), depending on a set of parameters α, and is coupled with all the other agents through a symmetric interaction potential U(x) = U(−x), with θ denoting the interaction strength.Furthermore, dW i , i = 1, . . ., N , are independent Brownian motions and σ(x) > 0 ∀x ∈ R is a multiplicative diffusion coefficient.The main assumption in this paper is that F (x), U(x) and the diffusion matrix Σ(x) = σ 2 (x) all have a polynomial functional form.We consider quadratic interactions, U(x) = x 2 2 , corresponding to cooperative interactions among the agents that attempt to synchronise them towards their common centre of mass x(t) = 1 We are interested in the thermodynamic limit N → +∞ of Eq. (1).It is known that the empirical measure [31][32][33] to the one particle distribution ρ(x, t) satisfying the nonlinear, nonlocal Fokker-Planck (McKean-Vlasov) PDE that, according to our setting, can be written as where ρ(x, 0) = ρ(x) and x = R yρ(y, t)dy represents the first moment of the distribution ρ(x, t) and Fα (x) = F α (x) + 1  2 σ(x)σ (x).Eq. ( 2) exhibits, at low temperatures, non-uniqueness of stationary solutions, that correspond to phase transitions [30,34].Stationary solutions of Eq. ( 2) can be written as a one parameter family of distributions where the parameter m satisfies the selfconsistency equation and Z(m) > 0 denotes the partition function.
Eq. ( 5) plays a major role in determining the stationary properties of the system.Solutions m of Eq. ( 5) correspond to stationary measures ρ 0 (x; m ) with first moment x = m , a suitable order parameter of the system for this type of quadratic interactions.Partial information on the stability of the invariant measures can be obtained by the investigation of the slope of the selfconsistency equation R (m ) = dR(m) dm | m .In particular, if R (m ) > 1, the stationary solution ρ(x; m ) is unstable.

II. REDUCED ORDER DYNAMICS
In order to construct the reduced order dynamics, we multiply Eq. ( 2) by x n , n ∈ N, and integrate over R. Given our assumptions on the drift and diffusion terms, this procedure results in an infinite hierarchy of equations for the moments M n = x n of ρ, see Appendix B for more details.In order to elucidate the procedure above, we will first consider model A defined by 2 is a double well potential if α > 0 and the diffusion matrix is Σ(x) = σ 2 +σ 2 m x 2 .This model was introduced in [29] to investigate the effect of multiplicative noise on spatially extended systems.We mention that, if σ m = 0, model A becomes the well-known Desai-Zwanzig model [35], a paradigmatic example featuring an equilibrium continuous phase transition.The state dependent noise arises as the parameter α is not known with infinite precision and is allowed to randomly fluctuate in time, namely α → α + σ m dξ where dξ is another uncorrelated Brownian motion.Model A shows a noise induced stabilisation phenomenon.When σ m = 0, the multiplicative noise has a rectifying effect, pushing, for strong enough coupling θ, the phase transition point to higher and higher σ, see Appendix A and in particular figure 3 for more details.We apply the procedure mentioned at the beginning of this section to model A and obtain the following equations for the moments M n with M 0 = 1, M −1 ≡ 0. Firstly, we observe that the global coupling among the agents gives rise to an interaction term between the order parameter x = M 1 and all the other moments M n , introducing a nonlinear term in the hierarchy for the moments.Secondly, the nonlinear features of the dynamics given by V α (x) introduce a (linear) dependence of lower moments on higher degree ones.The infinite hierarchy of moment equations ( 6) is equivalent to Eq. ( 2) and no reduction in the level of complexity of the mathematical description has been accomplished yet.The necessity of finding appropriate closure schemes for the hierarchy arises.Were we to truncate the system of Eqs. ( 6) at a specific level n, a closure scheme for M n+1 , M n+2 in terms of M n with n < n is needed.Truncated moment problems and closure schemes are not easily amenable to a mathematical investigation and are known to introduce statistical assumptions whose validity is difficult to justify, if not from an a posteriori perspective, see [36][37][38].
Following [35,39,40], we implement a cumulant truncation scheme [40][41][42].We introduce the cumulants k n as The truncation scheme consists of imposing the condition k n+1 = k n+2 = 0.This procedure provides a closure relations for Mn+1 = Mn+1 (M 1 , . . ., M n) and Mn+2 = Mn+2 (M 1 , . . ., M n).Alternatively, one can obtain from ( 7) and (2) an infinite hierarchy of equations for the cumulants where the explicit expression of the nonlinear function G n (•) is written appendix B. Eq. ( 8) indicates that the cumulant truncation scheme corresponds to a parametrization of the dynamics given by Eq. (1), in the limit N → +∞, in terms of a finite number n of cumulants.It is well known that such a scheme is inconsistent, since a function with a finite cumulant expansion cannot be positive if the order of the highest cumulant is larger than two [43].Heuristically, a parametrization in terms of cumulants is expected to perform better than parametrizations in terms of (central) moments based on the observation that a Gaussian distribution has vanishing cumulants k n = 0 for n > 2, while all (central) moments are nonzero.For non-Gaussian distributions, one expects that neglected higher-order cumulants will be smaller than the corresponding (central) moments.Moreover, the relevance of cumulants in the description of statistical properties of complex systems, especially in settings with athermal noise, has recently been highlighted, see [44,45] and references therein.We refer the reader to appendix C for the comparison between different parametrizations and the validation of the cumulant truncation scheme for the systems under investigation.For model A, Eq. ( 5) predicts that the stable solution x = 0 bifurcates when R (0) = 1 through a continuous phase transition in two symmetric, competing states with opposite order parameter.Panel (a) of Fig. 1 shows the continuous phase diagram for the state with positive order parameter, obtained with the exact selfconsistency equation and the reduced order dynamics, see Eq. ( 8).As soon as n = 4 cumulants (main panel) are introduced, the reduced dynamics provides a very good approximation of the phase diagram.The critical value of the parameter is underestimated by the reduced order dynamics, with such approximation getting progressively better as more cumulants are considered.The accuracy of the reduced dynamics has been quantitatively assessed in terms of the absolute error ∆ (shown in the inset) with respect to the selfconsistency approach.The reduced dynamics has also been compared to numerical simulations of an ensemble of N = 12000 agents described by Eqs.(1).We have used the Milstein scheme [48], which has strong order of convergence 1, with time step ∆t = 0.01 and estimated the order parameter as the time average, at stationarity, of the center of mass x(t).Moreover, the reduced order dynamics has been initialised with a Gaussian initial condition, such that (k 1 , k 2 ) = (0.1, 0.01) and all others cumulants set to zero.Very good agreement is observed.Close to the phase transition, finite size effects arise in the numerical simulations.Noise-induced transitions among the two symmetric solutions become a relevant feature and one should consider the rectified order parameter (shown in the figure), obtained as the time average of x(t) conditioned on the fact that the system is in the basin of attraction of the positive solution.We have also probed the validity of the cumulant based parametrization by investigating discontinuous phase transitions.We introduce model B, characterised by a tilted potential V α,µ = V α (x) + µx and additive noise Σ(x) = σ 2 .Stationary properties of the reduced dynamics, with Gaussian initial condition (k 1 , k 2 ) = (1, 0.01), are in very good agreement with the other two approaches, see Panel (b).The insets show the relative error ∆ rel between the reduced dynamics and the selfconsistency equation.The top one, referring to the top branch of the phase diagram, shows that in the very close proximity, represented as a shaded area, of the transition point, ∆ rel jumps to higher values, due to the fact that the reduced dynamics underestimates the critical value of the parameter and approaches it from below as n increases..The bottom inset shows that ∆ rel for the bottom branch of the phase diagram is instead a smooth function that is not affected by the transition.This confirms that the reduced dynamics is able to track, as σ is parametrically changed, the disappearing attractor until a transition occurs to the other stable, smoothly changing, attractor.Noise-induced transitions are observed close to the phase transition in the finite system.Due to the asymmetry between the two competing states, the metastable lifetime of the state with x > 0 decreases as the transition is approached and the system, after a short time, is driven to the other state of much longer lifetime.The above results confirm that the reduced order dynamics correctly retains information of the exact invariant measure of the system.Below, we show that the approximate dynamics also captures time-dependent properties, and, specifically, correlations, by investigating, in the spirit of the fluctuation dissipation theorem [46,47], its dynamical response to time-modulated external perturbations.We report linear response properties of the reduced dynamics for model A (see appendix A for response properties of model B).We perturb a stable stationary state with a homogeneous perturbation in the drift term F (x) → F (x) + εT (t), where ε is small.Such procedure results in a one-cumulant per-turbation k 1 + εT (t) for Eqs.(8), where k (0) 1 is the unperturbed order parameter.Following [14], we choose as temporal modulation for the forcing a Dirac's δ: T (t) = δ(t), which corresponds to a broad band forcing in frequency space.We then observe the Green function G(t), associated to the order parameter, defined by Convergence to the linear regime has been assessed evaluating the response for different values of ε.Panel (a) of Fig. 2 shows that, at the transition point (red lines), the Green function has an exponential decay (see inset) with an associated timescale that is order of magnitudes greater than what is observed in non-critical settings (blue and black lines).Moreover, such timescale is an increasing function of the level of truncation n of the reduced dynamics, whereas no dependence on n is observed for the noncritical Green functions -see panel (a) of Fig. 2. The critical behaviour is linked to the breakdown of linear response theory at the phase transition point, in the thermodynamic limit of Eq. ( 1) due to the agent-to-agent interactions, thus being associated with endogenous dynamical processes [13].As the number of agents N is increased, one observes an emerging singular behaviour in the susceptibility χ(ω), defined as the Fourier Transform of G(t), signalled by a development of a pole ω 0 on the real axis of the frequencies [14].The infinite hierarchy (7) corresponds to the thermodynamic limit of the ensemble of agents and one expects a diverging response in critical settings.However, we observe that the truncation scheme introduces a mollifying effect of the singular behaviour of the reduced response operators.The resonance of such operators can be investigated through the susceptibility of the reduced dynamics that can be written as χ(ω) = κ ω−ω0+iγ(n) +r(ω) where ω 0 = 0 and r(ω) is an analytic function in the upper complex ω plane.As the number of reaction coordinates increases, n → ∞, the regularising effect vanishes, γ(n) → 0, and the susceptibility develops a singular behaviour given by lim n→∞ χ(ω) = −iπκδ(ω − ω 0 ) + κP Panel (b) confirms the appearance of an emerging pole with an imaginary residue κ = i|κ|.The real part χ RE (main panel) of the susceptibility clearly shows the resonant δ−like behaviour for ω = ω 0 .Alternatively, the top inset shows that the primitive function of χ RE close to the pole (c = −0.01)converges accordingly to a Heaviside function.We observe that n = 4 does not show a resonant behaviour, even though it is associated with a longer timescale.The imaginary part χ IM (ω) of the susceptibility (bottom inset), behaving like a Cauchy principal value distribution, yields a quantitative estimate |κ| ≈ 1 for the residue of the pole.It is possible to obtain a formula for the amplitude of the residue as |κ| = 1 θτ x,A where (see appendix D) Numerical simulations on an ensemble of N = 16000 agents yield a value of τ x,A ≈ 0.25 and, since θ = 4, |κ| ≈ 0.99, validating thus our results.We remark that the existence of the pole ω 0 at the phase transition, as opposed to its residue κ, depends neither on the forcing T (t) nor on the choice of the observable and can be related to spectral properties of suitably defined evolution operators [14].This crucial property validates the use of our cumulant based reduced dynamics to settings where the order parameter is not known or cannot easily be written in terms of the cumulants.

III. CONCLUSIONS
In this paper, we considered a class of models describing an ensemble of N identical interacting agents subject to multiplicative noise that exhibits phase transitions in the thermodynamic limit.We derived a reduced lowdimensional system for the moments of the probability distribution function of the mean field dynamics.We showed that such approximate dynamics provides an accurate representation of the stationary phase diagram, even for a very low number (e.g. 4) of moments.This indicates that the cumulants act as effective reaction coordinates, which are able to capture the essential properties of the system with moderate loss of information due to the cumulant truncation.Additionally, the linear response properties of the projected dynamics agrees with that of the full system, and the breakdown of the corresponding linear response operators can be used to characterise the phase transition occurring in the system.Hence, our methodology seems useful for performing linear stability analysis for a large class of interacting multiagent systems, and for predicting their response to forcings of general nature.It is worth investigating how our dimension-reduction methodology compares with what one would obtain by applying variational autoencoders [49] to construct a surrogate, low dimensional representation of the system.On top of the detection of critical phenomena for high dimensional systems, a further application of our methodology relates to the issue of parameter estimation for interacting systems.Current parameter estimation techniques rely on suitable fitting procedures of the observational data to the infinite dimensional dynamics [50], whereas one could envision simpler settings where the reduced order dynamics is taken as the reference point.We expect that this complex reduction methodology will not prove to be as effective when the system does not exhibit a clear separation of time or phase space scales, see [51] and references therein for a review of systems that can be "effectively reduced" either from a theoretical or data-driven perspective.In this section we provide further details on the models we have studied in the paper.As specified in the main text, we investigate multi agent systems whose dynamics is given by the following equations where i = 1, . . ., N .The examples we have provided refer to a quadratic interaction potential U(x) = x 2 2 .This results in where x(t) = 1 is the common centre of mass of the system.Given that the interaction potential is convex, phase transitions of the system arise from non convexity features of the local vector field F (x). Model Awas introduced in [29] to study the effect of multiplicative noise on spatially extended systems.We consider the Desai-Zwanzig model [6,35,52] settings where the local dynamics F (x) = −V α (x) is given by a double well potential 2 and the noise is additive σ(x) = σ.The equations for motions are given by where the Ito convention is now used.The above equations describe a system at equilibrium.In the N → ∞ limit, it is useful to introduce the free energy functional The above equation describes the energy balance in the system: V[ρ] represents the internal energy associated to the local potential V α (x), W[ρ, ρ] is the energy given by the interaction among the agents and, lastly, S[ρ] is an entropic contribution.As explained in the main text, the empirical measure converges in the N → ∞ limit to a one agent distribution ρ(x, t) satisfying a non linear and non local Fokker Planck Equation.The corresponding non linear Fokker Planck Equation of equations (A3) can be written in terms of the Free Energy as Remarkably, this equation belongs to a rich class of dissipative PDEs, including the heat equation, the porous medium equation and the diffusion-aggregation equation, that are gradient flows with respect to the Wasserstein metric on the space of probability measure with finite second moment, see [53]  If an unique minimiser of the free energy exists, the dynamics converge exponentially fast, in relative entropy, to the unique stationary state and the rate of convergence to equilibrium can be established [55].However, the minimiser is not necessarily unique and multiple stationary solutions can coexist.Furthermore, convexity properties of the free energy functional provide a one-to-one characterisation of the stability properties of the stationary solutions.
The model we have investigated in the main text arises from the assumption that the parameter α is not known exactly but rather erratically fluctuates in time, that is α → α + σ m dξ where dξ is another, uncorrelated, Brownian motion.This results in a set of equations for the N interacting agents that reads where the symbol • stands for a generic (not necessarily Ito) prescription for the equations.It is convenient to write the above set of equations in the equivalent, in law, form where σ(x) = σ 2 + σ 2 m x 2 is a state dependent stochastic term.It is well known that the presence of multiplicative noise introduce a modelling issue, since it is not clear, a priori, what prescription should be given to the stochastic integral defining the stochastic equation [56][57][58]; see also discussion in [47].We interpret Equations (A8) as a generic one parameter family of stochastic integrals parametrised by a parameter ν ∈ [0, 1].Different values of ν correspond to different prescription of the SDEs.In particular, α = 0, 1/2, 1 correspond to the Ito, Stratonovich and Klimontovich prescription respectively.Different conventions of the stochastic integral lead to different stability properties of the SDE.Remarkably, the convention for a given system might also vary depending on the operational conditions [59].In the main text of the paper we always choose ν = 1 2 .It is known that a generic SDE can be transformed into an Ito-SDE by suitably modifying the drift coefficient as [56].Since it is more convenient to work with the Ito prescription, we apply this transformation to equations (A8) and obtain 2 .The introduction of a fluctuating parameter in the drift term corresponds to applying an external, statedependent noise that breaks the detailed balance condition, thus driving the N −particle system to an out of equilibrium state.Equation (3) in the main text yields in this setting The analysis of the self consistency equation ( 5) (main text) provides insightful information on the stationary phase diagram of the model.In particular, symmetries of the problem force the system to always have the trivial solution m = 0, corresponding to disordered state ρ 0 (x; 0) of vanishing order parameter.This can be easily shown by observing that R(−m) = −R(m) since stationary distributions satisfy ρ 0 (x; m) = ρ 0 (−x; −m), see equation ( 3) and ( 4) in the main text.Moreover, if m is a solution of the self consistency equation, so is −m .We thus expect that two symmetric branches of stable solutions will arise as soon as the disordered state loses stability.The disordered state becomes unstable as soon as R (0) = 1 which reads where the expectation value • 0 is taken with respect to the stationary distribution ρ 0 (x; 0).Since the order parameter vanishes at the transition point, the above equation yields, fixed all the other parameters, the critical value σ c = σ c (α, θ, σ m ) of the strength of the additive noise.Figure 3 shows the multiplicative noise induced stabilisation phenomenon we mentioned in the main text.Indeed, the multiplicative noise has a rectifying effect, pushing, for strong enough coupling θ, the transition point to higher and higher values of σ.Moreover, the amplitude of the order parameter gets magnified, since it exceeds the maximum value √ α, the minimum point of the potential V α (x), that is attained in the low noise regime (σ → 0) when σ m = 0. Model B we have investigated features a discontinuous phase transition and is obtained by breaking the symmetry x → −x through a tilted potential as V α,k = V α + µx, with µ > 0.Moreover, the system is subject to thermal noise σ(x) = σ.The pitchfork bifurcation of invariant solutions one obtains for µ = 0 disappears.In particular, there exists a smooth, stable branch of negative order parameter x for all values of the strength of the noise σ.However, decreasing σ, a pair of solutions appear through a saddle node bifurcation, yielding another branch of stable x > 0, with the other one being unstable, see Figure 1 in the main text.The saddle node bifurcation is characterised by the condition R (m c ) = 1 that reads where m c is the value of the positive order parameter at the transition point and the expectation value is taken with respect to the stationary distribution ρ 0 (x; m c ).Since m c is not known a priori and has to be evaluated numerically by solving the self consistency equation, the above equation does not directly provide the value of the critical noise σ c at which the saddle node bifurcation takes place.Nevertheless, it provides a criterion to assess how close the critical point evaluated numerically is to the exact one by evaluating the slope R (m c ) and comparing it to the exact value 1. Model B' most interesting feature is represented by the discontinuous phase transition and the jump from the top branch to the bottom one as the parameter σ is changed.Such analysis has been performed in the main text.Nevertheless one could study the dynamical response of the system as the transition point is approached from below on the top branch.Since it is associated with the loss of stability of the invariant measure, we expect similar results to hold for this model as well.We refer to the main text and to appendix D for the explanation (and for the notations) of the linear response investigation we have performed.Figure 4 shows that the Green function associated to the order parameter x and a time delta δ(t) homogeneous perturbation develops a timescale, for settings near the phase transition, that is orders of magnitude bigger than the timescale associated to non critical settings.We remark that such behaviour does not depend on the specific form of the forcing [13].The figure refers to a level of truncation of n = 22.One could also perform an analysis by looking at different values of n.We expect to obtain similar results to what is reported in the main text.However, such analysis is more complicated here by the discontinuous feature of the transition.Firstly, the reduced dynamics transition point depends on n and the analysis becomes increasingly hard very close to the transition point, see shaded area in panel (b) of Figure 1 in the main text.Secondly, Figure 4, clearly shows that the timescale associated to the Green function is highly sensitive to small deviations, such as δ = 0.1%, from the transition point.

Appendix B: Hierarchy of equations for the moments and cumulants
In this section we provide a few more details on how to obtain the dynamical evolution of the moments and cumulants of the distribution of the infinite system ρ(x, t).As explained in the main text, ρ(x, t) satisfies a non linear and non local Fokker Planck equation that we write here in an alternative way as If we multiply (B1) by x n and integrate on the phase space R, we obtain after performing some integration by parts where • represents the expectation value with respect to the probability distribution ρ and we have introduced the moments M n = x n .We observe that the main assumption in this paper, namely the fact that we assume that the local drift F α and the diffusion coefficient σ 2 (x) have a polynomial functional form, implies that both Fα x n−1 and x n−2 σ 2 (x) can be written in a closed form in terms of the moments M n .Indeed, let us explicitly carry out these calculations for model A. Similar results hold for model B. We recall that model A is de- fined by a diffusion coefficient is From (B2) one then obtains an infinite hierarchy of equations for the moments as The above calculations have been obtained for a quadratic interaction potential U(x) = x 2 2 , but we remark that infinite hierarchies of equations for the moments such as (B4) can be obtained for any polynomial interaction potential U(x).If the functions describing the dynamics are generic, as opposed to polynomials, it is not possible to find close equations for the moments.
However, one could potentially recur to a Taylor expansion to approximate, in a controlled way, these functions as polynomials and then construct the corresponding approximate hierarchy of equations for the moments.Of course, this would introduce another source of approximation on top of the one deriving from the truncation scheme of the hierarchy.
Following [35] one can alternatively obtain an infinite hierarchy of equations for the cumulants of the probability distribution ρ.We remark that the cumulants k n are defined through the cumulant generating function G(λ, t) = ln g(λ, t) as Equation (B1) yields an evolution equation for the cumulant generating function By separating the different powers of the variable x we can write the above equation in terms of G, its derivative G (λ, t) = ∂G ∂λ and higher order derivatives as Using the definition of the cumulants given in equation (B5) and comparing same powers of λ one finally obtains the equations for the cumulants 1 n This section is divided in two parts.In the first, we provide the algebra to perform a cumulant truncation scheme at any generic order n for the hierarchy of equations for the moments (B4).Secondly, we compare the performances of multiple truncation schemes and assess that the cumulant truncation scheme correspond to the best parametrisation choice for the thermodynamic limit of the interacting agents system.

Cumulant Truncation Scheme
The relationship between cumulants and moments of a probability distribution is where B nl (M 1 , . . ., M n−l+1 ) are partial (incomplete) Bell polynomials.In particular, these polynomials are given by where the sum is taken over all the sequences j 1 j 2 . . .j n−l+1 of non negative integers such that the following two conditions hold Moreover, we will make extensive use of the following two properties of the Bell polynomials The closure approximation Mn+1 can be easily found by separating the term l = 1 from equation (C1) and using (C3), In fact, evaluating the above equation for n = n + 1 and imposing the condition k n+1 = 0 results in The evaluation of Mn+2 requires more care since it involves Mn+1 as well.Let us first observe that the cumulant k n+2 can be written as, see equation (C1), Using equation (C4) we can write where we have separated the term k = 1 and k = n + 1 from the total sum.Finally, by imposing the condition k n+2 = 0 and consistently estimating M n+1 as Mn+1 we obtain the approximated value for M n+2 as In conclusion, the cumulant truncation scheme consists in the finite set of equations ( 6) with n = 1, . . ., n along with the boundary conditions M 0 = 1 and M n+1 = Mn+1 , M n+2 = Mn+2 as given by equations (C6) and (C9) respectively.

Comparison between different truncation schemes
The infinite hierarchy of equation for the moments (B4) or cumulants (B8) are equivalent to the McKean Vlasov equation (B1) describing the thermodynamic limit of the interacting agents system.For obvious practical reasons, it is necessary to find appropriate truncation schemes to the hierarchy resulting in a finite, preferably small, number of ordinary differential equations for the moments or cumulants.In particular, common truncation schemes include a moment truncation scheme (MT), a central moment truncation scheme (cMT) and a cumulant truncation scheme (CT).These schemes correspond to imposing ad hoc boundary conditions to the hierarchy of moments or cumulants.Following [35,40] we have implemented in the main text the CT scheme and proved that the cumulants act as effective reaction coordinates for the system.The low dimensional reduced order dynamics for a small number of cumulants, resulting from the CT scheme, is able to capture both stationary and time dependent properties of the thermodynamic limit of the interacting agents system.We recall that the CT scheme of order n is equivalent to imposing the condition k n+1 = k n+2 = 0 in equations (B8).This is equivalent, as explained in the previous section, to imposing the boundary conditions (C6) and (C9) to the hierarchy of equations for the moments (B4).Instead, the MT scheme at level n is obtained by imposing the condition M n+1 = M n+2 = 0 for equations (B4).Similarly, when the above vanishing condition is applied to the central moments one obtains the cMT scheme.Figure 5 provides a quantitative comparison between the three approaches and clarifies why the CT is preferable in our settings.Panel (a) shows the phase diagram of the system.The black solid line derives from solving numerically the self consistency equation and provides a reference point for the approximate results stemming from the reduced dynamics obtained from the CT (red dots) and the MT (lines with markers) schemes.It is clear that a parametrisation in terms of cumulants provides a better approximation, fixed the order n, of the dynamics of the system than a parametrisation in terms of moments.As shown in the main text too, a parametrisation in terms of as low as n = 4 cumulants yields a good approximation of the stationary dynamics, see also the bottom left inset showing the absolute error ∆ between the CT and the self consistency equation.In particular, as explained in the main text, near the phase transition point one needs to include a higher number of reaction coordinates to achieve a better performance.On the contrary, the MT scheme yields a reduced order dynamics that does not capture the stationary properties of the system in most of the range of values spanned by the strength of the noise σ.
In order to investigate in a quantitative way the difference between the three truncation schemes we introduce the metrics δ where we have denoted with M n the central moment of order n.These metrics provide a measure, at each order of truncation n, of the difference of the magnitudes of the moments and central moments with respect to the corresponding cumulant.Panel (b) shows that δ 1 and δ 2 are positive meaning that the cumulants k n are, in magnitude, always smaller than the corresponding (central) moments, validating a posteriori our choice of using a CT scheme.In this section we provide more details about the linear response properties of model A. The ultimate goal of this section is to prove the formula for the residue of the singular part of the susceptibility χ(ω) at the phase transition.The invariant measures ρ 0 (x) of the McKean Vlasov equation, see equation (2) in the main text, satisfy the eigenvalue problem L x 0 ρ 0 (x) = 0, where the linear differential operator L x 0 is defined by where ψ(x) is a smooth function and f x 0 (x) is defined in equation (3) in the main text.We now perturb the stationary state by applying a perturbation to the drift F α (x) → F α (x) + εX(x)T (t), where ε 1.We can observe the effect of the perturbation in terms of the measure of the system as ρ(x, t) = ρ 0 (x) + ερ 1 (x, t) + . . . .Alternatively, we can investigate the time dependent properties of any observable of the system after the perturbation.In the following we will observe the response of the order parameter x and write x = x 0 +ε x 1 (t) where • 1 represents the expectation value with respect to the measure ρ 1 (x, t).We define the Fourier Transform of any function f (t) as f (ω) = f (t)e iωt dt.The response of the order parameter in frequency space is given by [13,14] x where the susceptibility χ(ω) is written as The microscopic susceptibility Γ(ω) is related to microscopic correlation properties of the system in the unperturbed state described by ρ 0 .In particular, Γ(ω) is the Fourier Transform of the microscopic response function Γ(t) that can be written as a suitable correlation function as [13] Γ where the operator L † x 0 is the adjoint of L x 0 and can be interpreted as the generator of the Koopman operator of the stationary dynamics described by ρ 0 (x).For gradient systems with thermal noise, it is possible to write Γ(t) as a time derivative of suitable correlation properties.We remark that for general non equilibrium systems this is not always possible.However, given the structure of the problem, we are able find an analogous formula for Γ(t).As described in the main text, we evaluate the response of the system to a homogeneous perturbation X(x) = 1.where we have used the definition of the adjoint of an operator and equation ( 4) in the main text to evaluate the derivative of the stationary distribution.We now define the function g(x) = − 1 σσm arctan σm σ x such that its derivative is ∂g(x) ∂x = − 1 σ 2 (x) .We then evaluate the following expression where in the last line we have introduced the correlation function between observable x and observable A = arctan σm σ x defined as As σ m → 0, the above equations are compatible with the results of [52].We have numerically estimated the correlation function C x,A (t) by evaluating the one-agent correlation function c i (t) between x i and A(x i ) and then averaging over the whole ensemble of agents (N = 16000), thus yielding C x,A (t) = 1 N N i=1 c i (t).The integrated correlation time τ x,A has been estimated by imposing a cut off T = 1.5 on the time integral corresponding to the moment after which the noisy signal takes over the exponential decay of the correlation function (see inset of Figure 6).The resulting value is τ = 0.25091 with corresponding amplitude of the residue k = 0.99636, which agrees with what has been obtained through the reduced order dynamics, see Figure 2 in the main text.

FIG. 1 .
FIG. 1. Phase diagram, x = x (σ).The continuous blue line refers to the selfconsistency equation, the red dots to the reduced order dynamics (n = 4) and the magenta dots (with errorbars) to the numerical integration.Panel (a): continuous transition given by model A. The inset at the bottom shows the absolute error ∆ between the reduced order dynamics and the selfconsistency approach.The error ∆ for n = 4 is out of scale and peaks at a value ∆ ≈ 0.1.The vertical dashed line refers to the critical condition R (0) = 1.Fixed parameters are (α, θ, σm) = (1, 4, 0.8).Panel (b): discontinuous phase transition of model B. The insets show the relative error ∆ rel between the reduced order dynamics and the selfconsistency approach.The inset at the top (bottom) refers to the upper (lower) branch of the phase diagram.The vertical dashed line is obtained numerically through the selfconsistency approach and its value has been consistently checked to yield a slope R (m) such that R (m) − 1 ≈ 10 −4 .Fixed parameters are (α, θ, µ) = (1, 4, 0.02).

FIG. 2 .
FIG. 2. Green function G(t) (panel (a)) and susceptibility χ(ω) (panel (b)) for model A. The blue (black) lines refer to a non-critical setting 5% below (above) the transition point.Red lines refer to critical settings.Fixed parameters are as in Panel (a) of Fig. 1.The red color code and the arrows correspond to increasing values of n = 4, 6, 8, 10, 14, 18, 22.In panel (b) black and blue lines have been multiplied by a scaling factor for graphical purposes.

ACKNOWLEDGMENTS
VL acknowledges the support received by the European Union's Horizon 2020 research and innovation program through the project TiPES (Grant Agreement No. 820970).The work of GP was partially funded by the EP-SRC, grant number EP/P031587/1, and by J.P. Morgan Chase & Co through a Faculty Research Award 2019 and 2021.NZ has been supported by an EPSRC studentship as part of the Centre for Doctoral Training in Mathematics of Planet Earth (grant number EP/L016613/1) and by the Wallenberg Initiative on Networks and Quantum Information (WINQ).
Appendix A: The models and references therein.The free energy F [ρ] is a Lyapunov function for the dynamics and stationary solutions of the McKean Vlasov equation are critical points of the free energy functional.In fact, the time derivative of F [ρ] along solutions of equation (

FIG. 3 .
FIG. 3. Order parameter x as a function of (σ, θ) obtained via the self consistency equation analysis.The red dashed line and the continuous red line represent the exact transition curve for σm = 0 and σm = 0, see equation (A11).The other parameters of the model are fixed and equal to α = 1, σm = 1.5, ν = 1/2.

FIG. 4 .
FIG. 4. Green Function G(t) as a function of time for model B. δ represents the relative distance from the phase transition point.

FIG. 5 .
FIG. 5. Panel (a): phase diagram for model A. The continuous black line corresponds to the phase diagram as obtained from the self consistency equation, see main text.The red dots correspond to the CT scheme of order n = 4.The continuous lines with markers correspond instead to a MT schemes of increasing order.The bottom left inset shows the absolute error ∆ between the self consistency equation and the CT scheme.Panel (b): the top (bottom) panel shows the difference in magnitude between moments (central moments) and cumulants for increasing order of truncation.Moments, central moments and cumulants have been obtained from the known expression of the invariant distribution ρ0(x; m) where m has been evaluated through the self consistency equation, see main text.Here the parameters are (α, θ, σm, ν) = (1, 4, 0.2, 0.5).Moreover, in panel (b), σ ≈ 1.

L x 0 (FIG. 6 .
FIG. 6. Correlation function Cx,A(t) as a function of time.The orange line in the inset corresponds to an exponentially decaying function y = 0.1e −t/τ where τ = 0.25.The parameters of the model are the same as in Figure 2 of the main text.

Cθ
x,A (t) = x(t)A (x (0)) 0 − x 0 A 0 = dxx exp L x 0 t A(x)ρ 0 (x) − x 0 A 0 (D8)The microscopic susceptibility can thus be written asΓ(ω) = +∞ −∞ dte iωt Γ(t) = 2 σσ m C x,A (0) + iω Ĉx,A (ω) (D9)where Ĉx,A (ω) =+∞ 0 e iωt C x,A (t) is the (one-sided) Fourier transform of the correlation function C x,A (t).We can then show that the macroscopic susceptibility χ(ω) develops a singular behaviour for a real frequency ω 0 = 0 at the phase transition.Let us observe that equation (A11), that characterises the phase transition line, can be written as θ σσ m C x,A (0x 0 = 0 at the transition point.In conclusion, using all the above results, the susceptibility χ(ω) of the system isχ(ω) Ĉx,A (ω) (D11)Being related to the spectral properties of the operator L x 0 , the quantity Ĉx,A (ω) is an analytical function at the phase transition[13,14,52]. Consequently, the above equation shows that linear response theory breaks down at the phase transition, with the susceptibility χ(ω) developing a simple pole in ω = ω 0 = 0 with residue