Isotropic Brownian motions over complex fields as a solvable model for May-Wigner stability analysis

Introduction

A close link between stochastic calculus and random matrix theory exists since the seminal work of Dyson , who introduced a stochastic process on matrices which greatly simplifies the task to implement various canonical ensembles. Since its introduction, this method, now known as Dyson Brownian motion, has found a wide variety of applications in both physics and mathematics, see e.g. . Two important examples of matrix-valued stochastic models in physics concern the passive advection in fully developed turbulence , where one is interested in the stability of the flow, and the phase-coherent transport through disordered wires , where one encounters Anderson localization. While Dyson’s model has additive noise, these models have multiplicative noise, which gives rise to scaling regimes with statistics that differ from the usual Wigner–Dyson type. Not surprisingly, there are numerous close connections between these two classes of models with multiplicative noise .

In this work, we identify a similarly close relation starting from a paradigmatic model, sometimes referred to as isotropic Brownian motions , isotropic stochastic flows , or matrix-valued multiplicative diffusions . For motivation, we consider this model in the context of a May–Wigner-like stability analysis of large complex systems, as they occur in a wide variety of settings in physics, biology and beyond. We show that the finite-time Lyapunov exponents share the same statistics as those appearing in the transport through a quantum wire with chiral symmetry. This connection is surprising inasmuch as the exponents for the quantum wire are constrained by a condition related to flux conservation, which does not have a counterpart for the isotropic Brownian motion. The relation between the models allows us to find an exact solution for motions over complex fields, owing to the fact that the statistics in the corresponding transport problem is determined by a Hamiltonian of the Calgero-Sutherland type. Furthermore, we identify connections to products of random matrices, including Newman’s square law and biorthogonal ensembles , as well as Hermitian matrix models with an external source . These findings help us to establish a simple stability phase diagram, expressed in terms of an effective parameter that incorporates the rotational and gradient components of the flow, and defines a convenient starting point for more detailed investigations of the phase transition.

First, let us recall some basics about May–Wigner-like stability analysis. We consider a generic autonomous dynamical system given by

For sufficiently complex systems it is reasonable to assume (at least as a toy model) that $J=[J_{ij}]$ is a “maximally random” matrix, up to symmetries stemming from physical considerations. Such models were introduced in a seminal paper by May . The original paper was mainly focussing on the stability of large ecosystems, where it challenged the (at the time) common folklore that the stability of complex systems increases with the size of the system. Similar ideas as those introduced in have since been applied to many other systems, including e.g. machine learning , finance , and neural networks . There is also a close connection to what is known as random landscapes .

Rather than an expansion around a fixed point, we imagine to evaluate the gradient matrix along a trajectory, see e.g. . In this case, the gradient matrix depends on the position on the trajectory and therefore indirectly on time. With May’s model in mind, we write

with some time-dependent coupling matrix $J_{ij}(t)$, to be further specified below. This is such that we have an attractive trajectory with relaxation rate $\mu$ when the coupling strength vanishes, i.e. when $\sigma=0$. Our main assumption is that the time-dependent coupling $J_{ij}(t)$ can be treated as noise, which we let be white in time. Such an approximation is valid in the limit of strong chaoticity, where the correlation time is short.

We therefore consider the system (3) driven by Gaussian white noise, with the only constraint that the system must be isotropic. Originally, this type of process arose in models for passive advection in fully developed turbulence (here isotropy appears as an axiom stemming from K41 theory), see e.g. . The insight that temporal correlations may be neglected for turbulent velocity fields in the limit of high Reynolds numbers dates back to Kraichnan . This type of isotropic models has also been considered by the mathematical community, see e.g. . Later methods originating from random matrix theory and quantum field theory have been applied as well .

We recall that the white-noise limit of processes like (3) depends on the regularisation (the one-dimensional case, which is a geometric Brownian motion, provides a well-known example). Thus one needs to choose an interpretation, or alternatively impose additional constraints (see e.g. the appendix in ). For consistency with the existing literature we choose the Stratonovich interpretation as our starting point,

where $B(t)=[B_{ij}(t)]$ is a matrix-valued Brownian motion assumed to be statistically invariant under unitary similarity transformations. This is the type of process that we refer to as isotropic Brownian motion.

For the practical implementation, on the other hand, it will be convenient to work in Itô convention, in which the process (4) becomes

Before we go on to the main results of this paper, let us recall a few basic properties of isotropic Brownian motions. Given an initial condition $y(0)$, the solution to the system (5) may in matrix notation be written as

where $\Pi_{\mu}(t)$ is a random evolution matrix, which itself satisfies an Itô equation similar to (5),

with initial condition $\Pi_{\mu}(0)=\mathbf{I}$. Moreover, with $\Pi(t,s)$ denoting the evolution from time $s$ to time $t$, we have the almost-sure property $\Pi_{\mu}(t,s)\Pi_{\mu}(s,0)=\Pi_{\mu}(t,0)$ for all $t\geq s\geq 0$. Per definition, the solution is given as an ordered product,

where $\delta t:=t/n$ is the time regularisation and $\{X^{(k)}\}$ is a family of independent random matrices given as the Brownian increments, $X^{(k)}\delta t^{1/2}=B(k\delta t)-B((k-1)\delta t)$. We could, of course, also write the random evolution matrix as a time-ordered exponential, but we will not need that formulation here. It follows that finding the random evolution matrix (5) boils down to evaluating a product of independent random matrices; a topic which in the discrete setting has received considerable recent attention .

In this paper, we consider evolutions over real as well as complex vector fields; thus, the random matrix $X(=X^{(k)})$ may be either real or complex. Adopting notation from random matrix theory, the real and complex case will be denoted by an index $\beta=1$ or $\beta=2$, respectively, whenever the distinction is important.

Isotropy implies that $U^{-1}XU\stackrel{{\scriptstyle d}}{{=}}X$ for all rotations $U$, i.e. $U\in\operatorname{O}(N)$ for $\beta=1$ and $U\in\operatorname{U}(N)$ for $\beta=2$. Up to a shift proportional to the identity, the most general matrix $X$ consistent with the isotropy constraint and Gaussianity is a matrix from the Gaussian elliptic ensemble (see appendix A), i.e. a random matrix distributed with respect to the density

with $Z$ denoting the normalisation constant and $\tau\in(-1,1)$ denoting an interpolation parameter between Hermitian and skew-Hermitian matrices. In general, the matrix $X$ is non-Hermitian (Ginibre for $\tau=0$) and the random evolution matrix $\Pi_{\mu=0}(t)$ is therefore a diffusion on the general linear group, $\operatorname{GL}_{N}$. The special case is $\tau=-1$, where $X$ is skew-Hermitian and thus $\Pi_{\mu=0}(t)$ belongs to the unitary subgroup. Physically $\tau$ measures the rotational versus gradient nature of the flow, with $\tau=1$ representing a pure gradient flow.

Perhaps the simplest question we may ask about isotropic Brownian motions regards the large-time asymptotic of the strain matrix, $S(t):=\Pi_{\mu}^{\dagger}(t)\Pi_{\mu}(t)$. It is a well-known consequence of Oseledec’s multiplicative ergodic theorem that the matrix $(\log S(t))/2t$ stabilises for almost all realisations as $t$ tends to infinity. The eigenvalues of this limiting matrix are the so-called Lyapunov exponents and are used as a measure of stability. A system (4) is said to describe an attractive trajectory if

which holds true as long as the largest Lyapunov exponent is less than zero. Thus, the value of the largest Lyapunov exponent determines stability, while its fluctuations are essential for the stability-instability transition. We note that there exist a few general theorems claiming Gaussian fluctuations of the largest Lyapunov exponent, see e.g. . However, none of these are directly applicable in the limit of high dimensionality, $N\to\infty$. The heuristic understanding of this breakdown is that the theorems implicitly require a gap between the largest and second largest Lyapunov exponent, while the Lyapunov spectrum often becomes continuous when $N\to\infty$.

For isotropic Brownian motions over real vector fields, explicit formulae for the Lyapunov spectrum were obtained 30 years ago . However, these exponents do not capture the fluctuations in the dynamics at finite times. A more challenging problem is to study the statistical properties of the spectrum of the strain matrix as a function of time, and thereby gain information about the finite-time Lyapunov exponents. Such results can be used to further study the double scaling limit where both the time $t$ and the dimension $N$ tend to infinity. This adds interest since non-trivial scaling regimes are expected.

The remainder of this paper is organised as follows. In section 2 we show how to pass from the matrix-valued Itô equation (5) to the Fokker–Planck equation for the finite-time Lyapunov exponents, including the case of complex fields. In section 3 we show that the Fokker–Planck equation coincides with the one encountered in the phase-coherent transport problem, which is exactly solvable in the complex case, and provide an explicit, compact expression for the joint probability density function. Based on this expression, we also explain the connection to biorthogonal ensembles and Hamiltonian models with an external source. The final section is devoted to the discussion in terms of May–Wigner stability analysis, and a description of further open problems and possible applications. In the appendix we recall the relation between isotropy and the Gaussian elliptic ensemble.

Fokker–Planck equation for finite-time Lyapunov exponents

Our goal in this section is to find the Fokker–Planck equation for the finite-time Lyapunov exponents. To do this, we need to look at the time evolution of the eigenvalues of the strain matrix $S(t)$. The easiest way to proceed is to use the product formulation (8).

Let us denote the eigenvalues of the strain matrix by $s_{1}(t),\ldots,s_{N}(t)$. Since the strain matrix is Hermitian, it follows from ordinary perturbation theory that

If we consider evolutions over real vector fields we have $X^{*}=X$, but we keep the notation general to also account for evolutions over complex fields.

In order to find the Fokker–Planck equation for the eigenvalues of the strain matrix, we use the product formulation (8) and write down a recursive formula for the matrix density of the random evolution matrix $\Pi_{\mu}(n\delta t)$. Upon expansion in $\delta t$ and taking the limit $\delta t\to 0$ ($n\to\infty$), this recursion produces a differential equation for the matrix density. Finally, the Fokker–Planck equation is obtained by taking the expectation value of the empirical density with respect to the matrix density. This gives

where $\rho_{t}(s_{1},\ldots,s_{N})$ is the joint probability density function for the eigenvalues at time $t$. The next step is to note that (11) and (12) imply

which by insertion in (13), at least in principle, provides an explicit expression for the Fokker–Planck equation for the eigenvalues of our strain matrix.

With the random matrix $X$ distributed according to the density (9), we have covariances

Introducing these covariances into the above-given formulae, we obtain the Fokker–Planck equation

Here, the first line on the right-hand side represents a diffusive term with diffusion constant given as $\kappa=(1+\tau)\sigma^{2}/2$, while the second line represents a drift term.

The final step is a change of variables from the eigenvalues of the strain matrix to the exponents $\{\lambda_{k}:=\frac{1}{2}\log s_{k}\}$. After some standard manipulations, we find the Fokker–Planck equation for the exponents,

with repulsion term and initial condition given by

respectively. Here, the latter condition originates from the fact that the random evolution matrix is required to be equal to unity at $t=0$.

which more closely resembles the notation chosen by Dyson .

The Fokker–Planck equation (19) is our main result in this section. In the next section, we will show that this equation is exactly solvable for $\beta=2$. However, before this a few remarks are in order.

We first note that the diffusive term in (19) vanishes as $\tau\to-1$; recall that $\kappa=(1+\tau)\sigma^{2}/2$. This absence of diffusion occurs since the random evolution matrix $\Pi_{\mu=0}(t)$ is unitary when $\tau=-1$. It follows that the eigenvalues of the strain matrix $S(t)=\Pi_{\mu}^{\dagger}(t)\Pi_{\mu}(t)$ are all identical and equal to $-\mu t$.

They are independent of $\beta$ and equidistantly spaced over the interval $(-\kappa N-\mu,\kappa N-\mu)$. Evidently, if we take $\sigma^{2}=1/N$ and thus $\kappa=(1+\tau)/2N$, we have convergence of the global spectral density to a “square law” on an interval of length $1+\tau$, centred at $-\mu$. This applies to the iterated limit where $t\to\infty$ followed by $N\to\infty$, for which the square law was originally pointed out by Newman . More recent results for products of random matrices lead us to believe that this law is independent of the order of the limits .

The convergence of the global spectrum allows us to establish a stability phase diagram in the large $N$ and $t$ limit. With the aforementioned scaling, the system is stable if $(1+\tau)/2<\mu$ and unstable if $(1+\tau)/2>\mu$. The Lyapunov exponents themselves can however not be used to describe the finer structure of the phase transition; this requires information about their fluctuations. The Fokker–Planck equation (19) is a good starting point for a study of such fluctuations.

Heuristically, the emergence of the Lyapunov exponents (22) from the Fokker–Planck equation (19) may be understood by realising that for finite $N$, the eigenvalues separate exponentially fast compared with the eigenvalue repulsion. Thus, with the ordering $\lambda_{1}\ll\cdots\ll\lambda_{N}$, in the long-time limit

With this approximation the Fokker–Planck equation (19) turns into $N$ uncoupled heat equations. The exponents are seen to be independently Gaussian distributed, and in the long-time limit agree with (22).

The benefit of the Fokker–Planck equation (19) is that it provides information about the statistical properties of the exponents at all times $t$. This is in contrast to Newman’s method and its extensions, which only apply when $\kappa t\gg N$.

Solving the Fokker–Planck equation over complex fields

In this section, we show that the Fokker–Planck equation (19) is exactly solvable in the complex case, i.e. for $\beta=2$. We exploit that a specific version of the Fokker–Planck equation appears for a matrix model describing the phase-coherent transport properties of quasi-one-dimensional disordered wires with chiral symmetry . In this setting, one investigates the so-called transfer matrix $M$, which is a $2N\times 2N$ dimensional matrix that obeys the symplectic constraint $M^{\dagger}\sigma_{1}M=\sigma_{1}$, while chirality imposes $\sigma_{3}M\sigma_{3}=M$ (both conditions are expressed in terms of Pauli matrices $\sigma_{i}$). Due to the symplectic structure, the eigenvalues $\exp(2x_{n})$ of $M^{\dagger}M$ occur in reciprocal pairs, $(x_{n},-x_{n})$. Chirality enforces a block-diagonal structure $M={\rm diag}\,(A,(A^{\dagger})^{-1})$, so that the exponents $x_{n}$ arise from $AA^{\dagger}$ while the exponents $-x_{n}$ arise from $(AA^{\dagger})^{-1}$. The multiplicative law of transfer matrices with Gaussian statistics then results in the Fokker–Planck equation (19) with $\tau=\mu=0$, and $\lambda_{n}=x_{n}$ identified with one of the two sets of the eigenvalues (in this context, the Fokker–Planck equation is known as the DMPK equation). In our case the parameters $\tau$ and $\mu$ are finite, but this can be accounted for by rescaling and shifting. An important difference between the physics underlying the two models is that the transport is dominated by the transport exponent $x_{n}$ nearest zero, while the stability of the flow depends on the largest Lyapunov exponent $\lambda_{n}$.

Here, we obtain the exact solution following , and then bring the result into a compact form which more directly reveals the long time asymptotics of the Lyapunov exponents.

We first note that the drift $\mu$ only introduces an overall shift of the spectrum. This allows us to simplify notation by setting $\mu=0$ in the calculations, and then reintroduce the shift in the final result. After this simplification, the key idea to solve (19) is to parametrise the joint density as

where $\psi_{t}$ is some (wave) function, and $\nu=(\nu_{1},\ldots,\nu_{N})$ is a given initial condition. With this parametrisation and the notation from above, the Fokker–Planck equation turns into a Schrödinger equation in imaginary time,

which turns out to be of Calogero–Sutherland type . The evaluation of the potential is straightforward. By insertion of the definition of $\Omega(\lambda)$ from (20), we find

Here, the second term on the right-hand side is seen to be a constant by exploiting an identity for cyclic sums. For distinct $\lambda_{i},\lambda_{j},\lambda_{k}$, we have

Consequently the pair interaction vanishes for $\beta=2$, and the Hamiltonian in (25) becomes that of $N$ free particles. This is the feature which ensure solvability for $\beta=2$.

Now, we can return to the Schrödinger equation (25). For the rest of this section we will restrict our attention to the case $\beta=2$ only, and leave out the index for notational simplicity.

where $U=\kappa(N+1)N(N-1)/3$ is the constant contribution to the potential energy (29), while each $g_{j}^{t}(\lambda)$ satisfies a heat equation

This can be verified by inserting the wave function (30) into the Schrödinger equation (25). The solution to the heat equation (31) is a Gaussian,

Now, as a consequence of (24), we know that the joint density is

assuming a non-singular initial condition $\nu_{1}<\cdots<\nu_{N}$. We note that (33) reduces to a product of Dirac delta functions in the $t\to 0$ limit, as required.

We are interested in the singular initial condition $\nu_{1}=\cdots=\nu_{N}=0$, which may be obtained from (33) by successive use of l’Hôpital’s rule. Upon reordering of rows or columns, we find

To further simplify the expression for the joint density, we first need to make an observation about the Lyapunov exponents from the previous section. At zero drift ($\mu=0$), we have

we arrive at a surprisingly simple expression for the joint density,

The joint density (37) describes a special type of biorthogonal ensembles which have been coined polynomial ensembles . This type of ensembles has recently gained renewed attention partly due their prominent rôle in the study of the singular values of random matrix products, see e.g. . In the case of product ensembles, the Gaussian weights within the second determinant in (37) are replaced with weights given in terms of Meijer $G$-functions. In this way, isotropic Brownian motions fit neatly into the picture of recent developments regarding products of random matrices.

There is also a direct relation between our joint density (37) and another well-known matrix ensemble, the Gaussian Unitary Ensemble (GUE) with an external source (see with references and also ). This ensemble is based on an $N\times N$ random Hermitian matrix, $H$, distributed with respect to the measure

where $Z$ is a normalisation constant, $2\kappa t$ is the variance, and $A$ is an $N\times N$ Hermitian external source matrix which without loss of generality may be taken to be diagonal, $A:=\operatorname*{diag}(a_{1},\ldots,a_{N})$. The joint density for the eigenvalues of $H$ is obtained by integrating out irrelevant degrees of freedom using the Harish-Chandra–Itzykson–Zuber integral . One then finds

The similarity between (39) and (37) is immediately recognised. Thus the Lyapunov exponents in (37) may be reinterpreted as equidistantly spaced eigenvalues of an external source matrix, as considered e.g. in . While relations between ensembles with an external source and non-intersecting Brownian motions are not new (see and references within), the isotropic Brownian motions studied in this paper provide an example where the external source arises naturally; rather than from an imposed boundary condition.

Conclusions and open problems

In this paper we have studied a family of stochastic processes with isotropic matrix-valued multiplicative noise. We have shown that it is possible to formulate a Fokker–Planck equation for the finite-time Lyapunov exponents, and that this Fokker–Planck equation is exactly solvable for evolutions over complex fields, where they give rise to a biorthogonal ensemble. We motivated this stochastic process by its relation to a May–Wigner-like stability analysis for trajectories in an $N$-dimensional space, characterized by a mean relaxation rate $\mu$, a noise strength $\sigma$, and a parameter $\tau\in(-1,1)$ which characterizes the rotational nature of the flow ($\tau=1$ represents a pure gradient flow). For a noise variance $\sigma^{2}=1/N$, the infinite-time Lyapunov spectrum has compact support, and converges to a uniform distribution on an interval with length $1+\tau$ centred at $-\mu$. This allows us to establish a phase diagram for the system in the large-$N$ limit, according to which trajectories are stable if $(1+\tau)/2<\mu$ and unstable if $(1+\tau)/2>\mu$.

The heuristic argument at the end of section 2 suggests that the fluctuations of the largest Lyapunov exponent become Gaussian in the limit $\kappa t\gg N$, with more rigorous formulations of this statement following from . It is more challenging to investigate the fluctuations away from this limit, where they are expected to be non-Gaussian. Here, the Fokker–Planck equation established in this paper provides a good starting point. In particular, we note that $\kappa t\ll N$ results in a limit where adjacent stability exponents are close compared to their correlations length. With this in mind, we may substitute $\tanh(\lambda_{j}-\lambda_{i})\approx(\lambda_{j}-\lambda_{i})$ in (21), which transforms the Fokker–Planck equation into an ordinary Dyson diffusion. This suggests that the largest Lyapunov exponent follows the so-called Tracy–Widom law. Similar approximations may be made starting with the joint density function from section 3. We note that this conjectural transition from a classical random matrix law (Tracy–Widom distribution) to a Gaussian law is not completely unfamiliar; in it was shown that the largest eigenvalue (in terms of absolute value) for a product of complex Ginibre matrices undergoes a transition from a Gumbel distribution to a log-normal distribution.

Future work may be directed towards establishing this transition, and extending it to a detailed description of the stability-instability phase transition, including the critical exponents. From a more mathematical perspective, there are many other intriguing questions worth pursuing beyond the fluctuations of the largest Lyapunov exponent. This includes, but is not limited to, a study of the global spectrum as a function of time, and the local correlations in the bulk as well as near the edge at fixed time $t$.

We like to thank G. Akemann, P. J. Forrester, and M. Kieburg for useful discussions. JRI acknowledge financial support by ARC Centre of Excellence for Mathematical and Statistical Frontiers.

Appendix A Isotropic measures and elliptic ensembles

In this appendix we briefly recall the construction of isotropic Gaussian matrix measures and their relation to Gaussian elliptic ensembles. We will consider the real ($\beta=1$) and the complex ($\beta=2$) case separately.

for all $U\in\operatorname{O}(N)$. Consequently, any isotropic measure has a covariance tensor given by

with $a,b,c$ denoting constants unaffected by orthogonal similarity transformations. To understand the interpretation of these constants, we note that

Thus, we may write the random matrix $X$ as a sum

where $H$ is a symmetric matrix with standard Gaussian entries (i.e. GOE), $A$ is a skew-symmetric matrix with standard Gaussian entries (i.e. skew-GOE), and $\xi$ is a standard Gaussian random variable. Setting $a=0$, $b=1$, and $c=\tau\in(-1,1)$ result in the elliptic density (9) with $\beta=1$.

where $a,b,c,d$ are constants. Similar to the real case, we will look at correlations to obtain an interpretation of these constants. We then find

and can therefore write the random matrix $X$ as a sum

where $H$ is a Hermitian matrix with standard complex (real on the diagonal) Gaussian entries (i.e. GUE), $A$ is a skew-Hermitian matrix with standard complex (imaginary on the diagonal) Gaussian entries (i.e. skew-GUE), while $\xi$ and $\eta$ are standard real Gaussian random variables. Analogously to before we set $a=0$, $b=\tau\in(-1,1)$, $c=0$, and $d=1$, which results in a random matrix with density (9).

References