Spectroscopy of phase transitions for multiagent systems

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 behaviour of the susceptibility. We find clear evidence of the loss of analyticity due to a pole crossing the real axis of frequencies. Such behaviour has a degree of universality, as it does not depend on either the applied forcing nor 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.


I. INTRODUCTION
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 , synchronisation 4 , systemic risk 5 and consensus formation 6,7 . They are fundamental in developing algorithms for sampling and optimization 8 and they have also been used for the management of natural hazard 9 and climate change impact 10 .
Multiagent systems can often exhibit abrupt changes in their behaviour, 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 climate change, market crashes etc 11,12 . The importance of developing tools for predicting critical transitions has long been recognized. One of a) Electronic mail: n.zagli18@imperial.ac.uk 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][15][16] . Linear response theory has been successfully applied to classic problems of solid-state physics and optics 17 as well as plasma physics and stellar dynamics 18 (Ch.5); see some examples of application of the theory in both equilibrium and nonequilibrium systems [19][20][21][22][23][24] . Rigorous mathematical foundations for linear response theory have been provided for the case of Axiom A systems 25,26 (see e.g. 27 for further developments in the context of deterministic systems) and for diffusion processes, both in finite and in infinite dimensions 28,29 ; see also the interesting contributions 30 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 zero [31][32][33][34] as the Ruelle-Pollicott poles 35,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, multiagent system can also undergo a qualitative change of their properties through a different mathematical mechanism, namely phase transitions 3,37 , defined as exchange of stability of nonunique stationary distributions as the parameters of the systems vary; see a detailed analysis in 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 renormalised 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 emphasised in 12,43 . The link between phase transitions and slow decay of correlations for interacting particle systems is well established, see. e.g. 44 .
A. This Paper In this paper we focus in much greater detail on the relationship between the occurrence of phase transitions and the nonanalyticity 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 part 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, an 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 synchronisation transition of an infinite collection of coupled nonlinear oscillators. As the transition point is crossed, the order parameter (magnetization) acquires a non vanishing constant value for the DZ model and is periodically oscillating for the BCM model.

II. THE GENERAL FRAMEWORK
We investigate a system composed of N exchangeable interacting M dimensional sub-systems whose dynamics is determined by the following Itô stochastic differential equations (SDEs) (1) where x k ∈ R M and k = 1, . . . , N. F α : R M → R M is a smooth vector field, possibly depending on a parameter α, and W denotes a standard P−dimensional Brownian motion; S : R M → R M×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| 2 2 . 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 synchronising 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 equation 3,47-49 where D = SS T . 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 Eqn. 2 correspond to local minima of V α (y). In this case, the thermodynamic limit (2) can be written in the standard form as where Z N denotes the partition function of the N−particle system and F(ρ) denotes the free energy of mean field system 50 . A stationary state is characterised by the order parameter x 0 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 40 , the response of the order parameter in the frequency domain is written in terms of a macroscopic (or renormalised) susceptibilityΓ i (ω) as (Repeated indices are summed): where P i j (ω) = δ i j − θ Y i j (ω) and the susceptibilities Γ i (ω)and Y i j (ω) are respectively the Fourier Transform of the microscopic response functions that can be written as correlation functions in the unperturbed state as 40 represents the adjoint of L x 0 and · 0 is the expectation value on the unperturbed state ρ 0 (x). The Fokker-Planck operator L x 0 appears on the right hand side of (2) (evaluated at the stationary state ρ 0 (x)) and its adjoint L + x 0 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 L + The renormalisation of the susceptibility derives from the coupling among the subsystems; note thatΓ i (ω) inherits the poles of both Γ i (ω) and of the matrix P −1 i j (ω). Away from criticality, both the microscopic and macroscopic susceptibilities are analytic in the upper half of the ω complex plane.
As discussed in 40 , the critical behaviour of this class of multiagent systems, signified by the singular behaviour of the susceptibility, originates from two distinct physical phenomena that are associated to either the poles of Γ i (ω) or P −1 i j (ω). 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 P −1 i j (ω) 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 L + x 0 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) nor on the choice of the observable (x, in our case) and can be related to the spectral properties of a modified transfer operator 40 .

III. NUMERICAL RESULTS
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 (R 2 ). 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 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 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 estimateG i (τ; N) of the renormalised response function . Details on the numerical simulations are also reported in the Appendix. Figure 1 shows the response functionsG 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 towards the unperturbed state, while the one for the BCM combines the decay with an oscillatory behaviour taking place at the natural frequencyω = 1.
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 lim N→∞ κ ω−ω 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 38 , one might conjecture that for the Desai-Zwanzig model and related models the function γ would scale as 1/N. Instead, we have observed here that the behaviour 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 equibrium systems, the mean field limit N → ∞ and the limit T → T c , where T c 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 nor 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 behaviour 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 Figure 2 shows the function (ω − ω 0 )Im{Γ i (ω)}. As N → ∞, this function converges to k everywhere except for ω = ω 0 . An explicit expression for k is known 40 in the case of the DZ model 56 Using the statistics of the unperturbed runs we obtain k ≈ 0.89, which agrees within 2% with the one obtained from the limiting behaviour of the susceptibility, thus validating our results. In the case of the BCM model, our procedure allows one to derive a direct estimate k ≈ 0.44; in this case no expression for the residue is available in the literature and, following 40,46 , 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 dynam-  ical latent heat, that is observed for perturbations occurring a non-vanishing frequency. As discussed earlier, the singular behaviour 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 an additional examples.

IV. CONCLUSIONS
The study of how a large network of identical agents respond 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 behaviour, 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 behaviour 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.

DATA AVAILABILITY
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, reference number 57 .

ACKNOWLEDGMENTS
VL acknowledges the support received by the European Union's Horizon 2020 program through the project TiPES (Grant Agreement No. 820970) and from the EPSRC project EP/T018178/1. The work of GP was partially funded by the EPSRC, grant number 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. NZ has been supported by an EPSRC studentship as part of the Centre for Doctoral Training in Mathematics of Planet Earth (grant number EP/L016613/1).

Appendix A: The models
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 non conservative, leading to irreversible and dissipative processes that can exhibit complex behaviours, such as deterministic chaos, see Figure 3. An all-to-all coupling between the subsystems is given by a matrix L i j = ∇U (x i − x j ) where U (x) = U (−x) represents the interaction potential and x i is the state vector of the i-th subsystem. Weakly interacting diffusions are characterised by a coupling strength which 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 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 onedimensional subsystems x k whose dynamics is prescribed by the following equations (see main text for notation) where k = 1, . . . , N and the confining potential V α (x) = − α 2 x 2 + x 4 4 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 show 3 that the critical line is given by 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 selfsynchronisation transition. The time evolution of the network of oscillators is given by the following equations where the local force is not conservative, giving rise to the non equilibrium features of the system, and reads F α (x) = α − |x| 2 x + x + where x + = (−x 2 , x 1 ). The latter term corresponds to a rotation and makes the stationary state a non equilibrium one. The parameter α > 0 controls the amplitude of the oscillations of the individual non linear oscillators. In fact, when θ = σ = 0, each subsystem oscillates as x j (t) = √ α (cos(t + β j ), sin(t + β j )) where β j = tan(x j 2 (0)/x j 1 (0)). The coupling tries to synchronise the subsystems by attract- ing them towards the center of mass 1 N ∑ N j=1 x j . 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 non synchronised state where the order parameter vanishes. In particular, the transition happen at the surface of the (α, θ , σ ) parametric space defined by 46 where A = α θ − 1 and δ = √ 2σ 2 θ . In the following we have θ = α = 2. The colour code for the figures is given by: • non critical black : DZ σ ≈ 1 , BCM σ ≈ 2.

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 a Dirac's δ (τ = 0). The average of the response for the observable x over the n simulations gives an estimate ofG i (τ; N). The response functions away from the transitions are estimated on an ensemble of n = 10 5 simulations, while the critical response functions with n = 10 6 for DZ and n = 7 × 10 6 for BCM. Furthermore we investigate the response up to time τ = 5 × 10 3 . The corresponding sus-ceptibilityΓ i (ω; N) is simply defined as the Fourier Transform ofG i (τ; N). In the main text, we show the results for an additive perturbation X(x) ≡ X. However, the critical behaviour 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) = x 2 , see Figure 4. The response functionG 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 particle 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 behaviour and no evident dependence on N. At the phase transition, the susceptibility develops the expected singular behaviour κ ω−ω 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).