Response Theory Identifies Reaction Coordinates and Explains Critical Phenomena in Noisy Interacting Systems

From an operational point of view, the investigation of these critical phenomena heavily relies on the identification of a suitable set of observables, commonly referred to as order parameters (OPs), collective variables (CVs) or reaction coordinates (RCs), acting as thermodynamic variables by providing a low dimensional description of the macroscopic feature of the system. Whereas it is true that all the three notions of OP, CV and RC are related to a dimension reduction procedure provided by a projection from a high dimensional system to a reduced state space, they refer to different features of the system. Along the lines of [27], we define a CV as any function of the full phase space. An OP is then a suitable combination of CVs that is able to distinguish between macroscopic, possibly metastable, states of the sytem. As some relevant parameters of the system are varied, one wishes to identify critical points by observing asymptotic properties of the OPs, that is by constructing phase diagrams. On the other hand, the term RC refers to a suitable OP that is able to capture time dependent properties of the system, e.g. relaxation phenomena or transitions between metastable states.
A major challenge in the identification of critical phenomena is that, while order parameters or reaction coordinates can in many cases be deduced for equilibrium systems using, e.g., symmetry arguments, the identification of reaction coordinates in nonequilibrium settings is far less trivial [28, 29, 27] and is tailored to the specific problem at hand. Here we provide a constructive, as opposed to ad hoc, criterion to identify a set of reaction coordinates $\{A_{i}\}$ for a large class of interacting agent systems based on the microscopic interaction structure among the agents. Inspired by the success of response theory in explaining critical phenomena in a class of interacting systems [30, 31, 32, 33], we show that, according to this perspective, the $\{A_{i}\}$ not only represent OPs for the system but correspond to suitable RCs as they act as proper nonequilibrium thermodynamic variables. On the one hand, their stationary asymptotic properties identify the different macroscopic states of the system. On the other hand, they also carry information on correlation properties and resilience of the system to perturbations. Most importantly, we show that such observables represent the optimal observables to fully characterise the development of a critical phenomenon as the critical, strongly peaked, resonances in their linear response to external perturbations allow for a quantitative description of the critical behaviour.

Linear Response Theory provides a powerful framework to predict how statistical properties of a system change as a result of weak, yet arbitrary exogenous perturbations. Originally developed for statistical mechanical systems near equilibrium, the theory has much larger validity [34] with applications ranging from driven, strongly dissipative systems such as the climate [35, 36] to neural networks [37, 38, 39] and financial markets [40]. In brief, a linear response perspective corresponds to the following conceptual and mathematical framework. Given a system with a unique, physical, ergodic (not necessarily smooth) measure $\mu_{0}$, linear response theory aims at evaluating changes in the expectation value

where $\Phi$ is a generic observable of the system, $\mathbf{x}(t)$ is a typical trajectory and the last equality hold because of ergodicity. Linear response theory is valid for a vast class of systems ranging from deterministic dynamical systems [34, 41], Markov chain models [42] to continuous stochastic systems [43, 44] where it can be justified on rigorous mathematical grounds [45]. We here consider a finite dimensional system described by a set of Ito stochastic differential equations with deterministic drift $\mathbf{\mathcal{F}}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and noise law given by the matrix valued function $\boldsymbol{\Sigma}:\mathbb{R}^{d}\to\mathbb{R}^{d\times q}$

where $\mathrm{d}\mathbf{W}_{t}$ is a $q-$dimensional Wiener process. If we assume that the stochastic process is hypoelliptic [44], the invariant measure $\mu_{0}$ associated to the above stochastic process is absolutely continuous with respect to the Lebesgue measure with smooth probability density $\rho_{0}$, that is $\mu_{0}(\mathrm{d}\mathbf{y})=\rho_{0}(\mathbf{y})\mathrm{d}\mathbf{y}$, even when the deterministic drift $\mathbf{\mathcal{F}}$ induces a chaotic dynamics. We perturb the system by letting $\mathbf{\mathcal{F}}\to\mathbf{\mathcal{F}}+\varepsilon\mathbf{Y}(\mathbf{y})T(t)$ where the bounded function $T(t)$ represents the time modulation of the external perturbation and $Y(y)$ its phase space state dependence. In this setting, the statistical properties of the system change as

where $G_{\Phi}$ is the Green’s function associated to the observable $\Phi$. The Green’s function $G_{\Phi}$ encodes all the positive and negative feedbacks of the system and allows to evaluate the response of the system to general time dependent forcings. Response formulas provide a useful decomposition of any Green’s function as¹¹1In expression (4) the contribution of the continuous part of the spectrum is neglected, which is equivalent to assuming that the Markov semigroup describing the dynamics is quasi compact [46].

where the complex valued exponential rates $\lambda_{k}\in\mathbb{C}$ and $\mathbf{Re}\lambda_{k}<0$ are the eigenvalues of the generator of the Markov semigroup associated to the unperturbed dynamics (2) and are referred to as stochastic Ruelle Pollicott resonances and $g_{k}^{(l)}$ are suitable coefficients [47, 48, 49]. Moreover, if the eigenvalue $\lambda_{k}$ is not simple, that is its multiplicity is $m_{k}>1$, modulating polynomial terms arise in the expression for the Green’s function. We remark that the exponential rates $\lambda_{k}$ represent an intrinsic property of the unperturbed dynamics and correspond to natural relaxation and oscillatory timescales of the system. Green’s functions are thus the relevant quantities to be investigated in order to detect the onset of critical instabilities of general systems. In the context of interacting systems, it has been shown that long lasting and persistent modes with exponential rates of small real part are associated to poorly damped instabilities due to phase transitions in the thermodynamic limit [31].
It is important to stress that the choice of the observable under investigation remains a critical issue in detecting criticalities as the coefficients $g^{(l)}_{k}$ are not an intrinsic property of the system but depend both on the observable $\Phi$ and the state dependence $\mathbf{Y}(\mathbf{y})$ of the perturbation [47]. In particular, if the observable $\Phi$ has vanishingly small projection on the most unstable mode, its associated Green’s function $G_{\Phi}$ would not represent a good proxy for the detection of the critical behaviour of the system.

The main contributions of this paper can be summarised as follows:

• •

  We provide a constructive way to define a finite set of RCs $\{A_{i}\}$ for the thermodynamic limit of an ensemble of interacting agents. The RCs correspond to an effective dimension reduction of the system in terms of a few relevant observables. We provide mathematical conditions on the interaction structure among the agents for which the construction of RCs is exact. In the context of synchronisation phenomena for ensembles of phase oscillators with purely harmonic coupling, dimension reduction is usually obtained by invoking the Ott-Antonsen ansatz and generalisations thereof [50, 51, 52, 53]. Our results are instead valid for any coupling structure and even in settings where phase reduction techniques of nonlinear oscillators [54, 55] cannot be applied. In this context, this paper represents a generalisation of [6] to settings with (a) possibly multiplicative, white noise sources and (b) an interaction structure that is not a priori assumed to depend on (functions of) the empirical mean of the ensemble.

• •

  We derive response formulas for the thermodynamic limit of the ensemble of agents. We show that the response of general mean field observables of the system is mediated by these RCs. In particular, we show that the coupling among the agents manifest itself in the thermodynamic limit as a memory contribution, conveyed by the RCs, to the sensitivity properties of the system. We prove that the macroscopic Green’s function describing the response of the system at a macroscopic level encodes feedbacks stemming from both microscopic properties of each agent and the interaction structure among them.

• •

  We show that the investigation of criticalities arising from phase transitions are associated, in the thermodynamic limit, to the divergent behaviour of Green’s functions relative to the RCs. Phase transitions can be quantitatively identified as the singular part of the complex valued susceptibility, the Fourier Transform of the Green’s Function. These singularities can be characterised by a pole-residue pair $(\omega^{*},\kappa)$ with $\omega^{*}\in\mathbb{R}$ and $\kappa\in\mathbb{C}$. The pole $\omega^{*}$ corresponds to frequencies of the driving force that would critically destabilise the system, whereas the residue $\kappa$ determines the strength and detailed nature of the singularity.

• •

  We provide analytical and numerical evidence that the analysis of response properties of the RCs in a finite size ensemble of agents is able to capture the signature of a phase transition as a resonance around the singular thermodynamic behaviour. We show that both the pole $\omega^{*}$ and the residue $\kappa$ can be inferred in a finite system. In particular, we illustrate the physical role that the residue $\kappa$ plays on response properties of the system.

The paper is structured as follows. In section 2 we introduce the class of interacting systems under investigation and we describe when a dimension reduction in terms of RCs is possible. Section 3 is dedicated to the description of linear response properties of the system in terms of the RCs. In section 4 we investigate the critical response properties of a finite ensemble of Kuramoto oscillators. In the conclusions we mention future research directions of this linear response framework to investigate resilience features of general interacting systems.

where $\langle\langle\cdot\rangle\rangle$ represents the double expectation value with respect to $\rho(\mathbf{x},t;\mathbf{h})$ and $\mu$. Equation (6) can now be written as

where $\langle\langle\mathbf{A}(\mathbf{x})\rangle\rangle=\{\langle\langle a_{i;l}(\mathbf{x})\rangle\rangle\}^{M,d}_{i=1,l=1}\in\mathbb{R}^{p}$ is a finite vector ($p$ not necessarily equals to $d$ and $p\ll N$, see later discussion) of observables and $A_{j}(\mathbf{x})$ $j=1,\dots,p$ are linearly independent functions. We have also defined the nonlinear operator $\mathcal{L}_{\langle\langle\mathbf{A}\rangle\rangle}$ depending on the probability distribution through the vector $\langle\langle\mathbf{A}\rangle\rangle$. Following the terminology of [60], this scenario corresponds to a finite nonlinearity dimension equals to $p$ for equation (6) and it corresponds to a natural dimension reduction of the ensemble of agents in terms of the variables $\langle\langle\mathbf{A}\rangle\rangle$.
Separable kernels $\mathbf{K}(\mathbf{x},\mathbf{y})=\nabla\mathcal{U}(\mathbf{x}-\mathbf{y})$ arising from interaction potentials are routinely used in applications, with radial potentials of the form $\mathcal{U}(\mathbf{x})=u(|\mathbf{x}|)$ being a common choice. A paradigmatic example is represented by quadratic interactions $u(x)=\frac{x^{2}}{2}$, yielding a separable kernel $\mathbf{K}(\mathbf{x},\mathbf{y})=\mathbf{x}-\mathbf{y}$. This corresponds to linear forces among the particles trying to synchronise them towards their centre of mass $\bar{x}=\frac{1}{N}\sum_{j}{\mathbf{x}}_{j}(t)$ and has numerous applications, see our previous work [30, 31] on this. In this case, $\forall i=1,\dots,M$, the expansion (7) of the interaction kernel holds with $d=2$ and $(f_{i;1},g_{i;1},f_{i;2},g_{i;2})=(x,0,0,-y)$ resulting in a non linear dimension $p=M$ for equation (9) given by the observables $\langle\langle\mathbf{A}(\mathbf{x})\rangle\rangle=\langle\langle\mathbf{x}\rangle\rangle$. Similarly, higher order polynomials $u(x)=\sum_{k=2}^{n}u_{k}x^{k}$ with $n\geq 2$ give rise to separable kernels and a finite nonlinearity dimension. As a further example, we mention separable kernels provided by trigonometric interactions of any order with applications to synchronisation phenomena, see [63, 8] and the section below on the Kuramoto model. We refer the reader to the conclusions section for a discussion regarding non separable kernels.
We remark that equation (9) can be interpreted as a thermodynamic formulation of the system with the $\langle\langle\mathbf{A}\rangle\rangle$s being a generalisation of thermodynamic variables to nonequilibrium settings [64]. Firstly, stationary distributions $\rho_{0}=\rho_{0}(\mathbf{x};\mathbf{h})$ of (9) can be parametrised in terms of the stationary values of the thermodynamic variables $\langle\langle\mathbf{A}\rangle\rangle_{0}$ according to the self consistency equation given by the stationary condition $\mathcal{L}_{0}\rho_{0}=0$ where $\mathcal{L}_{0}=\mathcal{L}_{\langle\langle\mathbf{A}\rangle\rangle_{0}}$ is the operator defined in (9) evaluated at stationarity. This is better understood in the case of homogeneous equilibrium systems, characterised by a gradient dynamics $\mathbf{F}_{\alpha}(\mathbf{x})=-\nabla V_{\alpha}(\mathbf{x})$, an interaction potential $\mathcal{U}(\mathbf{x})$, thermal noise $\mathbf{s}(\mathbf{x})=\mathbf{1}_{M\times M}$ and no quenched disorder. In such systems, stationary solutions of (6) satisfy the self consistency equation [20]

where $Z$ is a normalisation constant. Assuming a separable kernel, equilibrium stationary solutions are parametrised by a finite number of observables $\langle\mathbf{A}\rangle_{eq}$ since they can be written as $\rho_{eq}\left(\mathbf{x};\langle\mathbf{A}\rangle_{eq}\rangle\right)=\frac{1}{Z}\exp\left(-\frac{2}{\sigma^{2}}\left(V_{\alpha}(\mathbf{x})-f(\mathbf{x},\langle\mathbf{A}\rangle_{eq})\right)\right)$ where $f$ is a suitable function deriving by the convolution product $\mathcal{U}\star\rho_{eq}$, similarly to equation (8).
Moreover, the $\langle\langle\mathbf{A}(\mathbf{x})\rangle\rangle$s represent proper RCs of the system as they govern time dependent properties of the ensemble of agents, such as its response, correlation and memory features, as we elucidate in the next section.

Following [43], we evaluate response properties of the finite system by performing $n$ simulations of equations (20) where the initial conditions are sampled according to $\rho_{0}=\frac{1}{2\pi}$ and where at time $t=0$ we apply a macroscopic perturbation $F(x;h)=h\to F(x;h)=h+\varepsilon X(x)\delta(t)$. The macroscopic Green’s function $\tilde{G}_{1}(t)$ is obtained as the average over the $n$ ensemble members of $\langle\langle A_{1}\rangle\rangle_{1}(t)$ which we estimate as $\frac{1}{N}\sum_{i=1}^{N}A_{1}(x_{i}(t))$. We remark that, since the unperturbed measure is uniform, a constant perturbation $X(x)=1$ would result in a vanishing response as $G_{1}(t)\equiv 0$, see (15a). We have here chosen a non trivial state dependent perturbation $X(x)\propto\sin(x)$. We note that such term has been used to model excitability features in Kuramoto-like models for brain dynamics [70]. The numerical experiments indicate that the finite system susceptibility can be written, at the phase transition, as

where $\psi(\omega)$ is an analytic function in the upper complex $\omega$ plane and the residue $\kappa$ characterises the strength of the divergence. The residue is related to the properties of $P(\omega)$ and can be found as $\kappa=\left(\theta P^{\prime}(\omega^{*})\right)^{-1}$, see Appendix for more details.

We remark that the finiteness of the system results in a mollifying effect of the singularity into a resonant behaviour, where the poles are slightly shifted as $\pm\omega^{*}\to\pm\omega^{*}+i\gamma(N)$, with $\gamma(N)\to 0$ as $N\to+\infty$. In the homogeneous case the residue turns out to be completely imaginary and equal to $\kappa=\frac{i}{2}$, whereas for the bimodal case it results in $\kappa=\kappa_{R}+i\kappa_{I}=-\frac{D}{4\omega^{*}}+\frac{i}{4}$. In the former case, one expects that, as $N\to+\infty$, the function $\omega\mathbf{Im}\tilde{\chi}_{1}(\omega)$ will be constant and equal to $|\kappa|=\frac{1}{2}$ in the proximity of the pole and to quickly drop to zero for $\omega=0$. A visual inspection of the bottom inset of panel (a) of Figure 1 provides a clear evidence of such theoretical prediction.
A fundamentally different situation occurs in the latter case as the residue has a non-vanishing real part leading to a distortion of the resonance at $\omega^{*}$ as $\mathbf{Im}\tilde{\chi}(\omega^{*})\neq 0$ and, thus, affecting the way the amplitude of the oscillations of the response track the forcing signal. This is more easily understood when considering a sinusoidal driving forcing $T(t)=\sin(\hat{\omega}t)$. The response (in the time domain) can be obtained by performing an inverse Fourier transform of (18) and results in

The amplitude of the response is given by the amplitude of the macroscopic susceptibility. We expect that driving frequencies around the pole $\omega^{*}$ would destabilise the system resulting in greatly amplified oscillations, whereas the system would be much more resilient for higher driving frequencies. Moreover, if the imaginary part of the susceptibility were to vanish, the response would feature either no phase shift as $\phi=0$ or phase reversal as $\phi=-\pi$. At the critical frequency $\hat{\omega}=\omega^{*}$ we have, from (26), that the response is given by $\langle\langle A_{1}\rangle\rangle_{1}(t;\omega^{*})=\frac{|\kappa|}{\gamma(N)}\sin\left(\omega^{*}t+\phi\right)$ with $\phi=\arctan(-\frac{\kappa_{I}}{\kappa_{R}})-\frac{\pi}{2}=\arctan(\frac{\omega^{*}}{D})-\frac{\pi}{2}$. The non vanishing real part of the residue yields a spontaneous phase shift of the response to perturbations as confirmed by Figure 2 (left panels). We remark that this provides a quantitative estimate, at the phase transition, of the phase shift phenomenon observed (in deterministic settings) in [71]. On the other hand, if the system is driven at higher frequency, e.g. $\hat{\omega}=2\omega^{*}$, we expect the real part of the susceptibility to decay faster than the imaginary part [72] and the response to be in quadrature with the forcing, i.e. $\langle\langle A_{1}\rangle\rangle_{1}(t;2\omega^{*})\propto\cos(2\omega^{*}t)$, as confirmed by Figure 2 (right panels).

In this paper we have considered the relevant issue of the identification of a set of reaction coordinates $\{A_{i}\}$ for the thermodynamic limit of ensembles of interacting agents in the general case of pairwise interaction and quenched disorder in the dynamics. Given a general class of interaction structure among the agents, we have provided a framework to define the $\{A_{i}\}s$ and construct a nonequilibrium thermodynamical formulation of the macroscopic description of the system. These reaction coordinates provide a parametrization of the stationary measure and fully characterise the linear response of the system to perturbations. Such a model reduction strategy is particularly useful as it allows to study not only smooth regimes of the response but also critical behaviours, characterised as an emergence of singular poles of the complex valued susceptibility of the system. The analysis in terms of the reaction coordinates of the paradigmatic Kuramoto model allows us to show (a) that our framework encompasses known results on the stability properties of the incoherent state and (b) that the time dependent properties of the onset of synchronisation can be quantitatively investigated through the susceptibility function describing the response of the system to general driving forces.
Our results are exact in the case of separable interactions. However, non separable kernels are also found in applications. If $\mathbf{K}(\mathbf{x},\mathbf{y})$ is a smooth function, an expansion as in (8) still holds by using a Schauder decomposition [73] but it contains in general an infinite number of terms, resulting in an infinite number of observables $\{A_{i}\}$ with a clear loss of dimension reduction of the system. However, there exist numerous rigorous techniques to approximate in a controlled manner a general kernel $\mathbf{K}(\mathbf{x},\mathbf{y})$ in terms of suitable separable kernels, see [62, Ch. 11][74]. Such approximating methods would be the basis for the development of phenomenological theories of interacting systems with suitable separable interactions for which the identification of a finite number of observables can be performed exactly.
On the one hand, the investigation of mean field models is very relevant as it provides a way of gaining a mathematical and analytical understanding of the smooth and critical features of ensembles of interacting agents. On the other hand, the lack of an underlying geometrical structure might pose a too restrictive limitation to the applications of these models. Future research in this direction will involve the generalisation of this framework to systems with a microscopic network structure [75] and in particular to all systems amenable to a description in terms of graphon [76, 77] and graphop theory [78, 79]. We envision that this framework can be extended to include the physically relevant cases where the dynamics is not assumed to be overdamped, in the direction of understanding nonequilibrium ensembles in statistical mechanics [80] and to non-Markovian interacting particles [81], where appropriate Markovian closures have to be considered [44, 82].