May-Wigner transition in large random dynamical systems

TL;DR

This paper investigates the stability of large, complex, nonlinear random dynamical systems, revealing a phase transition similar to the May-Wigner transition in linear models, using stochastic analysis of Lyapunov exponents.

Contribution

It introduces a stochastic framework for analyzing stability in high-dimensional nonlinear systems, extending the May-Wigner transition concept beyond linear models.

Findings

Identification of a stability-complexity phase transition in nonlinear systems

Derivation of a stochastic model for finite-time Lyapunov exponents

Comparison with existing linear models and known results

Abstract

We consider stability in a class of random non-linear dynamical systems characterised by a relaxation rate together with a Gaussian random vector field which is white-in-time and spatial homogeneous and isotropic. We will show that in the limit of large dimension there is a stability-complexity phase transition analogue to the so-called May-Wigner transition known from linear models. Our approach uses an explicit derivation of a stochastic description of the finite-time Lyapunov exponents. These exponents are given as a system of coupled Brownian motions with hyperbolic repulsion called geometric Dyson Brownian motions. We compare our results with known models from the literature.

Equations

$$\frac{d\mathbf{x}(t)}{dt}=\mathbf{v}(\mathbf{x}(t);t),$$

$$\frac{d\mathbf{r}(t)}{dt}=-\mu \mathbf{r}(t)+\mathbf{M}\mathbf{r}(t)$$

$$\frac{d\mathbf{x}(t)}{dt}=-\mu \mathbf{x}(t)+\mathbf{f}(\mathbf{x}(t)),$$

$$\frac{d\mathbf{r}(t)}{dt}=-\mu \mathbf{r}(t)+\mathbf{M}(t)\mathbf{r}(t),$$

$$\mathbf{a}=\left[\begin{array}{c}{a}_1\\ ⋮\\ a_n\end{array}\right]\text{and}\mathbf{A}=\left[\begin{array}{ccc}{A}_{11}& \cdots & A_{1n}\\ ⋮& & ⋮\\ A_{n1}& \cdots & A_{nn}\end{array}\right],$$

$$\frac{d\mathbf{x}(t)}{dt}=-\mu \mathbf{x}(t)+\mathbf{f}(\mathbf{x}(t);t)$$

$$\mathbf{x}(t)-\mathbf{x}(s)=-\mu {\int }_s^t\mathbf{x}(u)du+{\int }_s^t\mathbf{f}(\mathbf{x}(u);u)du.$$

$$\mathbb{E}[f_i(\mathbf{x};t)f_j(\mathbf{y},s)]=D_{ij}(\mathbf{x},\mathbf{y})\delta (t-s),i,j=1,\dots ,n,$$

$$\mathbf{f}(\mathbf{x};t)\stackrel{\text{law}}{=}\mathbf{f}(\mathbf{x}+\mathbf{a};t)\text{and}\mathbf{f}(\mathbf{Ux};t)\stackrel{\text{law}}{=}\mathbf{Uf}(\mathbf{x};t)$$

$$D_{ij}(\mathbf{r})\stackrel{{\scriptstyle\text{def}}}{{=}}D_{ij}(\mathbf{x},\mathbf{x}+\mathbf{r})=\Gamma_{1}\Big{(}\frac{\lVert\mathbf{r}\rVert^{2}}{2}\Big{)}\delta_{ij}+\Gamma_{2}\Big{(}\frac{\lVert\mathbf{r}\rVert^{2}}{2}\Big{)}r_{i}r_{j}$$

$${\Gamma }_1^{\prime }(0)\le 0\le {\Gamma }_1(0)\text{and}\frac{{\Gamma }_1^{\prime }(0)}{n+1}\le -{\Gamma }_2(0)\le -{\Gamma }_1^{\prime }(0).$$

$$\mathbf{r}(t)-\mathbf{r}(s)=-\mu {\int }_s^t\mathbf{r}(u)du+{\int }_s^t\mathbf{f}(\mathbf{y}(u);u)du-{\int }_s^t\mathbf{f}(\mathbf{x}(u);u)du.$$

$${\int }_s^t\mathbf{f}(\mathbf{y}(u);u)du={\int }_s^t\mathbf{f}(\mathbf{x}(u);u)du+{\int }_s^t\mathbf{J}(\mathbf{x}(u);u)\mathbf{r}(u)du+\text{error},$$

$$\mathbf{r}(t)-\mathbf{r}(s)=-\mu {\int }_s^t\mathbf{r}(u)du+{\int }_s^t\mathbf{J}(\mathbf{x}(u);u)\mathbf{r}(u)du.$$

$$d\mathbf{r}(t)=\big{(}-\mu\,dt\,\mathbf{I}+d\mathbf{B}(t)\big{)}\circ\mathbf{r}(t),\qquad\mathbf{B}(t)\stackrel{{\scriptstyle\text{def}}}{{=}}\int_{0}^{t}\mathbf{J}(\mathbf{x}(u);u)du.$$

$$\mathbb{E}[B_{ik}(t)B_{j\ell }(s)]$$

$$\displaystyle=-\big{(}\Gamma_{1}^{\prime}(0)\delta_{ij}\delta_{k\ell}+\Gamma_{2}(0)(\delta_{ik}\delta_{j\ell}+\delta_{i\ell}\delta_{jk})\big{)}\min\{t,s\}.$$

$${\Gamma }_1^{\prime }(0)=-{\sigma }^2\text{and}{\Gamma }_2(0)=-{\sigma }^2\tau .$$

$${\sigma }^2\ge 0\text{and}-\frac{1}{n+1}\le \tau \le +1$$

$$\mathbb{E}\bigg{[}\bigg{(}\sum_{i=1}^{n}B_{ii}(t)\bigg{)}^{\!2\,}\bigg{]}\geq 0\qquad\text{and}\qquad\mathbb{E}\big{[}\big{(}B_{ij}(t)\pm B_{ji}(t)\big{)}^{2}\big{]}\geq 0\quad\text{for}\ i\neq j.$$

$$d\mathbf{U}(t)=\big{(}-\mu\,dt\,\mathbf{I}+d\mathbf{B}(t)\big{)}\circ\mathbf{U}(t),\qquad\mathbf{U}(0)=\mathbf{I}.$$

$$\underset{t\to \infty }{\lim }(\mathbf{U}(t)\mathbf{U}(t{)}^T{)}^{1/2t}=\mathbf{H}$$

$$\Vert \mathbf{r}(t)\Vert \sim e^{{\mu }_{1}t}\Vert \mathbf{r}(0)\Vert \text{for}t\to \infty .$$

$${\lambda }_i(t)/t\to {\mu }_i\text{for}t\to \infty .$$

$$d{\lambda }_i(t)=(1+2\tau {)}^{1/2}\sigma d{B}_i(t)-\mu dt+\frac{(1+\tau ){\sigma }^2}{2}\sum _{j=1,j\ne i}^n\frac{dt}{\tanh ({\lambda }_i(t)-{\lambda }_j(t))},i=1,\dots ,n.$$

$$\frac{\partial\rho(\lambda;t)}{\partial t}=\sum_{i=1}^{n}\bigg{(}\mu\frac{\partial\rho(\lambda;t)}{\partial\lambda_{i}}-\frac{(1+\tau)\sigma^{2}}{2}\sum_{j=1,j\neq i}^{n}\frac{\partial}{\partial\lambda_{i}}\frac{\rho(\lambda;t)}{\tanh(\lambda_{i}-\lambda_{j})}+\frac{(1+2\tau)\sigma^{2}}{2}\frac{\partial^{2}\rho(\lambda;t)}{\partial\lambda_{i}^{2}}\bigg{)}$$

$$|{\lambda }_i(t)-{\lambda }_j(t)|\sim |{\mu }_i-{\mu }_j|t\text{for}t\to \infty \text{and}i\ne j.$$

$${\lambda }_i(t)={\lambda }_i(s)+(1+2\tau {)}^{1/2}\sigma (B_i(t)-B_i(s))-\mu (t-s)+\frac{(1+\tau ){\sigma }^2}{2}\sum _{j=1,j\ne i}^n{\int }_s^t\frac{du}{\tanh ({\lambda }_i(u)-{\lambda }_j(u))}.$$

$${\int }_s^t\frac{du}{\tanh ({\lambda }_i(u)-{\lambda }_j(u))}\sim \{\begin{array}{ll}t-s& \text{for}i<j\\ s-t& \text{for}i>j\end{array}$$

$${\lambda }_i(t)\sim {\lambda }_i(s)+(1+2\tau {)}^{1/2}\sigma (B_i(t)-B_i(s))+{\mu }_i(t-s)$$

$${\mu }_i=\frac{(1+\tau ){\sigma }^2(n-2i+1)}{2}-\mu .$$

Full text

Introduction

In many complex systems, time evolution is modelled mathematically by a system of $n$ coupled non-linear first-order ordinary differential equations. In vector notation, we may write such a system as

[TABLE]

where $\mathbf{x}$ is an $n$-dimensional vector-valued function and $\mathbf{v}$ is an $n$-dimensional non-linear vector field specified by the model in question. As a mental image, we might think of $\mathbf{x}$ as the location of a test particle moving in a velocity field $\mathbf{v}$.

While a few properties of the dynamical system (1) can be answered within a fairly general context, most questions of practical interest are highly dependent on the choice of the vector field $\mathbf{v}$. This dependence appears to be somewhat problematic because many (if not most) large and complex systems seem inaccessible (even numerically) and, even worse, our knowledge about the vector field $\mathbf{v}$ is often incomplete. So how can we make in any practical predictions whatsoever? Inspired by statistical physics, it seems natural to take the vector field $\mathbf{v}$ to be ‘maximally’ random under some physical constraints. The hope is that universal behaviour is present in the ‘thermodynamical’ limit $n\to\infty$.

A celebrated and highly influential model was introduced in a 1972 paper by May [20]. He imagined an autonomous dynamical system expanded to first order about an equilibrium point $\mathbf{x}_{e}$ and wrote

[TABLE]

with $\mathbf{r}(t)=\mathbf{x}(t)-\mathbf{x}_{e}$ denoting a small perturbation, $\mu$ a positive constant and $\mathbf{M}$ an $n$-by-$n$ random matrix. This linear equation is, of course, much simpler than the original non-linear monstrosity (1). In (2), we should think of $\mu$ as a decay rate such that if the matrix $\mathbf{M}$ vanishes, then any small perturbation dies out at that rate. The matrix $\mathbf{M}$ incorporates local structure near the equilibrium and should be taken random up to some physical constraints. May was mainly interested in the stability of large ecosystems with $\mathbf{M}$ being a community matrix incorporating interactions between species (whether cooperative, competitive, or predator-prey interaction); we refer to [2] for a review of May inspired models in ecology. For many physicists, it might seem more natural to impose a symmetry constraint such as isotropy (in distribution). If we additionally assume Gaussian disorder, then this symmetry constraint would imply that $\mathbf{M}$ is a matrix from so-called Gaussian elliptic ensemble [26].

It is well-known that the linear model (2) describes a stable (unstable) equilibrium point if the largest Lyapunov exponent is less (greater) than zero, or equivalently, if the eigenvalue of $\mathbf{M}$ with the greatest real part is less (greater) than $\mu$. The crucial observation made by May [20] was that (under quite general conditions) this stability-complexity transition turns into a phase transition in the limit of large systems, $n\to\infty$. This phase transition is known as the May–Wigner transition. The simplest situation arises when $\mathbf{M}$ is a symmetric matrix in which case the eigenvalues of $\mathbf{M}$ is identical to the Lyapunov exponents. Here, large deviation results for symmetric Gaussian matrices has been used to show that the May–Wigner transition is a third order phase transition, see e.g. [19].

The main drawback of May’s model (2) is that it is only physically meaningful locally. We can use May’s model to determine whether a given equilibrium point is stable, but the model can never tell us whether this equilibrium is globally or locally stable, neither can it tell us the number of equilibria. This is a significant problem, because large complex systems are notorious for having a plethora of metastable states. In order to address this issue, Fyodorov and Khoruzhenko [13] have recently suggested a non-linear generalisation of May’s model which is globally meaningful. The Fyodorov–Khoruzhenko model reads

[TABLE]

where $\mu$ is a positive constant (as in May’s model) and $\mathbf{f}$ is an $n$-dimensional zero-mean Gaussian vector field assumed to be homogeneous and isotropic (in distribution). It was shown in [13] that this model also has a phase transition in the limit of large systems, $n\to\infty$ (actually there is an additional phase transition [6] but this additional transition is irrelevant for our purposes). If the fluctuations of the Gaussian field are sufficiently small then the system (3) is expected to have a single equilibrium point (this equilibrium must be stable since $\mu>0$); on the other hand, if the fluctuations of the Gaussian field are too large then the system is expected to have exponentially many equilibria. The region of phase space with a single equilibrium is interpreted as a stable phase, while the region with exponentially many equilibria is interpreted as a chaotic (or unstable) phase. Thus, the phase transition between these two regions of phase space is considered as a non-linear generalisation of the May–Wigner transition. We note that the Fyodorov–Khoruzhenko model is highly inspired by previous studies of random landscapes [11, 14] (see [12] for a review) and it is intimately related to (mean field) spin glass theory [21, 7, 3].

In both May’s model (2) and the Fyodorov–Khoruzhenko model (3) it is assumed that the dynamical system is autonomous, i.e. the driving vector field has no explicit dependence on time, $\mathbf{v}(\mathbf{x};t)=\mathbf{v}(\mathbf{x})$. This assumption is reasonable in some models but not in others. Moreover, even in models where the assumption of time-independence is reasonable, this is often only true in an approximate sense. For example, in ecosystems we would expect that the inter-species interactions change over time due to environmental changes. Likewise, in glassy systems we would expect a time-dependence to be introduced by atomic fluctuations. Thus, it seems more than reasonable to attempt to study dynamical systems (1) with an explicit time-dependence and, in particularly, provide non-autonomous generalisations of May’s model and the Fyodorov–Khoruzhenko model. However, it is extremely challenging to include time correlations in a general setting. Thus, some simplification is needed also in the non-autonomous setting. A natural starting point for an analytic treatment is to consider systems driven by a vector field which is white-in-time. From a physical perspective, the assumption that our system is autonomous says that temporal fluctuations are slow compared to spatial fluctuations, while the white-in-time assumption says that spatial fluctuations are slow compared to temporal fluctuations.

It is, of course, not a new idea to consider dynamical systems (1) in which the vector field $\mathbf{v}$ is white-in-time. In fact, this type of dynamical systems plays a central role in the study of fully developed turbulence [22, 10], which was the inspiration for the mathematical theory of stochastic flows [17]. Surprisingly, these ideas have only very recently [15] been applied to the study of May–Wigner-like transitions. In [15], we considered a white-in-time version of May’s model (2),

[TABLE]

where $\mathbf{M}(t)$ is an isotropic Gaussian white-in-time process. It was shown that it was possible to switch from the matrix formulation (4) to a stochastic description of the finite-time Lyapunov exponents. This latter description was used to study stability. Interestingly, it was observed that the stability-complexity phase transition happens at a scale proportional to $n$, while the phase transition in May’s model (2) happens at a scale proportional to $n^{1/2}$. This is a manifestation of the (not so surprising) fact that the non-autonomous model (4) has a lower onset of complexity, than May’s autonomous model (2).

Naturally, the linear model (4) has the same drawback as May’s model, namely it is only physically meaningful in a neighbourhood of a stable trajectory. For this reason, it is the purpose of this paper to introduce and study stability in a white-in-time version of the Fyodorov–Khoruzhenko model which is globally meaningful. We will see that this new model has certain similarities with the linear model (4) but also some important differences. More precisely, we obtain a stochastic description for the finite-time Lyapunov exponents, which turns out to be identical to the linear model from [15] but with a restriction on parameter space. We use this description to draw a stability-complexity phase diagram. On a more heuristic level, we will also consider models with a more general time-dependence and discuss the possibility of persistent time-correlations. The rest of this paper is organised as follows:

— Section 2: we define the model which is the main subject of this paper.

— Section 3: we show how to obtain a stochastic description of the finite-time Lyapunov exponents.

— Section 4: we study stability and observe that special behaviour arises in the weak-noise limit.

— Section 5: we discuss the effect of time-dependence and auto-correlations in a more general setting.

— Section 6: we summarise our conclusions and state some open problems.

— Appendix A: we present some important results for geometric Dyson Brownian motions.

Notation**.**

Throughout this paper, lower case boldface symbols $\mathbf{a},\ldots$ refer to $n$-dimensional (column) vectors, upper case boldface symbols $\mathbf{A},\ldots$ refer to $n$-by-$n$ matrices, and $\mathbf{a}\cdot\mathbf{b}$ refers to the usual dot product. Moreover, we will always use the standard notation,

[TABLE]

for entries of vectors and matrices, and $\lVert\mathbf{a}\rVert=(\mathbf{a}\cdot\mathbf{a})^{1/2}$ for the (Euclidean) norm.

A model for large complex dynamical systems

Let $\mu>0$ be a positive constant, $\mathbf{x}$ be an $n$-dimensional vector-valued function dependent on an external time $t$, and $\mathbf{f}$ be an $n$-dimensional vector field depending explicitly on the time $t$. Consider a generic dynamical system

[TABLE]

or equivalently, in integral form,

[TABLE]

In this paper, we consider the case where $\mathbf{f}$ is a zero-mean Gaussian random vector field, which is smooth-in-space and white-in-time. We note that this system is not independent of stochastic regularisation, i.e. there will be an ‘Itô–Stratonovich dilemma’. In the following, we have chosen the Stratonovich interpretation, which corresponds to the assumption that the underlying physical phenomenon responsible for the noise is time-reversible. We emphasize that all computations in this paper can easily be performed using any other stochastic regularisation (e.g. Itô); this would merely introduce a spurious drift to the Lyapunov spectrum.

Let us consider our random vector field in greater detail. Since $\mathbf{f}$ is zero-mean Gaussian, it is uniquely determined by its covariance (or diffusion) tensor

[TABLE]

where $D_{ij}$ is some smooth function. We additionally assume that our vector field is homogeneous and isotropic (in distribution), i.e.

[TABLE]

for all translations $\mathbf{a}\in\mathbb{R}^{n}$ and all rotations $\mathbf{U}\in O(n)$. It is worth noting that the dynamical equation (6) preserves isotropy, while homogeneity is explicitly broken for $\mu>0$.

The statistical symmetries (9) imply that the covariance tensor only depends on the spatial separation. More precisely, we have

[TABLE]

with $\Gamma_{1}$ and $\Gamma_{2}$ denoting smooth functions. We also require that

[TABLE]

These latter constraints arises from obvious positivity conditions, see e.g. (19) in the next section.

It worth making an additional remark before we continue with a stability analysis of our model. It is seen from (6), that for a fixed time $t>0$, we can apply an analysis identical to that produced by Fyodorov and Khoruzhenko [13] in order to find the expected number of equilibria at time $t$. However, this analysis would not tell us anything about stability. The reason for this disconnection between the number of equilibria and stability is clear: in our non-autonomous setting neither the location nor the number of equilibria is preserved in time. Thus, even for systems with a single stable equilibrium there is no guarantee for stability.

Lyapunov exponents and geometric Dyson Brownian motion

The main purpose of this section is to show how our stochastic dynamical system (6) implies a stochastic description of the finite-time Lyapunov exponents. To do so, imagine that we introduce a small perturbation to the initial state of our system, $\mathbf{x}(0)\mapsto\mathbf{y}(0)=\mathbf{x}(0)+\mathbf{r}(0)$. In order to study the Lyapunov exponents, we must study the evolution of the spatial separation $\mathbf{r}(t)=\mathbf{y}(t)-\mathbf{x}(t)$. Taking the integral equation (7) as our starting point, we see that

[TABLE]

For time-intervals in which the separation $\mathbf{r}(t)$ remains sufficiently small, we can use the expansion

[TABLE]

where (by construction) the Jacobian matrix $\mathbf{J}=(J_{ij})=(\partial f_{i}/\partial x_{j})$ is zero-mean Gaussian and white-in-time. Using the regularity conditions imposed on $\mathbf{f}$, it can be seen that the error term in the expansion (13) is subleading for sufficiently small separations. Thus, in the limit of an infinitesimal initial perturbation, the integral equation (12) reduces to

[TABLE]

Writing the integral equation (14) using the usual Stratonovich convention, we get

[TABLE]

Here, $\mathbf{B}(t)$ is a zero-mean matrix-valued Brownian motion constructed such that the Jacobian matrix $\mathbf{J}(\mathbf{x}(t),t)$ is its formal time-derivative.

Before we continue our study of the stochastic equation (15), we need to take a closer look at the noise term. It follows from the covariance tensor (10) that

[TABLE]

We note that the fluctuations of the Brownian motion is independent of the spatial location $\mathbf{x}$. This spatial independence is a consequence of our vector field being both homogeneous-in-space and white-in-time; breaking either of these conditions will inevitable introduce spatial dependence, see section 5. For notational simplicity, let us henceforth use abbreviations

[TABLE]

These constants must satisfy inequalities

[TABLE]

due to the positivity conditions

[TABLE]

Note that the inequalities (18) are consistent with the constraints (11). With the above notation, $\sigma^{2}$ has an interpretation as an overall ‘variance’ of the Jacobian matrix, while $\tau$ is an interpolating parameter such that $\mathbf{f}$ is a conservative vector field for $\tau=+1$ and an incompressible vector field for $\tau=-1/(n+1)$.

We can now return to our stochastic equation (15). Given an initial perturbation $\mathbf{r}(0)$, we may write the general solution of (15) in the form $\mathbf{r}(t)=\mathbf{U}(t)\mathbf{r}(0)$. Here, the evolution matrix $\mathbf{U}(t)$ itself satisfies a stochastic equation

[TABLE]

We may use this equation to find the Lyapunov spectrum. For completeness, let us briefly recall the meaning of Lyapunov exponents in our setting. We know from Osledec’s multiplicative ergodic theorem [25] that the following limit exists,

[TABLE]

The matrix $\mathbf{H}$ is symmetric positive definite, hence we may denote its eigenvalues by $e^{\mu_{1}}\geq e^{\mu_{2}}\geq\cdots\geq e^{\mu_{n}}$. The exponents $\mu_{1}\geq\mu_{2}\geq\cdots\geq\mu_{n}$ are (by definition) called the Lyapunov exponents. If the largest Lyapunov exponent is unique ($\mu_{1}>\mu_{2}$) then we have the asymptotic property

[TABLE]

In other words, we have the familiar property that two points which are initially close tends to separate exponentially fast if the largest Lyapunov exponent is greater than zero, while their separation tends to decay exponentially fast when the Lyapunov spectrum is purely negative.

The Lyapunov exponents are invaluable tool for studies of asymptotic behaviour, but they neglect any finite-time effects which may be physically relevant (see e.g. [19]). Amazingly, it turns out that our model is simple enough to allow us to construct a stochastic description for the finite-time Lyapunov exponents, which includes all finite-time effects. Let $e^{\lambda_{1}(t)}\geq\cdots\geq e^{\lambda_{n}(t)}$ denote the singular values of the evolution matrix $\mathbf{U}(t)$, then we define the finite-time Lyapunov exponents as $\lambda_{1}(t)\geq\cdots\geq\lambda_{n}(t)$. We note that the finite-time Lyapunov exponents are defined such that

[TABLE]

The stochastic equation (20) is a matrix-valued generalisation of a geometric Brownian motion and it follows from Theorem A.1 in Appendix A that the finite-time Lyapunov exponents are given as a geometric Dyson Brownian motion,

[TABLE]

Furthermore, let $\rho(\lambda;t)$ denote the joint probability density function for the finite-time Lyapunov exponents, then it follows from Corollary A.2, that we have a Fokker–Planck–Kolmogorov equation given by

[TABLE]

with initial conditions $\lambda_{1}(0)=\cdots=\lambda_{n}(0)=0$. These two equations give us a full description of the finite-time Lyapunov exponents at any time $t>0$. Figure 1 shows a typical realisation of the finite-time Lyapunov exponents as described by the stochastic equations (24). We see that the finite-time Lyapunov exponents appear to diverge linearly in the long time limit. In the next section, we will look at the asymptotic behaviour of the finite-time Lyapunov exponents and use this to draw a stability-complexity phase diagrams.

May–Wigner transition and the weak-noise limit

In this section we will use the stochastic description of the finite-time Lyapunov exponents (24) and (25) to describe stability of our system. The first step is to look at the long-time behaviour of the finite-time Lyapunov exponents, i.e. Lyapunov exponents themselves. The Lyapunov exponents (23) can be obtained rigorously using several different approaches [18, 5, 23, 24] and can (up to a shift $\mu$) be found in [18] (see Proposition 1 and the remark below). We see no reason to repeat any of these methods here, rather let us provide a heuristic (but physically illuminative) derivation.

Our heuristic derivation is based on the justified (see e.g. [23]) assumption that the Lyapunov spectrum is non-degenerate, $\mu_{1}>\cdots>\mu_{n}$. Consequently, the finite-time Lyapunov exponents diverge as

[TABLE]

To proceed, we pick two times $t>s>0$ and assume that the finite-time Lyapunov exponents at time $s$ are known and pairwise distinct. The evolution of the exponents from time $s$ to time $t$ is given by our stochastic equations (24),

[TABLE]

Now, assuming that $s\gg 1$ and using the asymptotic formulae (26) together with the ordering $\mu_{1}>\cdots>\mu_{n}$ gives

[TABLE]

with exponentially suppressed correction terms. Thus, the asymptotic behaviour of the finite-time Lyapunov exponents is given by

[TABLE]

with

[TABLE]

We note that the asymptotic formula (29) is valid for any $n\geq 1$.

The approximation scheme given by (29) and (30) contains important information about our system. We see that at short times the finite-time Lyapunov exponents are strongly interacting (they repel), while in the asymptotic limit they become independent Brownian motions each with their own drift term $\mu_{i}$. We note that $\lambda_{i}(t)/t\to\mu_{i}$ for $t\to\infty$, hence the drift terms are the Lyapunov exponents. A numerical verification of this phenomenon is seen on figure 1.

Knowledge about the fluctuations of the largest Lyapunov exponent is important for studies of the May–Wigner transition (see e.g. [19]). However, in this paper we will limit ourselves to finding the location of the phase transition (i.e. we determine the phase diagram) and leave a detailed study of the phase transition itself for future work. Our system is said to be stable (unstable) if the largest Lyapunov exponent is less (greater) than zero, $\mu_{1}<0$ ($\mu_{1}>0$). Thus, it is immediate from (30) that our system is stable if $(1+\tau)(n-1)/2<\mu/\sigma^{2}$ and unstable otherwise. It is usual to introduce rescaled units $\hat{\mu}=\mu/\sigma^{2}n$, such that in the large-$n$ limit the May–Wigner transition happens at the critical value $\hat{\mu}_{c}=(1+\tau)/2$. It is worth to compare this result with the linear model (4) studied in [15]. We note that the two models predict a phase transition at the same location, namely $\hat{\mu}_{c}=(1+\tau)/2$. However, there is a crucial difference between the two models: the linear models allow an interpolation parameter $-1\leq\tau\leq+1$ (see [15]), while our present model restrict parameter space to $0\leq\tau\leq+1$. We note that the difference between the two models boils down to whether we first do a linearisation and then impose symmetry constraints (as in [15]) or we first impose symmetry constraints and then linearise (as in this paper).

The above given description of the phase diagram rests upon the assumption that the Lyapunov spectrum is non-degenerate. This is true for $\sigma>0$ but fails for $\sigma=0$ where all Lyapunov exponents are identical and equal to $\mu$. This suggest that there is special behaviour in the weak-noise limit $\sigma\to 0$.

Formally the weak-noise limit is obtained by introducing $\sigma=\varepsilon\tilde{\sigma}$ and $\lambda_{i}(t)=\varepsilon\tilde{\lambda}_{i}(t)$ followed by the limit $\varepsilon\to 0$. With this rescaling, the geometric Dyson Brownian motion (24) reduces to

[TABLE]

Likewise, the Fokker–Planck–Kolmogorov equation becomes

[TABLE]

We recognise this as an ordinary Dyson Brownian motion [9] (see [16] for a review). This stochastic system has a completely different asymptotic behaviour. Most importantly, the interaction between finite-time Lyapunov exponents does not die out in the long-time regime.

Let us consider the stability of the weak-noise system (31). In this case the Lyapunov spectrum becomes trivial $\mu_{1}=\cdots=\mu_{n}=\mu$. This is merely a consequence of the finite-time Lyapunov exponents separating slower than linearly with time, $\lvert\tilde{\lambda}_{i}(t)-\tilde{\lambda}_{j}(t)\rvert=o(t)$ for $t\to\infty$. Thus, we cannot use the Lyapunov exponents to say anything about stability. However, we have more knowledge than the Lyapunov exponents; we know that our system is stable if all the finite-time Lyapunov exponents become negative as time tends to infinity. It is well-known (see e.g. [16]) that the global spectrum of the Dyson process (31) is a growing semi-circle with drift and that the largest finite-time Lyapunov exponent grows as

[TABLE]

where we have introduced the convenient scale $\tilde{\mu}=\mu/\tilde{\sigma}n^{1/2}$. We see that the critical value for stability is $\tilde{\mu}_{c}=((1+\tau)/t)^{1/2}$ which tends to zero as time goes to infinity. We conclude that any weak-noise system (independently of the system size) with a positive decay rate $\mu$ is asymptotically stable. However, from a physical perspective it worth noting that the relevant time-scale for stability is set by $t\gg\tilde{\sigma}n^{1/2}$, hence we might need to wait a while.

Time-effects and persistent auto-correlations

Let us return to our dynamical system (6) described in Section 2. In this section, we look at the effect of weakening the condition that the correlations of our vector field must be white-in-time. More precisely, we will replace the correlation tensor (8) with

[TABLE]

where $D_{ij}(\mathbf{r})$ is given by (10) and $g$ is a smooth monotonically decreasing function. This new dynamical system has the same symmetries as our original system, i.e. the vector field is homogeneous and isotropic in space as well as stationary and time-reversible. However, we will see below that the non-trivial time-dependence given by the function $g$ might introduce persistent time-correlations and thereby have fundamental consequences for the dynamics.

We will proceed as in Section 3 and investigate evolution of a small initial perturbation. Following the same steps as in Section 3, we see that the perturbation $\mathbf{r}(t)$ once again satisfies an integral equation

[TABLE]

where the Jacobian matrix $(J_{ij})=(\partial f_{i}/\partial x_{j})$ is a Gaussian noise matrix. The Jacobian correlation tensor reads

[TABLE]

with $C_{ij;k\ell}$ being a fourth-order tensor given by differentiation of (34). Explicitly, we have

[TABLE]

As in Section 3, rather than working with the Jacobian matrix itself, it is convenient to introduce the integrated process

[TABLE]

By construction, this integrated process is zero-mean Gaussian with auto-correlations

[TABLE]

We see that in the white-in-time limit, $g(t)\to\delta(t)$, the correlations (39) with (37) for the integrated process (38) reduces to the matrix-valued Brownian motion (3). Thus, the integrated process becomes independent of the spatial location and we reclaim our previous results (as we should).

Intuitively, we expect that if the function $g(t)$ is sharply peaked about the origin then we may use the white-in-time limit as an approximation. In order to put a physical handle on this phenomenon, let us introduce the concept of an auto-correlation time for a trajectory. Consider a trajectory $\mathbf{x}(t)$ given as the solution to the differential equation (6) with initial condition $\mathbf{x}(0)=\mathbf{x}$. We define the corresponding auto-correlation time as

[TABLE]

The auto-correlation time (40) may be considered as a measure of ‘trajectory memory’ for a trajectory starting at $\mathbf{x}$. We note that spatial homogeneity in our dynamical system (6) is explicitly broken for $\mu>0$, thus trajectories and therefore the auto-correlation time (40) will in general depend on the starting point $\mathbf{x}$.

Let us take a closer look at the auto-correlation time (40). We know from the constraints (11) that $\Gamma_{1}^{\prime}(0)\leq 0\leq\Gamma_{1}(0)$. In other words, the function $\Gamma_{1}$ is decreasing in a neighbourhood of zero. From a physical perspective, it is natural to assume that this extends to a global property, i.e. spatial correlations decays with separation. With this assumption, we have the upper bound $\Gamma_{1}(x)\leq\Gamma_{1}(0)$ which ensures that the auto-correlation time is finite throughout space as long as $g$ decays sufficiently fast. The white-in-time approximation is then reasonable on time-scales $t\gg T$. On the other hand, if $g$ decays too slowly then the auto-correlation time might be divergent, which indicates persistent time-correlations. The worst case scenario is when $g$ is constant, which corresponds to the Fyodorov–Khoruzhenko model [13]. In this scenario, it is evident from (40) that the auto-correlation time will diverge if the trajectory remains in a small region of space for a long time, which is the case for trajectories starting in the neighbourhood of an equilibrium.

Conclusions and open problems

In this paper we have introduced a generic model for study of stability in large complex systems. Our model may be considered a non-autonomous version of a model recently introduced by Fyodorov and Khoruzhenko [13]. However, the behaviour of our model differs considerably from the behaviour of the Fyodorov–Khoruzhenko model due to its explicit time-dependence. Most importantly, we have seen that the stability-complexity phase transition (or May–Wigner transition) happens at a scale proportional to the system size $n$, while the same transition happens at scale proportional to $n^{1/2}$ in the Fyodorov–Khoruzhenko model. These types of scaling relations are important as they present one of the most promising avenues towards experimental tests. We note that a recent paper [4] claims to have experimentally confirmed the scaling $n^{1/2}$ from a closely related model [27] occurring in the description of neural networks.

From a more theoretical perspective, we have seen that it is possible to construct a stochastic description of the finite-time Lyapunov exponents for spatial homogeneous and isotropic dynamical systems under the assumption that the temporal fluctuations are sufficiently rapid. If the temporal fluctuations are too slow then the stochastic description breaks down and trajectories with persistent time-correlations become important for the dynamics.

Moreover, we have shown that there is a close connection between the non-linear dynamical system studied in this paper and the linear model (4) studied in [15]. That is a connection to matrix-valued stochastic equations of multiplicative type and geometric Dyson Brownian motions. In our opinion, these models are by themselves interesting. For instance, it would be intriguing to study the global spectrum of the finite-time Lyapunov exponents, which are related to the theory of free stochastic equation. This is of particularly interest since the short-time behaviour of the finite-time Lyapunov exponents must be very different from the asymptotic behaviour. It would also be worthwhile to study local correlations; in particular the largest exponent which determines the behaviour at the phase transition (see e.g. [19] where fluctuations of of the largest Lyapunov exponent was used to show that the May–Wigner transition for symmetric Gaussian community matrices is a third-order phase transtion). Another problem worth pursuing is a study of discrete time models, which would relate to products of random matrices [8, 1].

Acknowledgement:

We thank Peter Forrester for comments on a first draft of this paper, Yan Fyodorov for sharing an early version of the paper [6], and Andre Peterson for bringing the paper [4] to our attention. The author acknowledge financial support by ARC Centre of Excellence for Mathematical and Statistical frontiers (ACEMS).

Appendix A Geometric Dyson Brownian motion

Let $\mathbf{B}(t)$ be a $\mathfrak{gl}_{n}(\mathbb{R})$-valued isotropic Brownian motion, i.e. $\mathbf{B}(t)=[B_{ij}(t)]$ is an $n$-by-$n$ matrix those entries are Brownian motions such that [18, 15]

[TABLE]

with $a\pm c\geq 0$ and $a+nb+c\geq 0$. Consider the matrix-valued generalisation of a geometric Brownian motion,

[TABLE]

where $\mathbf{I}$ denotes the identity matrix and $\mathbf{U}(t)$ is a stochastic process which is $\textup{GL}_{n}(\mathbb{R})$-valued (a.s.) in virtue of the Stroock–Varadhan support theorem.

Let $e^{\lambda_{1}(t)},\ldots,e^{\lambda_{n}(t)}$ denote the singular values of $\mathbf{U}(t)$ at time $t>1$. We define the finite-time Lyapunov exponents as ${\lambda_{1}(t)},\ldots,{\lambda_{n}(t)}$ and we observe that these exponents are finite and pairwise distinct for any fixed time $t$ (a.s.) unless $a+c=b=0$. For $a+c=b=0$, the Brownian motion $\mathbf{B}(t)$ is anti-symmetric, hence $\mathbf{U}(t)$ belongs to orthogonal group (a.s.). Consequently, all the finite-time Lyapunov exponents are trivially equal to zero in this case.

Theorem A.1**.**

Let $\mathbf{U}(t)$ be a matrix-valued geometric Brownian motion defined as the solution to (42), then the corresponding finite-time Lyapunov exponents are solutions to the coupled stochastic equations

[TABLE]

with $B_{1}(t),\ldots,B_{n}(t)$ being independent standard Brownian motions.

Corollary A.2**.**

Let $\rho(\,\cdot\,;t)$ denote the joint probability density of the finite-time Lyapunov exponents at time $t>0$, then it satisfies the Fokker–Planck–Kolmogorov equation

[TABLE]

where

[TABLE]

Remark A.3**.**

We will refer to (43) as a geometric Dyson Brownian motion. Note that a geometric Dyson Brownian motion with $b=c=\mu=0$ was considered in [24], while we considered the case $b=0$ in [15] where also a complex-valued version was introduced and studied. Actually, the Fokker–Planck–Kolmogorov equation (44) is not explicitly mentioned in [24], while the stochastic formulation (43) is not explicitly mentioned in [15], but these two formulations are naturally closely related.

Proof of Theorem A.1.

We are interested in the finite-time Lyapunov exponents which relates to the singular values of $\mathbf{U}(t)$, thus it is convenient to introduce a Wishart-type matrix

[TABLE]

For a small time increment $dt>0$, it follows from (42) that

[TABLE]

As $\mathbf{S}(t)$ is a symmetric matrix, we may proceed as for ordinary Dyson Brownian motion [9] and use perturbation theory. If $s_{1}(t),\ldots,s_{n}(t)$ denote the eigenvalues of $\mathbf{S}(t)$, then we have

[TABLE]

We recall that the eigenvalues $s_{1}(t),\ldots,s_{n}(t)$ are pairwise distinct (a.s.). Furthermore, it follows from (47) and (41) that

[TABLE]

and

[TABLE]

Inserting (49) and (50) into the perturbative expansion (48) gives

[TABLE]

which after some standard manipulations of the sums may be rewritten as

[TABLE]

It remains to change variables from the eigenvalues $s_{1}(t),\ldots,s_{n}(t)$ to the finite-time Lyapunov exponents $\lambda_{1}(t),\ldots,\lambda_{n}(t)$. By construction, we have $s_{i}(t)=e^{2\lambda_{i}(t)}$, thus it follows by Itô’s formula that

[TABLE]

and consequently that

[TABLE]

The final observation is that due to (41), we may make the replacement $B_{ii}(t)\mapsto(a+b+c)^{1/2}B_{i}(t)$ where $B_{1}(t),\ldots,B_{n}(t)$ are independent standard Brownian motions. ∎

Corollary A.2 follows from the Theorem A.1 using standard arguments, see e.g. [28]. Perhaps the easiest approach is to first introduce a test function $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ which is smooth with bounded derivatives. By means of Taylor expansion, it is straightforward to show that for any such test function, we have

[TABLE]

where $L^{*}$ is the adjoint geometric Dyson operator, i.e.

[TABLE]

The corollary then follows using integration by parts.