In this paper, we study phase transitions for weakly interacting multiagent systems. By investigating the linear response of a system composed of a finite number of agents, we are able to probe the emergence in the thermodynamic limit of a singular behavior of the susceptibility. We find clear evidence of the loss of analyticity due to a pole crossing the real axis of frequencies. Such behavior has a degree of universality, as it does not depend on either the applied forcing or on the considered observable. We present results relevant for both equilibrium and nonequilibrium phase transitions by studying the Desai–Zwanzig and Bonilla–Casado–Morillo models.

Multiagent models feature in a very vast range of applications in natural sciences, social sciences, and engineering. We study here the Desai–Zwanzig (DZ) and Bonilla–Casado–Morillo (BCM) models, which are paradigmatic for equilibrium and nonequilibrium conditions, respectively. Phase transitions result from the coordination between the individual agents and are associated with the divergence of the linear response of the system. The occurrence of phase transitions is universal: it does not depend on the acting forcing and can be detected by looking at virtually any observable of the system. We showcase here how response theory is capable of providing a useful angle for understanding the universal properties of phase transitions in complex systems.

Agent based models are regularly employed to model various phenomena in the natural sciences, social sciences, and engineering.1,2 Multiagent systems are used to model diverse phenomena such as cooperation,3 synchronization,4 systemic risk,5 and consensus formation.6,7 They are fundamental in developing algorithms for sampling and optimization8 and they have also been used for the management of natural hazard9 and climate change impact.10 

Multiagent systems can often exhibit abrupt changes in their behavior, often corresponding to critical transitions that occur when a parameter, e.g., interaction strength or temperature, passes a certain threshold. Such transitions are often associated to cataclysmic events such as climate change, market crashes, etc.11,12 The importance of developing tools for predicting critical transitions has long been recognized. One of the main tools used in order to develop early warning signals for critical transitions is that of linear response theory.

Following the seminal contribution by Kubo,13 linear response theory represents a very powerful framework for studying the properties of statistical mechanical systems by investigating how they respond to external perturbations.14–16 Linear response theory has been successfully applied to classic problems of solid-state physics and optics17 as well as plasma physics and stellar dynamics (Chap. 5 of Ref. 18); see some examples of the application of the theory in both equilibrium and nonequilibrium systems.19–24 Rigorous mathematical foundations for linear response theory have been provided for the case of Axiom A systems25,26 (see, e.g., Ref. 27 for further developments in the context of deterministic systems) and for diffusion processes, both in finite and in infinite dimensions28,29 (see also the interesting contribution30 that bridges the deterministic and the stochastic viewpoints).

Critical transitions arise when the spectral gap of the transfer operator of the unperturbed system shrinks to zero31–34 as the Ruelle–Pollicott poles35,36 touch the real axis. Near criticality, the negative feedbacks of the system become increasingly ineffective, resulting in arbitrarily large, usually non-Gaussian, fluctuations and a divergence of correlation properties of the system.3,37,38

In the thermodynamic limit, the multiagent system can also undergo a qualitative change of its properties through a different mathematical mechanism, namely, phase transitions,3,37 defined as the exchange of stability of nonunique stationary distributions as the parameters of the systems vary (see a detailed analysis in Ref. 39).

In a previous paper,40 we derived linear response formulas for a system of weakly interacting diffusions described by an N-particle Fokker–Planck equation and have explicitly identified two qualitatively different scenarios for the breakdown of the linear response, associated with the previously mentioned critical transitions and phase transitions. We focus here on the latter case.

Phase transitions are a genuine thermodynamic phenomenon, where the divergence of the response stems from the coordination taking place, in suitable conditions, because of the coupling between the infinite number of agents composing the total system. The coupling among the subsystems results in a memory effect that leads to obtaining the macroscopic response function of the system as a renormalized version of its microscopic counterpart,40 with formal similarities with the well-known Clausius–Mossotti relation.17,41,42 The role of memory in determining criticality due to endogenous processes has been emphasized in Refs. 12 and 43. The link between phase transitions and slow decay of correlations for interacting particle systems is well established (see, e.g., Ref. 44).

In this paper, we focus in much greater detail on the relationship between the occurrence of phase transitions and the non-analyticity of the susceptibility of the system describing the frequency-dependent response of an observable to a given perturbation in the upper complex frequency plane. The singularity manifests itself as a pole that crosses the real axis of the frequency variable. We use a formalism that mirrors spectroscopic techniques that are used for investigating the frequency dependence of the optical properties of materials.17 By studying how the real and imaginary parts of the susceptibility of the systems depend on the number of agents, we are able to predict the position of the pole and the associated residue, which describe the emergence of the singularity in the thermodynamic limit. We verify that the position of the pole depends on the considered model, but, instead, that for a given model the loss of analyticity depends neither on the choice of the observable nor on the applied perturbation and is, in this sense, a universal feature of the system. Our numerical investigations are performed on the Desai–Zwanzig (DZ)45 and the Bonilla–Casado–Morillo (BCM)46 models. The DZ model exhibits a paradigmatic example of an equilibrium order–disorder phase transition, analogous to the Ising ferromagnetic transition,3,37 while the BCM model describes an out-of-equilibrium synchronization transition of an infinite collection of coupled nonlinear oscillators. As the transition point is crossed, the order parameter (magnetization) acquires a nonvanishing constant value for the DZ model and is periodically oscillating for the BCM model.

We investigate a system composed of N exchangeable interacting M dimensional sub-systems whose dynamics is determined by the following Itô stochastic differential equation (SDE):

dxk=Fα(xk)dtθNl=1NU(xkxl)dt+σS({xk})dW,
(1)

where xkRM and k=1,,N. Fα:RMRM is a smooth vector field, possibly depending on a parameter α, and W denotes a standard P-dimensional Brownian motion; S:RMRM×P is the volatility matrix and the parameter σ>0 controls the strength of the stochastic forcing, i.e., plays the role of the temperature. We consider a fully coupled system given by the quadratic (Curie–Weiss) interaction potential U(y)=|y|22. In this case, the order parameter is known and it is given by the first moment/magnetization.

The coefficient θ modulates the intensity of the coupling, which attempts at synchronizing all systems by attracting them to the center of mass. In the thermodynamic limit N, the one-particle distribution function converges to the distribution ρ(x,t) that satisfies a nonlinear and nonlocal Fokker–Planck equation3,47–49

tρ(x,t)=[ρ(x,t)(Fα(x)+θ(x(t)x))]+σ22(Dρ(x,t)),
(2)

where D=SST. This mean field partial differential equation might support multiple coexisting stationary measures at low temperatures/large interaction strengths. In particular, in a conservative system described by a confining potential Fα(y)=Vα(y) with additive noise such that S is the identity matrix, stationary solutions of Eq. (2) correspond to local minima of Vα(y). In this case, the thermodynamic limit (2) can be written in the standard form as limN+1NlogZN=infF(ρ), where ZN denotes the partition function of the N-particle system and F(ρ) denotes the free energy of mean field system.50 A stationary state is characterized by the order parameter x0 and the associated stationary distribution ρ0(x).

We now perturb the stationary state by setting Fα(x)Fα(x)+εX(x)T(t) and we study the response of the system by expanding the distribution function as ρ(x,t)=ρ0(x)+ερ1(x,t)+O(ε2). Following a tedious calculation reported in Ref. 40, the response of the order parameter in the frequency domain is written in terms of a macroscopic (or renormalized) susceptibility Γ~i(ω) as (repeated indices are summed),

xi1(ω)=Pij1(ω)Γj(ω)T(ω)Γ~i(ω)T(ω),
(3)

where Pij(ω)=δijθYij(ω) and the susceptibilities Γi(ω) and Yij(ω) are, respectively, the Fourier transform of the microscopic response functions that can be written as correlation functions in the unperturbed state as40 

Gi(τ)=Θ(τ)(ρ0(y)X(y))ρ0(y)exp(Lx0+τ)yi0,
(4)
Yij(τ)=Θ(τ)yjlogρ0(y)exp(Lx0+τ)yi0,
(5)

where Lx0+ represents the adjoint of Lx0 and 0 is the expectation value on the unperturbed state ρ0(x). The Fokker–Planck operator Lx0 appears on the right hand side of (2) [evaluated at the stationary state ρ0(x)] and its adjoint Lx0+ can be interpreted as the generator of the Koopman operator of the stationary dynamics.40,51 As such, correlation properties of the system are related to the spectrum of Lx0+.32,52,53 The renormalization of the susceptibility derives from the coupling among the subsystems; note that Γ~i(ω) inherits the poles of both Γi(ω) and of the matrix Pij1(ω). Away from criticality, both the microscopic and macroscopic susceptibilities are analytic in the upper half of the ω complex plane.

As discussed in Ref. 40, the critical behavior of this class of multiagent systems, signified by the singular behavior of the susceptibility, originates from two distinct physical phenomena that are associated to either the poles of Γi(ω) or Pij1(ω).

The case where Γi(ω) diverges pertains to the occurrence of critical transitions. It is of interest here the case where poles appear in the real ω axis for Γ~i(ω) because the matrix Pij1(ω) becomes singular. This corresponds to phase transitions originated by the coupling and do not show a divergence of correlation properties, because Γi(ω) is, instead, analytic in the upper complex ω plane. Equivalently, the spectral gap of Lx0+ remains finite at a phase transition. However, at the transition point, the usual dispersion relations need to be modified.17,40 The conditions underpinning the breakdown of linear response theory do not depend on the perturbation field X(x) or on the choice of the observable (x, in our case) and can be related to the spectral properties of a modified transfer operator.40 

Below, we present results for the Desai–Zwanzig (DZ)45 and the Bonilla–Casado–Morillo (BCM)46 models. These are composed by N interacting agents evolving according to Eq. (1). Each individual agent of the DZ (BCM) model evolves in R (R2). A detailed description of the two models is presented in the  Appendix. We repeat our experiments for various choices of N in order to detect the emergence of singularities for the combination of the parameters corresponding to phase transitions. Here, we keep fixed the values of the internal parameter α. Both models undergo a phase transition at the transition line σ~=σ(θ;α) in the parameter space (σ,θ) (see the  Appendix for the analytical evaluation of the transition line). Since one of the two parameters is redundant, we fix the coupling intensity θ and we vary, instead, the noise strength σ.

Following Ref. 54, we perform n simulations where the initial conditions are chosen according to the unperturbed invariant measure ρ0(x) and where at time τ=0, we apply a perturbation proportional to a Dirac δ function. The average of the response for the observable x over the n simulations gives an estimate G~i(τ;N) of the renormalized response function. Details on the numerical simulations are also reported in the  Appendix. Figure 1 shows the response functions G~i(τ;N) for an additive perturbation X(x)=1 for the DZ model (left panel) and X(x)=(1,0) for the BCM model (right panel). The two response functions are qualitatively different because, by and large, the one for the DZ model describes a monotonic decay, whereby the system relaxes toward the unperturbed state, while the one for the BCM combines the decay with an oscillatory behavior taking place at the natural frequency ω~=1.

FIG. 1.

Renormalized response functions as a function of time τ. Panel (a): response G~(τ;N) for the one-dimensional order parameter of the DZ model. Panel (b): response G~x(τ;N) of the first component of the bi-dimensional order parameter for the BCM model. Black and blue lines correspond to noncritical values of the strength of the noise σ. Red lines correspond to response functions at the transition point. For each value of σ, there are five lines corresponding to different values of N, namely, N=2k×103 with k=1,,5. The arrows and the color gradient indicate the direction of increasing N.

FIG. 1.

Renormalized response functions as a function of time τ. Panel (a): response G~(τ;N) for the one-dimensional order parameter of the DZ model. Panel (b): response G~x(τ;N) of the first component of the bi-dimensional order parameter for the BCM model. Black and blue lines correspond to noncritical values of the strength of the noise σ. Red lines correspond to response functions at the transition point. For each value of σ, there are five lines corresponding to different values of N, namely, N=2k×103 with k=1,,5. The arrows and the color gradient indicate the direction of increasing N.

Close modal

In the DZ model, the response functions initially undergo a fast and substantial decay, both far from and at the phase transition, associated with a time scale of order 1. However, at the phase transition, a new, much longer, timescale appears. This timescale increases monotonically with N. The same is observed in the case of the BCM model if one considers the envelope of the response function rather than the response function itself: at the transition, the decay of the oscillations becomes slower and slower as N increases.

The origin of the new timescales resides in the appearance of simple pole at ω=ω0 in the susceptibility Γ~i(ω;N), the Fourier transform of the response function. The pole is located at ω0=0 for the DZ model and at ω0=ω~=1 for the BCM model. When considering finite values of N, the susceptibilities describing the response of (virtually) any observable to (virtually) any external perturbation have a contribution of the form κωω0+iγ(N), where γ(N)0+ as N+ and κ represents the residue of the pole, because limNκωω0+iγ(N)=iπκδ(ωω0)+κP(1/(ωω0)). The quantity κC depends on the choice of observable and of the perturbation. We remark that the asymptotic property does not depend on how fast the function γ(N) vanishes for increasing values of N. Following Ref. 38, one might conjecture that for the Desai–Zwanzig model and related models, the function γ would scale as 1/N. However, we have observed here that the behavior of γ is different. This is an issue of fundamental importance that we will explore in future work, also in the case of nonequilibrium systems. Note also that, in the case of equilibrium systems, the mean field limit N and the limit TTc, where Tc is the critical temperature, do not commute.55 

We next investigate the phase transitions by looking at the properties of the susceptibilities (see Fig. 2). When σσ~, the susceptibilities do not show any singularity or any remarkable dependence on N, thus indicating that the thermodynamic limit has been reached to a good approximation. As N increases, for both the DZ model (left panel) and the BCM model (right panel), the resonance at ω=ω0 of the real part of the susceptibility approaches the limiting Dirac function πkδ(ωω0) with coefficient k>0. This singular behavior is clear from the plot of the primitive function of the real part of the susceptibility (bottom inset) that tends to step function. For both models, κ=ik is an imaginary number. Indeed, the imaginary part of the susceptibility behaves exactly as the Cauchy principal value distribution and can be used to get easily a quantitative estimate of k. The top insets of Fig. 2 show the function (ωω0)Im{Γ~i(ω)}. As N, this function converges to k everywhere except for ω=ω0.

An explicit expression for k is known40 in the case of the DZ model:56k=x20θ0+dtx(t)x(0)0. Using the statistics of the unperturbed runs, we obtain k0.89, which agrees within 2% with the one obtained from the limiting behavior of the susceptibility, thus validating our results. In the case of the BCM model, our procedure allows one to derive a direct estimate k0.44; in this case, no expression for the residue is available in the literature, and following Refs. 46 and 40, its evaluation seems cumbersome. We here observe that, by evaluating the susceptibility for finite values of N, we are able to predict the residue of the pole at ω=ω0, which appears, instead, only in the thermodynamic limit. The residue plays the role of a latent heat of phase change in classical thermodynamics. Our results, though, allow one to deal with the case of a dynamical latent heat, which is observed for perturbations occurring in a non-vanishing frequency. As discussed earlier, the singular behavior of the susceptibility has some degree of universality. By this, we mean that while for a given model, the value of the residue is forcing- and observable-dependent, its position is a fundamental property of the model itself (see the  Appendix for additional examples).

FIG. 2.

Macroscopic susceptibilities as a function of the frequency ω. Panel (a): susceptibility Γ~(ω) for the one-dimensional order parameter of the DZ model. Panel (b): susceptibility Γ~x(ω) for the first component of the two-dimensional order parameter for the BCM model. The parameter in the lower extreme of the integral is for both cases with c=0.05. Blue and black lines in panel (b) have been multiplied by a factor of 5 for visualization purposes. Color code and plotting conventions as in Fig. 1.

FIG. 2.

Macroscopic susceptibilities as a function of the frequency ω. Panel (a): susceptibility Γ~(ω) for the one-dimensional order parameter of the DZ model. Panel (b): susceptibility Γ~x(ω) for the first component of the two-dimensional order parameter for the BCM model. The parameter in the lower extreme of the integral is for both cases with c=0.05. Blue and black lines in panel (b) have been multiplied by a factor of 5 for visualization purposes. Color code and plotting conventions as in Fig. 1.

Close modal

The study of how a large network of identical agents responds to exogenous perturbations is of the uttermost importance in different areas of science. One might be interested not only in the smooth response of the system, where its properties change ever so slightly but also in the critical, nonsmooth regime, where small perturbations can lead to large and possibly undesired changes. Multiagent systems modeled as weakly interacting Itô diffusions represent a rich class of models exhibiting such critical behavior, and for which rigorous analysis and careful numerical investigations can be carried out. Usually, these phenomena are accompanied by a large spatial (among sub-systems) and temporal restructuring of the system where correlations get highly magnified. The critical behavior due to the emergence of a phase transition is a genuine thermodynamic phenomenon arising from the complex interactions among the infinite number of agents. Nevertheless, we have shown in this paper that linear response theory provides a powerful framework for detecting and anticipating phase transitions by investigating the response of a finite particle system to external perturbations. We have been able to predict the appearance of poles in the susceptibility, which describes the frequency-dependent response of the system, as well as to obtain a correct estimate of critical thermodynamic properties, such as the residue of the poles, based on the knowledge of the response for the finite particle system in two paradigmatic models describing equilibrium and nonequilibrium phase transitions. This is an encouraging starting point for improving our ability to understand and predict transitions in more complex multiagent systems.

V.L. acknowledges support received from the European Union’s Horizon 2020 program through the project TiPES (Grant Agreement No. 820970) and from the EPSRC (Project No. EP/T018178/1). The work of G.A.P. was partially funded by the EPSRC (Grant No. EP/P031587/1) and by J.P. Morgan Chase & Co. Any views or opinions expressed herein are solely those of the authors listed, and may differ from the views and opinions expressed by J.P. Morgan Chase & Co. or its affiliates. This material is not a product of the Research Department of J.P. Morgan Securities LLC. This material does not constitute a solicitation or offer in any jurisdiction. N.Z. was supported by an EPSRC studentship as part of the Centre for Doctoral Training in Mathematics of Planet Earth (Grant No. EP/L016613/1).

The data that support the findings of this study are openly available in “Spectroscopy of Phase Transitions” at https://figshare.com/projects/Spectroscopy_of_phase_transitions/101846, Ref. 57.

Weakly interacting diffusions represent a rich class of agent based models, describing a network of interacting subsystems. The local dynamics of each subsystem is determined by a smooth vector field Fα(x). The local force is, in general, nonconservative, leading to irreversible and dissipative processes that can exhibit complex behaviors, such as deterministic chaos (see Fig. 3).

FIG. 3.

Weakly interacting diffusions. The local dynamics of each subsystem, given by Fα(x), is in general dissipative and can support a wide range of complex behaviors, including deterministic chaos, as depicted in this figure.

FIG. 3.

Weakly interacting diffusions. The local dynamics of each subsystem, given by Fα(x), is in general dissipative and can support a wide range of complex behaviors, including deterministic chaos, as depicted in this figure.

Close modal

An all-to-all coupling between the subsystems is given by a matrix Lij=U(xixj), where U(x)=U(x) represents the interaction potential and xi is the state vector of the ith subsystem. Weakly interacting diffusions are characterized by a coupling strength that is inversely proportional to the number of subsystems N. As N increases, the interaction structure gets more and more intricate, while the intensity becomes weaker and weaker (see Fig. 3). As mentioned in the main text, the DZ model has been introduced, and thereafter commonly used, as a paradigmatic example of an equilibrium continuous phase transition reminiscent of the Ising-like ferromagnetic transitions in spin systems.3,37 The DZ model describes a network of one-dimensional subsystems xk whose dynamics is prescribed by the following equations (see main text for notation):

dxk=Vα(xk)dtθNl=1N(xkxl)dt+σdWk,
(A1)

where k=1,,N and the confining potential Vα(x)=α2x2+x44 has a double well shape for α>0. Without loss of generality, we here consider α=1. In the absence of coupling, θ=0, the above equations describe the simple motion of a particle in a double well potential, subject to additive noise. The presence of the coupling allows for a long range coordination of the system that in the thermodynamic limit N+ results in a proper phase transition. In this regime, by varying the parameters (α,θ), the order parameter x undergoes a continuous order–disorder transition, similar to the pitchfork bifurcation diagram for the Ising model. It is possible to show3 that the critical line is given by D3/2(θ1σ)D1/2(θ1σ)=σθ, where Dν(z) is a parabolic cylinder function. Here, the coupling θ is kept fixed (θ=0.55) and we vary σ to approach the transition point.

The BCM model describes an ensemble of bi-dimensional nonlinear oscillators undergoing an out of equilibrium self-synchronization transition. The time evolution of the network of oscillators is given by the following equations:

dxk=Fα(xk)dtθNl=1N(xkxl)dt+σdWk,
(A2)

where the local force is not conservative, giving rise to the nonequilibrium features of the system, and reads Fα(x)=(α|x|2)x+x+ where x+=(x2,x1). The latter term corresponds to a rotation and makes the stationary state a nonequilibrium one. The parameter α>0 controls the amplitude of the oscillations of the individual nonlinear oscillators. In fact, when θ=σ=0, each subsystem oscillates as xj(t)=α(cos(t+βj),sin(t+βj)), where βj=tan(x2j(0)/x1j(0)). The coupling tries to synchronize the subsystems by attracting them toward the center of mass 1Nj=1Nxj. In the thermodynamic limit and for sufficiently low values of the noise, the order parameter x=x(t) exhibits a periodic time evolution, resulting from the subsystem oscillating in a coherent way. On the other hand, high values of the noise correspond to a nonsynchronized state where the order parameter vanishes. In particular, the transition happens at the surface of the (α,θ,σ) parametric space defined by46 

A=δ22[11δexp(A2δ2)[Aδer2dr]1],
(A3)

where A=αθ1 and δ=2σ2θ. In the following, we have θ=α=2. The color code for the figures is given by

  • noncritical black: DZ σ1 , BCM σ2,

  • noncritical blue: DZ σ0.87 , BCM σ1.8, and

  • critical red: DZ σ~0.75 , BCM σ~1.59.

1. Numerical linear response experiments

As mentioned in the main text, we perform n simulations of the response given by Eqs. (A1) and (A2) where the initial conditions are chosen according to the respective unperturbed invariant measure ρ0(x) and where at time τ=0, we apply a perturbation proportional to Dirac’s δ(τ=0). The average of the response for the observable x over the n simulations gives an estimate of G~i(τ;N). The response functions away from the transitions are estimated on an ensemble of n=105 simulations, while the critical response functions with n=106 for DZ and n=7×106 for BCM. Furthermore, we investigate the response up to time τ=5×103. The corresponding susceptibility Γ~i(ω;N) is simply defined as the Fourier transform of G~i(τ;N). In the main text, we show the results for an additive perturbation X(x)X. However, the critical behavior of the response does not depend on the type of perturbation, modulo, a potential degenerate class of perturbations that have zero projection on the invariant measure ρ0(x). We have thus decided to investigate the response of the DZ model for a spatially dependent perturbation X(x)=x2 (see Fig. 4). The response function G~i(τ;N), both away and at the phase transition, has a rapid initial decay with a timescale that is different from the response function shown in the main text. As a matter of fact, the timescale associated to the dominant mode of the response function for τ0 does, in general, depend on the applied perturbation.40 As expected, the response function at the phase transition develops a much longer timescale that increases as the number of particles increases. A more accurate comparison with the result shown in the main text can only be performed in the frequency domain. Figure 4 (right panel) shows that, away from the transition, the susceptibilities have a smooth behavior and no evident dependence on N. At the phase transition, the susceptibility develops the expected singular behavior κωω0+iγ(N), where γ(N)0+ as N+ due to the appearance of a simple pole ω0=0. The residue κ is purely imaginary and its magnitude k can be inferred by visual inspection of the top inset representing the function(ωω0)Im{Γ~(ω)} to be just less than 0.29. As mentioned in the main text, the residue depends both on the observable and on the perturbation X(x).

FIG. 4.

Panel (a): response function G~(τ;N) and panel (b): susceptibility Γ~(ω;N) for a spatially dependent perturbation X(x)=x2 for the one-dimensional order parameter for the DZ model. Color code and plotting conventions as in Fig. 1. The lower extreme of the integral of the bottom panel is c=0.05.

FIG. 4.

Panel (a): response function G~(τ;N) and panel (b): susceptibility Γ~(ω;N) for a spatially dependent perturbation X(x)=x2 for the one-dimensional order parameter for the DZ model. Color code and plotting conventions as in Fig. 1. The lower extreme of the integral of the bottom panel is c=0.05.

Close modal
1.
G.
Naldi
,
L.
Pareschi
, and
G.
Toscani
,
Mathematical Modeling of Collective Behavior in Socio-Economic and Life Sciences
(
Birkhäuser
,
Basel
,
2010
).
2.
L.
Pareschi
and
G.
Toscani
, Interacting Multiagent Systems: Kinetic Equations and Monte Carlo Methods (Oxford University Press, 2013).
3.
D. A.
Dawson
, “
Critical dynamics and fluctuations for a mean-field model of cooperative behavior
,”
J. Stat. Phys.
31
,
29
85
(
1983
).
4.
J. A.
Acebrón
,
L. L.
Bonilla
,
C. J.
Pérez Vicente
,
F.
Ritort
, and
R.
Spigler
, “
The Kuramoto model: A simple paradigm for synchronization phenomena
,”
Rev. Mod. Phys.
77
,
137
185
(
2005
).
5.
J.
Garnier
,
G.
Papanicolaou
, and
T.
Yang
, “
Large deviations for a mean field model of systemic risk
,”
SIAM J. Financial Math.
4
,
151
184
(
2013
).
6.
C.
Wang
,
Q.
Li
, and
W. E.
B. Chazelle
, “
Noisy Hegselmann–Krause systems: Phase transition and the 2R-conjecture
,”
J. Stat. Phys.
166
,
1209
1225
(
2017
).
7.
J.
Garnier
,
G.
Papanicolaou
, and
T.
Yang
, “
Consensus convergence with stochastic effects
,”
Vietnam J. Math.
45
,
51
75
(
2017
).
8.
A.
Garbuno-Inigo
,
N.
Nüsken
, and
S.
Reich
, “
Affine invariant interacting Langevin dynamics for Bayesian inference
,”
SIAM J. Appl. Dyn. Syst.
19
,
1633
1658
(
2020
).
9.
J.
Simmonds
,
J. A.
Gómez
, and
A.
Ledezma
, “
The role of agent-based modeling and multi-agent systems in flood-based hydrological problems: A brief review
,”
J. Water Climate Change
11
,
1580
1602
(
2019
). jwc2019108,
10.
S.
Geisendorf
, “
Evolutionary climate-change modelling: A multi-agent climate-economic model
,”
Comput. Econ.
52
,
921
951
(
2018
).
11.
M.
Scheffer
, Critical Transitions in Nature and Society, Princeton Studies in Complexity (Princeton University Press, Princeton, 2009).
12.
D.
Sornette
, “Endogenous versus exogenous origins of crises,” in Extreme Events in Nature and Society, edited by S. Albeverio, V. Jentsch, and H. Kantz (Springer, Berlin, Heidelberg, 2006) pp. 95–119.
13.
R.
Kubo
, “
The fluctuation–dissipation theorem
,”
Rep. Progress Phys.
29
,
255
284
(
1966
).
14.
U. M. B.
Marconi
,
A.
Puglisi
,
L.
Rondoni
, and
A.
Vulpiani
, “
Fluctuation–dissipation: Response theory in statistical physics
,”
Phys. Rep.
461
,
111
195
(
2008
).
15.
M.
Baiesi
and
C.
Maes
, “
An update on the nonequilibrium linear response
,”
New. J. Phys.
15
,
013004
(
2013
).
16.
A.
Sarracino
and
A.
Vulpiani
, “
On the fluctuation–dissipation relation in non-equilibrium and non-Hamiltonian systems
,”
Chaos
29
,
083132
(
2019
).
17.
V.
Lucarini
,
J. J.
Saarinen
,
K.-E.
Peiponen
, and
E. M.
Vartiainen
,
Kramers–Kronig Relations in Optical Materials Research
(
Springer
,
New York
,
2005
).
18.
J.
Binney
and
S.
Tremaine
,
Galactic Dynamics
, 2nd ed. (
Princeton University Press
,
Princeton
,
2008
).
19.
C. E.
Leith
, “
Climate response and fluctuation dissipation
,”
J. Atmos. Sci.
32
,
2022
(
1975
).
20.
G.
North
,
R.
Bell
, and
J.
Hardin
, “
Fluctuation dissipation in a general circulation model
,”
Clim. Dyn.
8
,
259
(
1993
).
21.
H.
Öttinger
,
Beyond Equilibrium Thermodynamics
(
Wiley
,
Hoboken
,
2005
).
22.
V.
Lucarini
,
F.
Ragone
, and
F.
Lunkeit
, “
Predicting climate change using response theory: Global averages and spatial patterns
,”
J. Stat. Phys.
166
,
1036
1064
(
2017
).
23.
B.
Cessac
, “
Linear response in neuronal networks: From neurons dynamics to collective response
,”
Chaos
29
,
103105
(
2019
).
24.
G. A.
Gottwald
, “
Introduction to focus issue: Linear response theory: Potentials and limits
,”
Chaos
30
,
020401
(
2020
).
25.
D.
Ruelle
, “
Nonequilibrium statistical mechanics near equilibrium: Computing higher-order terms
,”
Nonlinearity
11
,
5
18
(
1998
).
26.
D.
Ruelle
, “
A review of linear response theory for general differentiable dynamical systems
,”
Nonlinearity
22
,
855
870
(
2009
).
27.
V.
Baladi
, “Linear response, or else,” in ICM Seoul 2014, Proceedings (KYUNG MOON SA Co. Ltd., 2014), Vol. III, pp. 525–545.
28.
A.
Dembo
and
J.-D.
Deuschel
, “
Markovian perturbation, response and fluctuation dissipation theorem
,”
Ann. Inst. Henri Poincaré Probab. Stat.
46
,
822
852
(
2010
).
29.
M.
Hairer
and
A. J.
Majda
, “
A simple framework to justify linear response theory
,”
Nonlinearity
23
,
909
922
(
2010
).
30.
C. L.
Wormell
and
G. A.
Gottwald
, “
Linear response for macroscopic observables in high-dimensional systems
,”
Chaos
29
,
113127
(
2019
).
31.
C.
Liverani
and
S.
Gouëzel
, “
Banach spaces adapted to Anosov systems
,”
Ergodic Theory Dyn. Syst.
26
,
189
217
(
2006
).
32.
M. D.
Chekroun
,
J. D.
Neelin
,
D.
Kondrashov
,
J. C.
McWilliams
, and
M.
Ghil
, “
Rough parameter dependence in climate models and the role of Ruelle–Pollicott resonances
,”
Proc. Natl. Acad. Sci. U.S.A.
111
,
1684
1690
(
2014
).
33.
V.
Lucarini
, “
Response operators for Markov processes in a finite state space: Radius of convergence and link to the response theory for Axiom A systems
,”
J. Stat. Phys.
162
,
312
333
(
2016
).
34.
A.
Tantet
,
V.
Lucarini
, and
H. A.
Dijkstra
, “
Resonances in a chaotic attractor crisis of the Lorenz flow
,”
J. Stat. Phys.
170
,
584
616
(
2018
).
35.
M.
Pollicott
, “
On the rate of mixing of Axiom A flows
,”
Inventiones Math.
81
,
413
426
(
1985
).
36.
D.
Ruelle
, “
Resonances of chaotic dynamical systems
,”
Phys. Rev. Lett.
56
,
405
407
(
1986
).
37.
M.
Shiino
, “
Dynamical behavior of stochastic systems of infinitely many coupled nonlinear oscillators exhibiting phase transitions of mean-field type: H theorem on asymptotic approach to equilibrium and critical slowing down of order-parameter fluctuations
,”
Phys. Rev. A
36
,
2393
2412
(
1987
).
38.
M. G.
Delgadino
,
R. S.
Gvalani
, and
G. A.
Pavliotis
, “
On the diffusive-mean field limit for weakly interacting diffusions exhibiting phase transitions
,”
Arch. Ration. Mech. Anal.
2021
,
1
58
(
2021
).
39.
J. A.
Carrillo
,
R. S.
Gvalani
,
G. A.
Pavliotis
, and
A.
Schlichting
, “
Long-time behaviour and phase transitions for the McKean–Vlasov equation on the torus
,”
Arch. Ration. Mech. Anal.
235
,
635
690
(
2020
).
40.
V.
Lucarini
,
G. A.
Pavliotis
, and
N.
Zagli
, “
Response theory and phase transitions for the thermodynamic limit of interacting identical systems
,”
Proc. R. Soc. A
476
,
0688
(
2020
).
41.
J. D.
Jackson
,
Classical Electrodynamics
, 2nd ed. (
Wiley
,
New York
,
1975
).
42.
E.
Talebian
and
M.
Talebian
, “
A general review on the derivation of Clausius–Mossotti relation
,”
Optik
124
,
2324
2326
(
2013
).
43.
D.
Sornette
and
A.
Helmstetter
, “
Endogenous versus exogenous shocks in systems with memory
,”
Phys. A: Stat. Mech. Appl.
318
,
577
591
(
2003
).
44.
N.
Yoshida
, “
Phase transition from the viewpoint of relaxation phenomena
,”
Rev. Math. Phys.
15
,
765
788
(
2003
).
45.
R. C.
Desai
and
R.
Zwanzig
, “
Statistical mechanics of a nonlinear stochastic model
,”
J. Stat. Phys.
19
,
1
24
(
1978
).
46.
L. L.
Bonilla
,
J.
Casado
, and
M.
Morillo
, “
Self-synchronization of populations of nonlinear oscillators in the thermodynamic limit
,”
J. Stat. Phys.
48
,
571
591
(
1987
).
47.
A.
Sznitman
, “Topics in propagation of chaos,” Ecole d’Eté de Probabilités de Saint-Flour XIX—1989, Lecture Notes in Mathematics Vol. 1464, edited by P. L. Hennequin (Springer, Berlin, 1989).
48.
K.
Oelschlager
, “
A martingale approach to the law of large numbers for weakly interacting stochastic processes
,”
Ann. Probab.
12
,
458
479
(
1984
).
49.
D. A.
Dawson
and
J.
Gärtner
, “
Large deviations from the McKean–Vlasov limit for weakly interacting diffusions
,”
Stochastics
20
,
247
308
(
1987
).
50.
B.
Helffer
, Semiclassical Analysis, Witten Laplacians, and Statistical Mechanics, Series in Partial Differential Equations and Applications Vol. 1 (World Scientific Publishing Co., Inc., River Edge, NJ, 2002), p. x+179.
51.
T.
Frank
, “
Fluctuation–dissipation theorems for nonlinear Fokker–Planck equations of the Desai–Zwanzig type and Vlasov–Fokker–Planck equations
,”
Phys. Lett. A
329
,
475
485
(
2004
).
52.
M. D.
Chekroun
,
A.
Tantet
,
H. A.
Dijkstra
, and
J. D.
Neelin
, “
Ruelle–Pollicott resonances of stochastic systems in reduced state space. Part I: Theory
,”
J. Stat. Phys.
179
,
1366
1402
(
2020
).
53.
A.
Lasota
and
M. C.
Mackey
, Chaos, Fractals, and Noise, 2nd ed., Applied Mathematical Sciences Vol. 97 (Springer-Verlag, New York, 1994), p. xiv+472.
54.
U. M. B.
Marconi
,
A.
Puglisi
,
L.
Rondoni
, and
A.
Vulpiani
, “
Fluctuation–dissipation: Response theory in statistical physics
,”
Phys. Rep.
461
,
111
(
2008
).
55.
P.-H.
Chavanis
, “
The Brownian mean field model
,”
Eur. Phys. J. B
87
,
120
(
2014
).
56.
There is a typo in the formula given in Ref. 40. Furthermore, the convention for the Fourier transform we use here has the opposite sign, ditto the residue.
57.
N.
Zagli
(2021). “Spectroscopy of phase transitions,” Figshare. https://figshare.com/projects/Spectroscopy_of_phase_transitions/101846