Response theory and phase transitions for the thermodynamic limit of interacting identical systems

We study the response to perturbations in the thermodynamic limit of a network of coupled identical agents undergoing a stochastic evolution which, in general, describes non-equilibrium conditions. All systems are nudged towards the common centre of mass. We derive Kramers–Kronig relations and sum rules for the linear susceptibilities obtained through mean field Fokker–Planck equations and then propose corrections relevant for the macroscopic case, which incorporates in a self-consistent way the effect of the mutual interaction between the systems. Such an interaction creates a memory effect. We are able to derive conditions determining the occurrence of phase transitions specifically due to system-to-system interactions. Such phase transitions exist in the thermodynamic limit and are associated with the divergence of the linear response but are not accompanied by the divergence in the integrated autocorrelation time for a suitably defined observable. We clarify that such endogenous phase transitions are fundamentally different from other pathologies in the linear response that can be framed in the context of critical transitions. Finally, we show how our results can elucidate the properties of the Desai–Zwanzig model and of the Bonilla–Casado–Morillo model, which feature paradigmatic equilibrium and non-equilibrium phase transitions, respectively.


Introduction
Multi-agent systems are used routinely to model phenomena in the natural sciences, social sciences, and engineering. In addition to the standard applications of interacting particle systems to, e.g. plasma physics and stellar dynamics, phenomena such as cooperation [1], synchronization [2], systemic risk [3], consensus opinion formation [4,5] can be modelled using interacting multi-agent systems. Multi-agent systems are finding applications also in areas like management of natural hazards [6] and of climate change impacts [7]. We refer to [8] for a recent review on interacting multi-agent systems and their applications to the social sciences, and to [9] for a collections of articles showcasing their application in many different areas of science and technology. Additionally, multi-agent systems are also used as the basis for algorithms for sampling and optimization [10].
In this paper, we focus on a particular class of multi-agent systems, namely weakly interacting diffusions, for which the strength of the interaction between the agents is inversely proportional to the number of agents. Under the assumption of exchangeability, i.e. that the particles are identical, it is well known that one can pass to the limit as the number of agents goes to infinity, i.e. the mean field limit. In particular, in this limit the evolution of the empirical measure is described by a nonlinear, non-local Fokker-Planck equation, the McKean-Vlasov Equation [11,12]. We refer to [13] for a comprehensive review of the McKean-Vlasov equation from a theoretical physics viewpoint. The class of multi-agent models considered in this paper is sufficiently rich to include models for cooperation, systemic risk, synchronization, biophysics, and opinion formation.
An important feature of weakly interacting diffusions is that in the mean field (thermodynamic) limit they can exhibit phase transitions [1,4,[14][15][16]. Phase transitions are characterized in terms of exchange of stability of non-unique stationary states for the McKean-Vlasov equation at the critical temperature/interaction strength.
In the case of equilibrium systems, such stationary states are associated with critical points of a suitably defined energy landscape. For example, for the Kuramoto model of nonlinear oscillators, at the critical noise strength the uniform distribution (on the torus) becomes unstable and stable localized stationary states emerge (phase-locking), leading to synchronization phase transition [17]. A complete theory of phase transitions for the McKean-Vlasov equation on the torus, that includes the Kuramoto model of synchronization, the Hegselmann-Krause model of opinion formation, the Keller-Segel model of chemotaxis. etc., is presented in [18], see also [19]. The effect of (infinitely) many local minima in the energy landscape on the structure of the bifurcation diagram was studied in [20]. Phase transitions for gradient system with local interactions were studied in [21,22]. Synchronization has been extensively discussed in the scientific literature, see [17,[23][24][25][26][27][28].

(a) Linear response theory
One of the main objectives of this paper is to investigate phase transitions for weakly interacting diffusions by looking at the response of the (infinite dimensional) mean field dynamics to weak external perturbations. We associate the nearing of a phase transition with the setting where a very small cause leads to very large effects, or, more technically, to the breakdown of linear response in the system, as described below.
Linear response theory provides a general framework for investigating the properties of physical systems [29]. Well-known applications of linear response theory include solid-state physics and optics [30] as well as plasma physics and stellar dynamics ( [31], ch. 5). Furthermore, the range of systems for which linear response theory is relevant is very vast, see e.g. [32][33][34][35][36]. Recently, many new areas of applications of linear response theory are emerging across different disciplinary areas-see, e.g. a recent special issue [37]-and new formulations of the problem are being presented, where the conceptual separation between acting forcing and observed response is blurred [38]. In particular, recent applications of linear response theory include the prediction of climate response to forcings [39][40][41][42][43][44]. In modern terms, the goal is to define practical ways to reconstruct the measure supported time-dependent pullback attractor [45] of the climate by studying the response to perturbations of a suitably defined reference climate state [46].
The mathematical theory of linear response for deterministic systems was developed by Ruelle in the context of Axiom A chaotic systems [47,48]. He provided explicit response formulae and showed that, in the case of dissipative systems, the classical fluctuation-dissipation theorem does not hold, and, as a result, natural and forced fluctuations are intimately different [49]. Ruelle's results have then been re-examined through a more a functional analytic lens by studying the impacts of the perturbations to the dynamics on the transfer operator [50] and then extended to a more abstract mathematical framework [51][52][53]. The direct implementation of Ruelle's formulae is extremely challenging, because of the radically different behaviour of the system along the stable and unstable manifold [54], which is related to the insightful tout court criticism of linear response theory by Van Kampen [55], so that alternative strategies have been devised [56,57]. Very promising progresses have been recently obtained in the direction of using directly Ruelle's formulae thanks to adjoint and shadowing methods [58][59][60].
Linear response theory and fluctuation-dissipation theorems have long been studied in detail for diffusion processes ( [61], [62, ch. 7], [63], ch. 9), and, more recently, rigorous results have been obtained in this direction [64,65]. An interesting link between response theory for deterministic and stochastic systems has been proposed in [66]. The results presented in [64,65] can be applied to the McKean-Vlasov equation in the absence of phase transitions to justify rigorously linear response theory and to establish fluctuation-dissipation results. See also [67] for formal calculations. This is not surprising, since it is well known that, in the absence of phase transitions, fluctuations around the mean field limit are Gaussian and can be described in terms of an appropriate stochastic heat equation [1,68].

(b) Critical transitions versus phase transitions
Critical transitions appear when the spectral gap of the transfer operator [51] of the unperturbed system becomes vanishingly small, as a result of the Ruelle-Pollicott poles [69,70] touching the real axis. Since there is a one-to-one correspondence between the radius of expansion of linear response theory and the spectral gap of the transfer operator [51,71], near critical transitions the linear response breaks down and one finds rough dependence of the system properties on its parameters [72,73]. Systems undergoing critical transitions appear often in the natural and social sciences [74] and a lot of effort has been put in the development of early warning signals for critical transitions [75][76][77]. Early warning signals include an increase in variance and correlation time as the system approaches the transition point.
In the deterministic case, at the critical transition the reference state loses stability and the system ends up in a possibly very different metastable state. Indeed, the presence of critical transitions is closely related to the existence of regimes of multi-stability [78,79]. Transition points for finite-dimensional stochastic systems correspond to points where the topological structure of the unique invariant measure changes ( [80], [63,Sect. 5.4]). Contrary to this, more than one invariant measure can exist in the mean field (thermodynamic) limit, e.g. in equilibrium systems, when the free energy is not convex [18]. The loss of uniqueness of the stationary state at the critical temperature/noise strength corresponds to a phase transition.
Phase transitions are usually defined by (a) identifying an order parameter and (b) verifying that in the thermodynamic limit, for some value of the parameter of the system, the properties of such an order parameter undergo a sudden change. It should be emphasized, however, that, for the mean field dynamics it is not always possible to identify an order parameter. The way we define phase transitions in this work comes from a somewhat complementary viewpoint, which aims at clarifying analogies and differences with respect to the case of critical transitions.
Sornette and collaborators have devoted efforts at separating the effects of endogeneous versus exogenerous processes in determining the dynamics of a complex system and, especially in defining the conditions conducive to crises [81], and proposed multiple applications in the natural (e.g. [82]) as well as the social (e.g. [ response of the system to exogeneous perturbations and the decorrelation due to endogenous dynamics is interpreted as resulting from a fluctuation-dissipation relation-like properties. Finally, Sornette and collaborators have also emphasized the importance of memory effects especially in the context of endogenous dynamics [84,85]. While our viewpoint and methods are different from theirs, what we pursued here shares similar goals and delves into closely related concepts.

(c) This paper: goals and main results
The main objective of this paper is to perform a systematic study of linear response theory for mean field partial differential equations (PDEs) exhibiting phase transitions. Indeed, it has been shown that, for nonlinear oscillators coupled linearly with their mean, the so-called Desai-Zwanzig model [86], the fluctuations at the phase transition point are not Gaussian [1], see also [19] for related results for a variant of the Kuramoto model (the Haken-Kelso-Bunz model). Indeed, the fluctuations are persistent, non-Gaussian in time, with an amplitude described by a nonlinear stochastic differential equation, and associated with a longer timescale [1]. At the transition point, the standard form of linear response theory breaks down [14]. More general analyses performed using ideas from linear response theory of how a system of coupled maps performs a transition to a coherent state in the thermodynamic limit can be found in [87,88].
Here, we consider a network of N identical and coupled M-dimensional systems whose evolution is described by a Langevin equation. We then study the response to perturbations in the limit of N → ∞. We investigate the conditions determining the breakdown of the linear response and separate two possible scenarios. One scenario pertains to the closure of the spectral gap of the transfer operator of the mean field equations, and can be dealt with through the classical theory of critical transitions. A second scenario of breakdown of the linear response results from the coupling among the N systems and is inherently associated with the thermodynamic limit. We focus on the second scenario of breakdown of the linear response, which we interpret as corresponding to a phase transition. The main results of this paper can be summarized as follows: -the derivation of linear response formulae for the thermodynamic limit of a network of coupled identical systems and of Kramers-Kronig relations and sum rules for the related susceptibilities; -the statement of conditions leading to phase transitions as opposed to the classical scenario of critical transitions; -the explicit derivation of the corrections to the standard Kramers-Kronig relations and sum rules occurring at the phase transition; -the clarification, through the use of functional analytical arguments, of why one does not expect divergence of the integrated autocorrelation time of suitable observables in the case of phase transitions, whereas the opposite holds in the case of critical transitions; -the re-examination, also through numerical simulations, of classical results on phase transitions in the Desai-Zwanzig model [86] and in the Bonilla-Casado-Morillo model [89].
The rest of the paper is organized as follows. In §2, we introduce our model and present the linear response formulae for the mean field equations as well as for the renormalized macroscopic case. In §3, we discuss the properties of the frequency-dependent susceptibility, present the Kramers-Kronig relations connecting their real and imaginary parts, and find explicit sum rules. In §4 we discuss under which conditions the response diverges, and clarify the fundamental difference between the case of critical transitions and the case of phase transitions, which can take place only in the thermodynamic limit. Section 5 is dedicated to finding results that specifically apply to the case of gradient systems, corresponding to reversible Markovian dynamics. In §6 we re-examine the case of phase transitions for the Desai-Zwanzig and Bonilla-Casado-Morilla models, which are relevant for the case of equilibrium and non-equilibrium dynamics, respectively. Finally, in §7, we present our conclusions and provide perspectives for future investigations.

Linear response formulae: mean field and macroscopic results
We consider a network of N exchangeable interacting M-dimensional systems whose dynamics is described by the following stochastic differential equations: where F α is a smooth vector field, possibly depending on a parameter α. Additionally, dW i , i = 1, . . . , N are independent Brownian motions (the Ito convention is used); s ij is the volatility matrix, and the parameter σ > 0 controls the intensity of the stochastic forcing. Additionally, the N systems undergo an all-to-all coupling through the Laplacian matrix given by the derivative of the potential U(y) = |y| 2 /2. We emphasize the fact that the linear response theory calculations presented below are valid for arbitrary choices of the interaction potential. We choose to present our results for the case of quadratic interactions since in this case the order parameter is known; furthermore, the stationary state(s) are known and are parametrized by the order parameter [1]. The coefficient θ modulates the intensity of such a coupling, which attempts at synchronizing all systems by nudging them to the centre of mass 1/N N k=1 x k . If θ = 0, the N systems are decoupled. We remark that the theory of synchronization says that for this choice of the coupling, if dx/dt = F α (x) has a unique attractor and is chaotic with λ 1 > 0 being the largest Lyapunov exponent, the N nodes undergo perfect synchronization for any N ≥ 2 in the absence of noise (σ = 0) if θ > λ 1 [17,28,90,91].
If F α (y) = −∇V α (y), we interpret V α as the confining potential [63]. In some cases, equation (2.1) describes an equilibrium statistical mechanical system, in particular if F α = −∇V α (y) and s ij is proportional to the identity. More generally, equilibrium conditions are realized when the drift term-the deterministic component on the right-hand side of equation (2.1)-is proportional to the gradient of a function defined according to the Riemannian metric given by the diffusion matrix C ij = s ik s jk [92].
We now consider the empirical measure ρ (N) , which is defined as ρ (N) = 1/N N k=1 δ x k (t) . Following [11,18,93], we investigate the thermodynamic limit of the system above. As N → ∞, we can use martingale techniques [1,11,93,94] to show that the one-particle density converges to some measure ρ(x, t) satisfying the following McKean-Vlasov equation, which is a nonlinear and non-local Fokker-Planck equation : where we have separated the linear operator L 0 α,θ and the nonlinear operator Λ θ ({ρ(x, t)}) = θ ∇ · (ρ(x, t) x (t)), with x (t) = d M yρ(y, t) and denotes the convolution. Additionally, we have that˜ is a linear diffusion operator such that˜ ρ( t)), which coincides with the standard M-dimensional Laplacian (˜ = ) if the diffusion matrix C ij is the identity matrix. If σ = 0, we are considering a nonlinear Liouville equation. We assume that, if σ > 0, equation (2.1) describes a hypoelliptic diffusion process, so that ρ [63, ch. 6]. In what follows, we refer to the case σ > 0. Conditions detailing the well-posedness of this problem can be found in [95].
Let us define ρ (0) α,θ (x) as a reference invariant measure of the system such that L 0 α,θ (ρ (0) α,θ (x)) + θΛ θ ({ρ (0) α,θ (x)}) = 0. Since we are considering a system with an infinite number of particles, such an invariant measure needs not be unique [1,18,96,97]. Specifically, if s ij is proportional to the identity and F α (y) = −∇V α (y) and V α (y) is not convex, thus allowing for more than one local minimum, for a given value of θ the system undergoes a phase transition for sufficiently weak noise; see discussion in §5.
We remark that the invariant measure depends on the values of α and θ, and, in particular, α,θ = x 0 is a constant vector, where in the last identity we have dropped the lower indices to simplify the notation. As a result, we have that: is the eigenvector with vanishing eigenvalue of the linear operator M 0 α,θ, x 0 . Taking inspiration from [98,99], we now study the impact of perturbations on the invariant measure ρ (0) α,θ (x). We follow and extend the results presented in [67]. We modify the right-hand side of equation (2.2) by setting F α (x) → F α (x) + X(x)T(t) and we study the linear response of the system in terms of the density ρ(x, t). We then write ρ( ( 2 ) and obtain the following equation up to order : We remark that the linear operatorM 0 α,θ, x 0 acting on ρ (1) α,θ (x, t) on the right-hand side of the previous equation is not the operator whose zero eigenvector is the unperturbed invariant measure. The correction proportional to θ emerges as a result of the nonlinearity of the McKean-Vlasov equation. We will discuss the operatorM 0 α,θ, x 0 in §2a. One then derives: We now evaluate the response of the observable x i . This is sufficient for our purposes, since we know that, for this model, the order parameter (in the mean field limit) is the mean position x (magnetization). By definition, we have that α,θ (y, t)Φ(y) for a generic observable Φ. We obtain: where we have defined the following operator: where O + is the adjoint of O. Following [67], we can interpret this as the Koopman operator for the unperturbed dynamics; see later discussion. We can rewrite the previous expression as: where the Green function is causal. Note also that if X(x) =v k , wherev k is the unit vector in the kth direction, then Notwithstanding the Markovianity of the dynamics, the second term on the right-hand side of equation (2.8) describes a memory effect in the response of the observable x. Such a term emerges in the thermodynamic limit and effectively imposes a condition of self-consistency between forcing and response; see different yet related results obtained by Sornette and collaborators [81,84,85].
If σ > 0 the invariant measure is smooth, so that we can perform an integration by parts of the previous expressions and derive the following Green functions: where the Green functions are written as correlation functions times a Heaviside distribution enforcing causality. We remark that we can, at least formally, write: where {λ j } ∞ j=1 are the eigenvalues (point-spectrum) of M 0,+ α,θ, x 0 and Π j is the spectral projector onto the eigenspace spanned by the eigenfunction ψ j , and in particular, Π 0 projects on the invariant measure. Then, the operator R(t) is the residual operator associated with the essential spectrum. The norm of R(t) is controlled by the distance of essential spectrum from the imaginary axis.
We then have: and α,θ (y)). Note that the j = 0 term vanishes because the corresponding scalar product Φ X ψ 0 0 has nil value for any choice of the vector field X. We now apply the Fourier transform F to equation (2.8) and obtain: where we have used a (standard) abuse of notation in defining the Fourier transform of T(t) and x j 1 (t) and have defined and We remark that the susceptibilities given in equations (2.17) and (2.18) are holomorphic in the Note that all susceptibilities, regardless of the observable considered, share the same poles located at ω j = iλ j , j = 1, . . . , ∞. Additionally, if ω j is a pole, so is also −ω * j (and, correspondingly, λ j comes together with λ * j ). By introducing the inverse matrix Π α,θ = P −1 α,θ , we obtain from equation (2.16) our final result: The previous expression generalizes previous findings presented in [87]. We will discuss below the invertibility properties of the matrix P ij,α,θ (ω). If the coupling is absent, so that θ = 0, we obtain the same result as in the case of a single particle N = 1 system: Additionally, we trivially get Γ i,α,θ=0 (ω) =Γ i,α,θ=0 (ω). The effect of switching on the coupling and taking θ > 0 is twofold in terms of response: -First, the function Γ i,α,θ (ω) is modified, because the unperturbed evolution operator M 0 α,θ, x 0,θ (see equation (2.4)) and the unperturbed invariant measure ρ (0) α,θ (x) depend explicitly on θ. Indeed, changes in the value of θ impact expectation values and correlation properties. From the definition of M 0 α,θ, x 0 , we interpret Γ i,α,θ (ω) as the mean field susceptibility. -More importantly, the presence of a non-vanishing value of θ introduces a non-trivial correction with respect to the identity to the matrix P ij,α,θ (ω). We can interpret the functionΓ i,α,θ (ω) as the macroscopic susceptibility, which takes fully into account, in a self-consistent way, the interaction between the systems. Equation (2.19) generalizes the frequency-dependent version of the well-known Clausius-Mossotti relation [30,100,101], which connects the macroscopic polarizability of a material and the microscopic polarizability of its elementary components.
The integration by parts used for deriving equations (2.11) and (2.12) from equations (2.9) and (2.10) amounts to deriving a variant of the fluctuation-dissipation relation [29,33], as the Green functions are written as the causal part of a time-lagged correlation of two observables as determined by unperturbed dynamics. In other terms, the poles ω j , j = 1 . . . , ∞ of the susceptibilities above correspond to the Ruelle-Pollicott poles [69,70] of the unperturbed system, just as in the case of systems described by the standard Fokker-Planck equation [73,102]. This establishes a close connection between forced and free variability or, using a different terminology, between the properties of response to exogenous perturbations and endogenous dynamics [81].

(a) Another expression for the macroscopic susceptibility
A somewhat unsatisfactory aspect of the previous derivation resides in the fact that we are dealing with the operator exp(M 0,+ α,θ, x 0 t), which is associated with the mean field approximation. We can instead proceed from equation (2.5) using the operator exp(M 0 α,θ, x 0 t) introduced above and derive directly the following results: and We can rewrite the previous expression as: where the Fourier transform of: is the macroscopic susceptibility introduced in equation (2.19). Note thatM 0,+ α,θ, x 0 cannot be interpreted as the generator of time translation for smooth observables.
Clearly, the benefit of deriving the expression ofΓ i,α,θ (ω) as done in the previous section lies in the possibility of bypassing the space-integral operator included in the definition ofM 0,+ α,θ, x 0 . Similar to equation (2.13), we can write: where the corresponding symbols are used. We then have: We now apply the Fourier transform to equation (2.26) and obtain: Comparing equations (2.27) and (2.20), it is clear that the polesω j ofΓ i,α,θ (ω) are those of Γ i,α,θ (ω) plus those of the matrix Π ij,α,θ (ω), see earlier comments by Dawson [1] for the case of the Desai-Zwanzig model [86] (see also §6a).

Dispersion relations far from criticalities
We assume that all λ j , j = 1, . . . , ∞ have negative real part. As discussed above, since G Let us now consider the short-time behaviour τ → 0 + of the response functions G i,α,θ (τ ). Using equations (2.7) and (2.9), we derive: As a result, the high-frequency behaviour of the susceptibility Γ i,α,θ (ω) can be written as: The causality of G i,α,θ (τ ) implies that, using an abuse of notation, . By performing the Fourier transform of both sides of this identity, we obtains the following identity where indicates the convolution product andΘ(ω) = −iP(1/ω) + πδ(ω) is the Fourier transform of Θ(τ ), with P indicating the principal part. By separating the real (Re) and imaginary (Im) parts of Γ i,α,θ (ω), the previous relation can be written as: Since Γ i,α,θ (τ ) is a real function of real argument τ , its Fourier transform obeys the following conditions: We derive an alternative form of the Kramers-Kronig relations [30]: and It is then possible to derive the following sum rules: where is a measure of the decorrelation of the system, see a related result in [103] on the Desai-Zwanzig model [86] discussed below. Note that Im{Γ i,α,θ (ω)} is an odd function of ω. Additionally, if X i (x) 0 = 0, so that the imaginary part of the susceptibility decreases asymptotically at least as fast as ω −3 , the following additional sum rules holds: Let us now look at the asymptotic properties for large values of ω of the matrix P ij,α,θ (ω). We proceed as above and consider the short time behaviour of Y {i,j},α,θ (τ ): (3.10) As a result, for large values of ω, we have that so that P ij,α,θ (ω) = δ ij (1 − i(θ/ω)) + o(ω −2 ) and Π ij,α,θ (ω) = δ ij (1 + i(θ/ω)) + o(ω −2 ), so that: where we note a correction in the asymptotic behaviour with respect to the case of the mean field susceptibility given in equation (3.2). Nonetheless, if P ij,α,θ has full rank for all values of ω in the upper complex ω-plane, the Kramers-Kronig relations (3.5) and (3.6) and the sum rules (3.7)-(3.9) apply also for the macroscopic susceptibilitiesΓ i,α,θ (ω).

Criticalities
We remark again that the dispersion relations presented above apply for the mean field susceptibilities for values of α and θ such that (i) the real part of all the eigenvalues of M 0,+ α,θ, x 0 is negative; and for the macroscopic susceptibility if, additionally, (ii) the matrix P ij,α,θ is invertible and, additionally, has no zeros in the upper complex ω-plane. Conditions (i) and (ii) correspond to the case where the real part of all the eigenvalues ofM 0,+ α,θ, x 0 is negative. The breakdown of condition (i) for, say, (α, θ) = (ᾱ,θ) is due to the presence of a vanishing spectral gap for the operator M 0,+ α,θ, x 0 , and, a fortiori, for the operatorM 0,+ α,θ, x 0 . In such a scenario, the functions Γ i,ᾱ,θ (ω) andΓ i,ᾱ,θ (ω) feature one or more poles in the real ω-axis. In other terms, linear response blows up for forcings having non-vanishing spectral power |T(ω)| 2 at the corresponding frequencies.
In this case, because of the link discussed above between the poles of the mean field susceptibilities and the Ruelle-Pollicott poles of the unperturbed system, the blow-up of the linear susceptibilities corresponds to an ultraslow decay of correlations leading to a singularity in the integrated decorrelation time. In other terms, in this case the results conform to the classic framework of the theory of critical transitions [46,72,73,76,104]. We remark that the presence of a divergence does not depend on the specific functional form of the perturbation field X, while the properties of the response do depend in general from it.
The breakdown of condition (ii) for, say, (α, θ) = (α,θ) is associated with the fact that the spectral gap of the operatorM 0,+ α,θ, x 0 vanishes, while the spectral gap of the operator M 0,+ α,θ , x 0 remains finite. In this latter case, only the functionsΓ i,α,θ (ω) have one or more poles for real values of ω, whereas the functions Γ i,α,θ (ω) are holomorphic in the upper complex ω-plane. We remark that the non-invertibility of the P matrix depends on the presence of sufficiently strong coupling between the systems, which leads to them being coordinated, as discussed in detail in §6.
The nonlinearity of equation (2.2) emerges as a result of the thermodynamic limit N → ∞. Therefore, we interpret the singularities in the linear response resulting from the breakdown of condition (ii) as being associated to a phase transition of the system, yet not a standard one. Indeed, the blow-up of the linear susceptibilities does not correspond to a blow-up of the integrated correlation time (see §6a).

(a) Phase transitions
In what follows, we focus on the criticalities associated with condition (ii) only, which emerge specifically from effects that cannot be described using the mean field approximation.
Let's then assume that for some reference values for α = α 0 and θ = θ 0 the system is stable. This corresponds to the fact that the inverse Fourier transform ofΓ i,α 0 ,θ 0 (ω), which defines a renormalized linear Green function that takes into account all the interactions among the identical systems, has only positive support. Correspondingly, the macroscopic susceptibilitiesΓ i,α 0 ,θ 0 (ω), just like the mean field ones, are holomorphic in the upper complex ω-plane. This implies that the entries of the matrix Π ij,α 0 ,θ 0 (ω) do not have poles in the upper complex ω-plane.
We have that P ij,α,θ does not have full rank for ω = ω l , l = 1, . . . , R. For such value(s) of ω, the macroscopic susceptibilities diverge. Indeed, we remark that the invertibility conditions of the matrix P ij,α,θ (ω) is intrinsic and does not depend on the applied external forcing X, which enters, instead, only in the definition of the mean field susceptibility Γ i,α,θ (ω). We interpret this as the fact that the divergence of the response is due to eminently endogenous, rather than exogeneous, processes.
We also remark that P ij,α,θ (ω) = δ ij − Υ {i,j},α,θ (ω), where Υ {i,j},α,θ (ω) can be seen as mean field susceptibility for the expectation value of x i associated with an infinitesimal change of the value of the j th component of x 0 , see equations (2.4) and (2.10). This supports the idea that x is a appropriate order parameter for the system.
We assume, for simplicity, that only simple poles are present. We then decompose the matrix Π ij,α,θ (ω) in the upper complex ω-plane as follows: where we have separated the holomorphic component Π h ij,α,θ (ω) from the singular contributions coming from the poles ω l , l = 1, . . . , R; note that Res(f (ω)) ω=ν indicates the residue of the function f for ω = ν. Note that if ω l is a pole on the real axis, −ω l is also a pole. Additionally, Res(f (ω)) ω=ω l = −Res(f (ω)) * ω=−ω l , so that if ω l = 0 the residue has vanishing real part. Building on equation (4.1), the macroscopic susceptibility can then be written as: where the Kramers-Kronig relations given in equation (3.3) are then modified as follows, taking into account the extra poles along the real ω-axis: Res(Π ij,α,θ (ω)) ω=ω l ω l − ω Γ i,α,θ (ω l ). (4.3) By taking the limit ω → ∞ we can generalize the sum rule given in equation (3.7): Instead, by taking the limit ω → 0 we can generalize the sum rule given in equation (3.8) as follows: where we note that the zero-frequency poles do not contribute to the second term on the righthand side.

(b) Two scenarios of phase transition
In the discussion above, we are assuming that for (α, θ) = (α,θ) we have that det(P ij,α,θ (ω)) vanishes for R real values of ω. Since P ij,α,θ (ω) = (P ij,α,θ (−ω * )) * , we have that det(P ij,α,θ (ω)) = (det(P ij,α,θ (−ω * ))) * . Therefore, the solutions to the equation det(P ij,α,θ (ω)) = 0 come in conjugate pairs if they are complex. Generically, we can assume that as we tune the parameter s to the critical values such that (α˜s, θ˜s) = (α,θ ) either one real solution or the real part of one pair of solutions crosses to positive values. We then consider the following two scenarios for the poles ω l , l = 1, . . . , R: - Indeed, we wish to consider the two qualitatively different cases of either (i) a single pole with zero frequency; or (ii) a pair of poles with non-vanishing and opposite frequencies emerging at (α, θ ) = (α,θ). Of course, more than two poles could simultaneously emerge (α, θ) = (α,θ), but we consider this as a non-generic case.
-If ω l = 0 is a pole, then we have a static phase transition, associated with a breakdown in the linear response describing the parametric modulation of the measure of the system, see §6a. While such a statement applies for rather general systems and perturbations, this situation can be better understood by considering the specific perturbation X(x) = x 0 − x with T(t) = 1, which amounts to studying, within linear approximation, how the measure of the system changes as the value of θ is changed to θ + . This phase transition corresponds to a insulator-metal phase transition in condensed matter, because the electric susceptibility χ (1) ij (ω) of a conductor diverges as iσ ij /ω for small frequencies, where σ is a real tensor and describes the static electric conductivity, which is vanishing for an insulator [30].
-If, instead, we have a pair of poles located at ±ω l = 0, we have a dynamic phase transition activated by a forcing with non-vanishing spectral power at the frequency ±ω l . In this case, a limit cycle emerges corresponding to self-sustained oscillation, which is made possible by the feedback encoded in the nonlinearity of the McKean-Vlasov equation, see e.g. [89] and §6b.
In §6, we will present examples of phase transitions occurring according to the two scenarios above.

Equilibrium phase transitions: gradient systems
When the local force can be written as a gradient of a potential F α (y) = −∇V α (y) and the diffusion matrix is the identity matrix s ij = δ ij , equations (2.1) describe an equilibrium system. In particular, the N particles system has a unique ergodic invariant measure when the potential satisfies suitable confining properties [16,63] (see later discussion). Equivalently, the generator of the finite particle stochastic process has purely discrete spectrum, a non-zero spectral gap and the system converges exponentially fast to the unique equilibrium state, both in the L 2 space weighted by the invariant measure and in relative entropy. In the limit N → ∞, the system is described by the McKean-Vlasov equation (2.2) whose stationary measures are solutions of the Kirkwood-Monroe equation [105]: When the confining and interaction potentials are strongly convex and convex, respectively, then it is well known that equation (5.1) has only one solution, corresponding to the unique steady state of the McKean-Vlasov dynamics [106]. In addition, the dynamics converges exponentially fast, in relative entropy, to the stationary state and the rate of convergence to equilibrium can be quantified [106]. However, when the confining potential is not convex, e.g. is bistable, then more than one stationary states can exist, at sufficiently low noise strength (equivalently, for sufficiently strong interactions). A well-known example where the non-uniqueness of the invariant measure is that of the Desai-Zwanzig model [1,86,103], where the interaction potential is quadratic (see §6a for more details). In this framework, the loss of uniqueness of the invariant measure can be interpreted as a continuous-phase transition, similar to some extent to the phase transition for the Ising model. For a quadratic interaction potential, the equilibrium stationary measure (5.1) can be written as where we have introduced the modified potentialV( It is well known [63,Sect. 4.5] that, if the modified potentialV satisfies the property then the operator M 0 α,θ, x 0 in (5.3) has a spectral gap in L 2 (ρ 0 α,θ ), the space of square integrable functions weighted with by the invariant density ρ 0 α,θ . In particular, condition (5.4) prevents the system from undergoing a phase transition via scenario (i). When detailed balance holds, the mean field susceptibility G i,α,θ (τ ) relative to a uniform spatial forcing X = const. can be written as the time derivative of suitable correlation functions. In fact, from equation (2.11), the mean field susceptibility can be written as α,θ (y)X(y) . (5.5) Without loss of generality, let us consider an uniform forcing X =v k , withv k being the unit vector in the kth direction. The mean field susceptibility thus becomes Since the system is at equilibrium and the stationary probability density can be written as in (5.1), α,θ ∂ y kV , physically representing the fact that the probability current associated with the invariant measure vanishes at equilibrium. Furthermore, using (5.3) it is easy to verify the following identity M 0 α,θ, x 0 (y k ρ α,θ ∂ y kV . The mean field susceptibility can then be written as Note that this time scale differs from the one introduced in equation (3.8), which in this case can be written as By comparing the expressions of τ G i and τ ij and by considering equation (2.15), one understands that τ G i and τ ij correspond to two differently weighted averages of the timescales associated with each subdominant mode of the operator M 0 α,θ, x 0 . Usually, the singular behaviour of correlation properties has been used as an indicator of critical transitions [74]. However, let us remark again that, being related to the spectrum of the operator M 0 α,θ, x 0 , in our case neither τ G i nor τ ij show any critical behaviour at transitions occurring according to the scenario (ii), while they both diverge in the case of critical transitions corresponding to the scenario (i).

Examples
In what follows we re-examine the linear response of two relevant models that have been extensively investigated in the literature. Using the framework developed above, we investigate the phase transitions occurring in the Desai-Zwanzig model [86] and the Bonilla-Casado-Morilla model [89], which are taken as paradigmatic examples of equilibrium and non-equilibrium systems, respectively. We also provide the result of numerical simulations for both models.

(a) Equilibrium phase transition: the Desai-Zwanzig model
The Desai-Zwanzig model [86] has a paradigmatic value as it features an equilibrium thermodynamic phase transition (pitchfork bifurcation) arising from the interaction between systems [107] and has been used also as a model for systemic risk [3]. Each of the systems can be interpreted as a particle, moving in one dimension (M = 1) in a double well potential V α (x) = −(α/2)x 2 + x 4 /4, interacting with the other particles via a quadratic interaction U(x). The N− particle system is described by where k = 1, . . . , N. The local force is F α = −V α , the interaction potential is U(x) = x 2 /2 and the volatility matrix is the identity matrix s ij = δ ij . Furthermore, V α is double-well-shaped when α > 0, otherwise it has a unique global minimum. In the thermodynamic limit N → ∞, the one particle density satisfies the McKean-Vlasov equation (2.2) and it has been proven [1,103] that the infinite particle system undergoes a continuous phase transition, with x being a suitable order parameter. The Desai-Zwanzig model can be seen as a stochastic model of key importance for elucidating order-disorder phase transitions [107].
We have studied the Desai-Zwanzig model also through numerical integration of equation (6.1) by adopting an Euler-Maruyama scheme [108]. We have tested the convergence of our results in the thermodynamic limit N → ∞ by looking at increasing values of the number N of particles. We present in figures 1a-c the results obtained with N = 5000 for 0. equation (2.8). Given the simplicity of this model, it is possible to explicitly evaluate all the relevant quantities that characterize a phase transition relative to scenario (ii). Indeed, if we consider a perturbation such that X(x) = 1, one has Y α,θ (τ ) = θG α,θ (τ ) and equation (2.16) can be written as P(ω) x 1 (ω) = Γ α,θ (ω)T(ω) (6.2) where the 1 × 1 matrix is P(ω) = 1 − θΓ α,θ (ω). The macroscopic susceptibility is then obtained as Furthermore, this is a gradient system satisfying all the assumptions that have been made in §5, so that the mean field susceptibility can be written as (see also [103]) where z(t) = x − x 0 . Taking the Fourier transform results in where γ (ω) = ∞ 0 dte −iωt z(t)z(0) 0 is the Fourier transform of the correlation function. As previously mentioned, Γ α,θ (ω) can be written in terms of the spectrum of the operator M 0,+ α,θ, x 0 which in this specific example reads (see equation (5.3)) where the modified potential isV = V α − θ (x 2 /2 − x 0 x). It can be proven [1] that the above operator is self-adjoint and has a pure point spectrum {λ μ } with 0 = λ 0 > λ 1 > λ 2 > . . . , with the vanishing eigenvalue corresponding to the stationary distribution ρ (0) α,θ . In fact, it is easy to show that condition (5.4) holds. The operatorM 0,+ α,θ, x 0 is instead Dawson [1] proved that, away from the transition point-in particular, above it, where x 0 = 0the nonlinear operatorM 0,+ α,θ, x 0 has similar spectral properties to M 0,+ α,θ, x 0 . At the transition, though,M 0,+ α,θ, x 0 shows a vanishing spectral gap, with the operator developing a null eigenvalue. This situation corresponds to the breakdown of the aforementioned condition (ii) in which the mean field susceptibility Γ α,θ (ω)-and thus γ (ω)-is holomorphic in the upper complex ω-plane, while the macroscopicΓ α,θ develops a pole, arising from the non-invertibility of P(ω). Let us observe again that this implies that at the transition there is no divergence of the integrated autocorrelation time τ , because the spectral gap of the operator M 0,+ α,θ, x 0 does not shrink to zero. This is clearly shown in the two-dimensional map shown in figure 1c and in the two sections shown in figure 2a,b. We can fully characterize the singular behaviour of the macroscopic susceptibilityΓ α,θ at the transition. As a matter of fact, the transition point is characterized [103] by the condition so that the macroscopic susceptibility becomes As previously discussed in relation to equation (4.2), the above expression shows that at the transition pointΓ α,θ develops a simple pole in ω = 0, with residue (b) Non-equilibrium phase transition: the Bonilla-Casado-Morilla model In this section, we will study the Bonilla-Casado-Morrillo model [89] and elucidate the properties of a non-equilibrium self-synchronization phase transition, by looking at the divergence of the macroscopic susceptibilityΓ α,θ . We anticipate that the susceptibility develops a pair of symmetric poles ω 1 = −ω 2 > 0 at the transition point, thus following the scenario (ii) discussed above. The model consists of N two-dimensional nonlinear oscillators x k = (x k 1 , x k 2 ), interacting via a quadratic interaction potential U(x) = |x| 2 /2 and subjected to thermal noise The local force is not conservative, giving rise to a non-equilibrium process, and reads F α (x) = (α − |x| 2 )x + x + where x + = (−x 2 , x 1 ). This term corresponds to a rotation, which is divergencefree with respect to the (Gibbsian) invariant measure and, therefore, does not change the stationary state, but it makes it a non-equilibrium one [109][110][111]. The systematic study of linear response theory for such non-equilibrium systems is an interesting problem that we leave for future study. In the thermodynamic limit, the system is described by a McKean-Vlasov equation whereF = F α − θ x, the last term representing the mean field contribution of the coupling to the local force. The authors in [89] prove that the infinite particle system undergoes a phase transition, with a stationary measure ρ 0 (x) losing stability to a time-dependent probability measureρ = ρ(x, t). Physically, this phenomenon can be interpreted as a process of synchronization. In fact, ρ (0) (x) represents a disordered state, with the oscillators moving out of phase, whileρ describes a state of collective organization with the oscillators moving in an organized rhythmic manner. The transition can be investigated via the order parameter x which vanishes in the asynchronous state, x 0 = 0 , and is different from zero and time dependent in the synchronous state. In particular, the stationary measure ρ 0 (x) can be written as and satisfies the stationary McKean-Vlasov equation with M 0 α,θ, x 0 (g) = −∇ · [Fg] + (σ 2 /2) g, being the Fokker-Planck operator describing the stationary state ρ (0) (x). Note that x 0 = 0. We can perform a linear response theory around this stationary state ρ 0 by replacing F α → F α + εX(x)T(t) and studying the perturbation ρ (1) of the measure defined via ρ(x, t) = ρ (0) (x) + ερ (1) (x, t). As previously outlined, ρ (1) (x, t) satisfies equation (2.4) from which the whole linear response theory follows. However, to conform to the notation in [89] we will here define ρ (1) (x, t) = (ρ (0) ) 1/2 q(x, t) and write the corresponding equation for q(x, t). After some algebra, it is possible to write that ∂ t q(x, t) = M α,θ,0 (q) − T(t)(ρ (0) ) −1/2 ∇ · X(x)ρ (0) + θ(ρ (0) ) 1/2 (ρ (0) ) 1/2 y, q(y, t) · ∇φ(x) =M α,θ,0 (q) − T(t)(ρ (0) ) −1/2 ∇ · X(x)ρ (0) (6.15) where we have defined We mention that operator has the structure of a Schrödinger operator in a magnetic field ( [63], Sec. 4.9). Furthermore,M α,θ,0 (q) = M α,θ,0 (q) + θ(ρ (0) ) 1/2 (ρ (0) ) 1/2 y, q(y, t) · ∇φ(x) with f , g = dyf (y)g(y) being the usual scalar product. In particular, let us observe that (ρ (0) ) 1/2 y, q(y, t) = (ρ (0) ) 1/2 yq(y, t)dy = yρ (1) (y, t)dy = y 1 . A formal solution of the above equation is +θ (ρ (0) ) 1/2 (ρ (0) ) 1/2 y, q(y, s) · ∇φ(x) (6.17) which is the analogous of equation (2.5). Using the above expression, we can evaluate the response of the observable x i as × −T(s)(ρ (0) ) −1/2 ∇ · X(x)ρ (0) + θ(ρ (0) ) 1/2 (ρ (0) ) 1/2 y, q(y, s) · ∇φ(x) . In particular, their spectrum is related to the Fourier transform of the mean field susceptibility Γ α,θ and macroscopic susceptibilityΓ α,θ (respectively) through equations similar to (2.17) and (2.27). The authors in [89] study the spectrum of both these operators in order to perform a stability analysis of the stationary distribution ρ 0 (x). In particular, they observe that the operator M α,θ,0 can be written as M α,θ, and vanishes, cannot happen in this setting. In particular, correlation properties will never diverge. Phase transition can, instead, take place according to the scenario (ii) above. Indeed, the authors in [89] show that the spectral gap of the operatorM α,θ,0 vanishes at surface in the (α, σ , θ) parametric space defined by the following equation: where A = α/θ − 1 and δ = √ 2σ 2 /θ . In particular, they are able to prove that the eigenvalues associated with eigenfunctions ofM α,θ,0 which are orthogonal to the subspace of L 2 (R 2 ) spanned by (ρ (0) ) 1/2 and n · x(ρ (0) ) 1/2 , n ∈ R 2 being any unit vector, are always negative. Nevertheless, M α,θ,0 can become unstable from eigenfunctions which are not orthogonal to n · x(ρ (0) ) 1/2 . It is possible to identify the eigenfunctions that at the transition yield eigenvalues with vanishing real part. In particular, at the transition line (6.21), the eigenfunction Ω(x) = (0, 1) · x(ρ (0) ) 1/2 + i(1, 0) · x(ρ (0) ) 1/2 gives an eigenvalueλ j = i, with Ω(x) * corresponding to the complex conjugate eigenvalueλ * j = −i. The macroscopic susceptibility (2.27) consequently develops a pair of symmetric poles in ω = ±1, corresponding to a dynamic phase transition, giving rise to a Hopflike bifurcation yielding the time-dependent stateρ(x, t) that defines the synchronized state. As a result, near the transition, the order parameter x , where the expectation value is computed using the measureρ(x, t), oscillates at frequency ω = 1 with amplitude A 1 (α, σ , θ). Instead, since it is a quadratic quantity, the rescaled variance θ/σ 2 z 2 , where z = x − x , oscillates at frequency ω = 2 with amplitude A 2 (α, σ , θ) around the value B 2 (α, σ , θ).   (panel (a)), of the amplitude A 2 of the oscillations of the re-scaled variance θ/σ 2 z 2 (panel (b)), and of the time mean value of θ/σ 2 z 2 (panel (c)). The red dotted line represents the transition line given by equation (6.21); see [89]. See text for details. (Online version in colour.) We have investigated this non-equilibrium transition through numerical integration of equation (6.11) via an Euler-Maruyama scheme [108]. The convergence of our results to the thermodynamic limit has been tested by looking at increasing values of the number of agents. We display here the results by taking N = 5000 and choosing α = 2. Figure 3 shows the value of A 1 (α = 2, σ , θ ) (panel a), A 2 (α = 2, σ , θ ) (panel b) and B 2 (α = 2, σ , θ) (panel c) in the parametric region 0.2 ≤ σ ≤ 3, 0.5 ≤ θ ≤ 6 of the two dimensional parameter space (σ , θ). For the sake of clarity, we also provide in figure 4 a snapshot of a horizontal and vertical section of the heat plots.
These numerical experiments confirm that the system indeed undergoes a continuous phase transition, with a collective synchronization stemming from a disordered state as the system passes through the transition line given by equation (6.21) for α = 2. Let us remark again that the fluctuations, being related to the spectrum of M 0,+ α,θ, x 0 , are always finite, see figure 4.

Conclusion
The understanding of how a network of exchangeable interacting systems responds to perturbations is a problem of great relevance in mathematics, natural and social sciences, and  technology. One is in general interested in both the smooth regime of response, where small perturbations result into small changes in the properties of the system, and in the non-smooth regime, which anticipates the occurrence of critical, possibly undesired, changes. Often, critical phenomena, which can be triggered by exogeneous or endogenous processes, are accompanied by the existence of a large-scale restructuring of the system, whereby spatial (i.e. across systems) and temporal correlations are greatly enhanced. The emergence of a specific spatial structure is especially clear when considering order-disorder transitions. Spatial-temporal coordination becomes evident when studying the multi-faceted phenomenon of synchronization. Finally, slow decay of temporal correlations-the so-called slowing down-indicates that nearby critical transitions the negative feedback of the system become ineffective. This paper is the first step in a research programme that aims at developing practical tools for better understanding and predicting-in a data-driven framework-critical transitions in complex systems. We have here developed a fairly general theory of linear response for such a network in the thermodynamic limit of an infinite number of identical interacting systems undergoing deterministic and stochastic forcing. Our approach is able to accommodate both equilibrium and non-equilibrium stationary states, thus going beyond the classical approximation of gradient flows. We remark that the existence of equilibrium stationary (Gibbs) states, the gradient structure (in a suitable metric) and the self-adjointness of the Fokker-Planck operator are equivalent. The presence of interaction between the systems leads to McKean-Vlasov evolution equation for the one-particle density, which reduces to the classical Fokker-Planck equation if the coupling is switched off.
We find explicit expressions for the linear susceptibility and are able to evaluate its asymptotic behaviour, thus allowing for the derivation of a general set of Kramers-Kronig relations and related sum rules. The susceptibility, in close parallel to the classic Clausius-Mossotti expression of macroscopic electric susceptibility for condensed matter, is written in a renormalized form as the product of a matrix describing the self-action of the system times the mean field susceptibility. This allows for further clarifying the relationship between endogenous and exogenous processes, which generalized the fluctuation-dissipation theorem for this class of systems.
Linear response breaks down when the susceptibility diverges, i.e. it develops poles in the real axis. We separate two scenarios of criticality-one associated with the divergence of the mean field susceptibility, and another one associated with singularities of the matrix describing the self-action of the system. The first case pertains to the classical theory of critical transitions.
The second case is here for us of greater interest and can be realized only in the thermodynamic limit. We interpret such a second scenario as describing phase transitions for the system. We define two scenarios of phase transition-a static one, and a dynamic one, where a pole at vanishing frequency and two poles at opposite frequency appear in the linear susceptibility, respectively. At the phase transition, the Kramers-Kronig relations and sum rules valid in the smooth regime of response break down and a detailed study of the poles allows one to find the correction terms. Again, one can establish a link with results from condensed matter physics, as the correction terms resemble those appearing when studying frequency-dependent optical properties of a material at the insulator-metal phase transition, where the static conductivity becomes non-vanishing. We prove that, against intuition, a phase transition is-as opposed to the case of critical transitions-not accompanied by a divergence in the autocorrelation properties of the system, i.e. no critical slowing down is observed. Our interpretation is supported by the use of the formalism developed in this paper to revisit through analytical and numerical tools the classical results for phase transitions occurring in the Desai-Zwanzig model on the Bonilla-Casado-Morrillo model, for which it is easy to define appropriate order parameters. The criticalities in the these two models conform to the scenario of static and dynamic phase transition, respectively.
We remark that studying the linear response of the order parameter is the optimal choice for detecting the phase transition but not the only one. In fact, we expect that a broader class of observables can be used in order to identify the critical behaviour. This is especially important in non-equilibrium cases, where the identification of such order parameter can be extremely nontrivial.
The work reported in this paper opens up several avenues for future research. Four natural next steps are: (a) to investigate in greater detail multi-dimensional reversible (equilibrium) McKean-Vlasov dynamics exhibiting phase transitions; for such systems the self-adjointness of the linearized McKean-Vlasov operator enables the systematic use of tools from spectral theory for selfadjoint operators in appropriate Hilbert spaces. (b) To use the analytical tools developed in this paper to design early warning signals for phase transitions, as opposed to critical transitions for which there exists an extensive literature. (c) To better define the class of observables for which the divergence of the linear response can be used to define and detect phase transitions. In particular, we aim at developing systematic analytical and data-driven methodologies for identifying order parameters in agent based models. These tools will enable us to move beyond the quadratic interaction between subsystems that was considered in this paper. (d) To use the framework developed in this paper in order to revisit phenomena such as synchronization, cooperation and consensus in multi-agent systems, and more generally the emergence of coherent structures in complex systems, both in natural and social sciences as well as technology.
Data accessibility. The codes used to run the simulations and the data used for producing the figures are available on figshare.com at https://figshare.com/projects/Response_theory_and_phase_transition_for_ thermodynamic_limit_of_interacting_identical_systems/89516.