What drives transient behaviour in complex systems?

TL;DR

This paper investigates the transient dynamics of complex systems modeled by non-linear ODEs, revealing new regimes and statistical distributions of transient trajectories through an extended May-Wigner random matrix framework.

Contribution

It introduces a novel stable-transient regime, extends the May-Wigner model to include transient behavior, and derives exact distributions for transient trajectories using advanced random matrix theory.

Findings

Identification of a stable-transient regime based on initial response.

Exact calculation of transient trajectory distributions, including Gaussian and Tracy-Widom.

Connection of transient behavior to eigenvector non-orthogonality and an extended May-Wigner model.

Abstract

We study transient behaviour in the dynamics of complex systems described by a set of non-linear ODE's. Destabilizing nature of transient trajectories is discussed and its connection with the eigenvalue-based linearization procedure. The complexity is realized as a random matrix drawn from a modified May-Wigner model. Based on the initial response of the system, we identify a novel stable-transient regime. We calculate exact abundances of typical and extreme transient trajectories finding both Gaussian and Tracy-Widom distributions known in extreme value statistics. We identify degrees of freedom driving transient behaviour as connected to the eigenvectors and encoded in a non-orthogonality matrix $T_0$. We accordingly extend the May-Wigner model to contain a phase with typical transient trajectories present. An exact norm of the trajectory is obtained in the vanishing $T_0$ limit where it describes a normal matrix.

Equations

$$\frac{d{x}_i}{dt}=f_i(x_1...x_N),i=\mathrm{1...}N,$$

$$x_i=x_i^{\ast }+y_i,$$

$$\frac{d}{dt}{y}_i(t)=\sum _{j=1}^N{M}_{ij}{y}_j(t),$$

$$\{\begin{array}{l}\dot{x_1}=-x_1+x_2^3\\ \dot{x_2}=\alpha x_1-2x_2-x_1{x}_2-x_2^4\end{array},$$

$$M_{|x^{\ast }}=\left(\begin{array}{cc}-1& 0\\ \alpha & -2\end{array}\right),$$

$$\{\begin{array}{ll}\dot{y_1}& =-y_1\\ \dot{y_2}& =\alpha y_1-2y_2\end{array}.$$

$$M=-\mu +X,$$

$$P_{\beta }(X)[dX{]}_{\beta }=c_{\beta }\exp \left(-\frac{\beta N}{2{\sigma }^2}\text{Tr}{X}^{†}X\right)[dX{]}_{\beta },$$

$$\overline{O}={⟨O(X)⟩}_X=\int [dX{]}_{\beta }{P}_{\beta }(X)O(X).$$

$$A=\underset{t\ge 0}{\max }\frac{|y(t){|}^2}{|y(0){|}^2},$$

$$R=\frac{1}{|y(0){|}^2}\underset{t\to 0}{\lim }\frac{d|y(t){|}^2}{dt},$$

$$|y(t)⟩=e^{Mt}|y_0⟩,$$

$$R=\frac{⟨y_0|(M^{†}+M)|y_0⟩}{⟨y_0|y_0⟩},$$

$$O_{\text{av}}$$

$$O_{\max }$$

$${\rho }_M(x,y)=\frac{1}{\pi {\sigma }^2}\theta \left({\sigma }^2-(x+\mu {)}^2-y^2\right),$$

$$R_{\max }=\underset{y_0}{\max }\frac{⟨y_0|(M^{†}+M)|y_0⟩}{⟨y_0|y_0⟩}={\lambda }_{\max }\left(M^{†}+M\right).$$

$$\underset{N\to \infty }{\lim }{⟨R_{\max }⟩}_X=-2\mu +2{\mu }_T,$$

$$\overline{g_{\max }}(r)={⟨\underset{y_0}{\max }g(r)⟩}_X,\overline{g_{\text{av}}}(r)={⟨g(r)⟩}_{X,y_0}.$$

$$\overline{g}=c_{\beta }\int [dX{]}_{\beta }{P}_{\beta }(X)\delta (r-R(y_0)),$$

$$\overline{g}(r)=\frac{1}{2\pi }\int dp{e}^{ip(r+2\mu )}{e}^{-\frac{2{\sigma }^2{p}^2}{N}},$$

$$\overline{g}(r)=\frac{1}{\sqrt{2\pi {\sigma }_{\beta ,R}^2}}{e}^{-\frac{(r+2\mu {)}^2}{2{\sigma }_{\beta ,R}^2}},$$

$$\overline{g_{\max }}(r)$$

$$=\frac{d}{dr}{⟨\theta \left(\mu +\frac{r}{2}-{\lambda }_{\max }\left(\frac{X^T+X}{2}\right)\right)⟩}_X,$$

$$\overline{g_{\max }}(r)=\frac{d}{dr}{F}_{N,\beta }\left(\frac{\sqrt{N}}{\sigma }\left(\mu +\frac{r}{2}\right)\right).$$

$$\overline{N_{\max }}={\int }_0^{\infty }\overline{g_{\max }}(r)dr,\overline{N_{\text{av}}}={\int }_0^{\infty }\overline{g_{\text{av}}}(r)dr,$$

$$\overline{N_{\max }}(\mu )$$

$$\overline{N_{\text{av}}}(\mu )$$

$$\overline{N_{\text{av}}}\sim e^{-\beta N{\left(\frac{\mu }{{\mu }_T}\right)}^2}.$$

$$F_{N,\beta }\left(\sqrt{2N}-\frac{\sqrt{N}\Delta }{\sigma }\right)\sim e^{-\beta N^2\varphi (\Delta /\sigma -\sqrt{2})},$$

$$F_{N,\beta }\left(\sqrt{2N}-\frac{\sqrt{N}\Delta }{\sigma }\right)\sim F_{\beta }\left(-N^{\frac{2}{3}}\frac{\sqrt{2}\Delta }{\sigma }\right),$$

$$\displaystyle\overline{N_{\max}}(\mu_{T}-\Delta)\sim 1-\left\{\begin{tabular}[]{ll}$e^{-\beta N^{2}\phi(\Delta/\sigma-\sqrt{2})}$&, $\Delta\sim\mathcal{O}(1)$\\ $F_{\beta}\left(-N^{\frac{2}{3}}\frac{\sqrt{2}\Delta}{\sigma}\right)$&, $\Delta\sim\mathcal{O}(N^{-\frac{2}{3}})$\end{tabular}\right.,$$

Full text

One of the key problems in studying complex systems is answering the question of stability. The standard linearization approach relies heavily on the large time asymptotics and can be misleading for intermediate times. This is especially pronounced when systems develop transient behaviour – the analysis based on eigenvalues loses its significance and different approach is needed. This shortcoming in physics literature can be traced back to the work of Orr Orr (1907) in the hydrodynamical context. Since then, similar ideas were revived in the context of fluid dynamics Trefethen et al. (1993); Schmid (2007), plasma physics Camporeale et al. (2009); Ratushnaya and Samtaney (2014), diffusion in porous media Rapaka et al. (2008) or pattern formation Ridolfi et al. (2011); Biancalani et al. (2017). Further motivation for this work is rooted in the ecological literature on biological networks May (1972); Neubert and Caswell (1997); Tang and Allesina (2015).

The physical mechanism of transient behaviour is both relatively simple and quite general. It needs the system’s components to interact asymmetrically and be stabilized by an effective dissipation mechanism. Asymmetry is indispensable since only then inter-eigenmodes fluxes of "energy" can be formed. Such an unbalanced flow renders particular eigenmodes overpopulated or amplified. Crucially, this mechanism does not break the overall stability – the dissipation eventually wins over and the amplification effect is only temporary or transient.

In the paper we focus on such transient trajectories for a system of non-linear ODE’s. On an elementary example we show how eigenvalue-based linearization technique becomes misleading and how simultaneously the transient property is developed. We utilize the May-Wigner model to include this mechanism and inspect its generic features. We find a new transient regime where transient trajectories are present although uncommon. Based on these findings, we identify relevant degrees of freedom driving the transient behaviour and propose a natural extension of the May-Wigner model. Such a modification leads to a generic transient behaviour arising as a robust transient phase.

.1 Stability of systems

We focus on a typical complex dynamical system of a set of $N$ first-order non-linear ordinary differential equations:

[TABLE]

where $x_{i}$ are the relevant degrees of freedom (neurons, concentrations of chemical compounds, species, etc.) and the non-linear functions $f_{i}$ encode the interactions (f.e. the Lotka-Volterra competitive predator-prey model for $f_{i}=x_{i}(1-\sum_{j}\alpha_{ij}x_{j})$).

If the form of functions $f_{i}$ is known, a question of stability is answered by a standard argument revised here briefly. In present analysis we ignore chaotic attractors or limit cycles and restrict to a simple binary notion of stability, the latter was addressed recently in Ipsen and Schomerus (2016). As a first step, we find all the points $f_{i}(x^{*})=0$ at which the solutions remain constant in time. Next, we expand Eq. (1) around a certain point $x^{*}$ from that set:

[TABLE]

and find a linearized system of equations

[TABLE]

where $M_{ij}=\partial_{x_{j}}f_{i}(x)_{|x=x^{*}}$. According to the Hartman-Grobman theorem (H-G theorem) Teschl (2012), the chosen point $x^{*}$ is stable if the real parts of the eigenvalues of $M$ are all strictly negative and unstable otherwise. The main assumption in the H-G theorem is that of locality – the perturbation $y$ around a stable point should be small. We address its importance in an example considered in the following and presented in Fig. 2.

To proceed, we define the norm of the solution $y_{i}$ of Eq. (3) as $|y(t)|^{2}=\sum_{i}|y_{i}(t)|^{2}$ and group them into into three groups:

stable non-transient (or non-transient) when $|y(t)|^{2}\overset{t\to\infty}{\longrightarrow}0$ and $\underset{t}{\max}|y(t)|^{2}=|y(0)|^{2}$, 2. 2.

stable transient (or transient) when $|y(t)|^{2}\overset{t\to\infty}{\longrightarrow}0$ and $\underset{t}{\max}|y(t)|^{2}\neq|y(0)|^{2}$, 3. 3.

unstable when $|y(t)|^{2}\overset{t\to\infty}{\longrightarrow}\infty$.

Instances of these types are shown in Fig. 1.

We present an example demonstrating the importance of locality assumption and simultaneously motivating this study. We define an $N=2$ dimensional non-linear system:

[TABLE]

which has two relevant stable points: $x^{*}=(x_{1}^{*},x_{2}^{*})=(0,0)$ and $x^{**}$ (given implicitly). We linearize the system around $x^{*}$ and find the matrix:

[TABLE]

so that Eq. (3) reads

[TABLE]

The resulting matrix $M_{|x^{*}}$ is in a triangular form, eigenvalues $-1,-2$ are strictly negative, the point $x^{*}$ is stable and so is the full system given by Eq. (4). We can solve the Eq. (6) explicitly and find the norm of its solution $|y(t)|^{2}=\left(y^{(0)}_{2}-\alpha y^{(0)}_{1}\right)^{2}e^{-4t}+2y^{(0)}_{1}\alpha e^{-3t}\left(y^{(0)}_{2}-\alpha y^{(0)}_{1}\right)+e^{-2t}\left(y^{(0)}_{1}\right)^{2}(1+\alpha^{2})$ depending on the initial value vector $y^{(0)}=\left(y^{(0)}_{1},y^{(0)}_{2}\right)$. From this formula one readily computes that for $\alpha>\frac{(y^{(0)}_{1})^{2}+2(y^{(0)}_{2})^{2}}{2y^{(0)}_{1}y^{(0)}_{2}}$, the transient behaviour of the norm $|y(t)|^{2}$ is present and absent otherwise.

In Fig. 2 we inspect the trajectories of both full and linearized system given by Eqs (4) and (6) as we vary the $\alpha$ parameter. We put emphasis on two features emerging in a correlated fashion – shrinkage of the $x^{*}$ related basin of attraction of the full system (top plots of Fig. 2) and simultaneous development of transient dynamics in the linearized system (bottom plot of Fig. 2). In the spirit of previous studies on the linearized dynamics Trefethen et al. (1993); Schmid (2007), we hypothesize that $\alpha$ parameter (representing all of the non-eigenvalue degrees of freedom for $N>2$) drives both mechanisms. Therefore, extracting transient behaviour becomes particularly important when either the non-linear solution or the structure of the phase space is not known and so only the linearized information is accessible. Then, transient dynamics can be seen as a herald of (non-linear) instability present already in the linear regime. This observation is closely related to the study of basins of attraction by the so-called Lyapunov functions.

.1.1 Randomness and complexity

Our aim is to study statistical features of transient phenomena highlighting its average features. To this end we chose a framework of random matrices as a unique insight into generic behaviour of such systems and is often treated as a first approximation or null-model analogous to a Gaussian distribution in univariate statistical analysis.

A matrix $M$ of size $N\times N$ introduced in Eq. (3) is taken to be

[TABLE]

where $\mu$ is understood as a diagonal matrix with entries equal to $\mu>0$. It is used since the linearization procedure is computed at a stable point by assumption. In the following we consider both real and complex matrices $X$ denoted by an index $\beta=1$ and $\beta=2$ respectively. The matrix $X$ is random and drawn from a joint pdf:

[TABLE]

where $\sigma^{2}$ is the variance, the real matrix is decomposed as $X_{kl}=x_{kl}$ whereas the complex matrix reads $X_{kl}=x_{kl}+iy_{kl}$. The joint measure for the real case reads $[dX]_{\beta=1}\equiv\prod_{i,j=1}^{N}dx_{ij}$, for the complex case is $[dX]_{\beta=2}\equiv\prod_{i,j=1}^{N}dx_{ij}y_{ij}$ and the normalization constant $c_{\beta}^{-1}=\int P_{\beta}(X)[dX]_{\beta}$. A notation for the average over $X$ is given by

[TABLE]

Lastly, we describe indicators of transient trajectories $y(t)$ as introduced in Neubert and Caswell (1997). By inspecting definitions of trajectories depicted in Fig. 1 one readily finds a good description of transient behaviour as the maximal possible amplification of the norm

[TABLE]

which identifies trajectory as stable non-transient if $A=1$, stable transient if $1<A<\infty$ and unstable if $A\to\infty$. Since for arbitrary systems it is hard to compute explicitly, instead a reactivity parameter $R$ was proposed:

[TABLE]

as a measure of the initial response of the system. By restricting to stable trajectories, we define transient behaviour if $R>0$ and lack thereof if $R<0$. The reactivity is an imperfect indicator – truly transient trajectories can be misidentified as a non-transient. Since the opposite cannot occur, it systematically over-counts non-transient trajectories and numerical results suggest it is a small effect.

We compute reactivity by writing down a formal solution to Eq. (3):

[TABLE]

where we introduce an identification between the vectors $y_{i}$ and kets $(\ket{y})_{i}$ and denote the initial vector as $\ket{y_{0}}=\ket{y(0)}$. We plug Eq. (12) into Eq. (11) and find:

[TABLE]

where the braket notation dictates that $\left<y_{0}|y_{0}\right>=|y(0)|^{2}$. The two measures are related as the reactivity $R$ is a linear term in the expansion of the amplification $A$ around $t=0$, $A=\max_{t}\left(1+Rt+\mathcal{O}(t^{2})\right)$. It therefore takes into account only the initial amplification.

Lastly, we address the treatment of initial conditions $\ket{y_{0}}$. A priori, we consider two scenarios – an extreme case where we chose a particular vector $y_{0}$ to maximize the quantity in question (f.e. reactivity) or a typical case where we average over all initial conditions. We introduce a notation to designate both scenarios:

[TABLE]

with a flat measure $[dy_{0}]_{\beta}$ over real $\beta=1$ or complex $\beta=2$ initial vectors and $p_{0}$ denoting a prescribed pdf for the initial conditions.

.2 Transient phase in the May-Wigner model

The model defined by Eqs (7) and (8) was introduced in the seminal work of May May (1972) to answer a key question in biological systems about the interplay between stability and complexity. The main finding is that there is an inherent (linear) instability of the system as we increase the complexity (matrix size $N$). Firstly, we review briefly this classic result and show how to extend the model to also include transient dynamics.

To recreate classic stability regimes we inspect the eigenvalue spectrum of $M$ as the matrix size grows to infinity $N\to\infty$. The asymptotic spectral density $\rho_{M}(x,y)=\lim\limits_{N\to\infty}\frac{1}{N}\left<\text{Tr}\delta(x+iy-M)\right>_{X}$ is given by the circular law Girko (1985):

[TABLE]

where $\theta$ is the Heaviside theta function confining the eigenvalues inside a circle of radius $\sigma$ centered around $(-\mu,0)$. In the large $N$ limit, the result (16) is valid for both values of $\beta=1,2$. The standard stability criterion based on the H-G theorem means that all eigenvalues of $M$ have real parts less than zero. In geometric terms, we keep the circular support of $\rho_{M}$ from crossing the $x=0$ line to stay in the stable regime. The stability/instability transition thus occurs along with the crossing and there are two equivalent ways of achieving that – by increasing the radius $\sigma$ or by moving the center point $\mu$. In further discussion we focus on the latter formulation and restrict to modifying the $\mu$ parameter. Thus, we identify two regimes – stable for $\mu>\mu_{S}$ and unstable for $\mu<\mu_{S}$ with $\mu_{S}=\sigma$ and depict this transition in Fig. 3.

By inspecting transient character of trajectories, stable regime is additionally split into transient and non-transient parts. A boundary is defined by the maximal reactivity of Eq. (11) with $\left<R_{\max}\right>_{X}>0$ for stable transient and $\left<R_{\max}\right>_{X}<0$ for stable non-transient regime. A similar boundary was studied in Tang and Allesina (2014) and in our context it is an example of an extreme scenario in the sense of Eq. (14). As reactivity of Eq. (11) is a Rayleigh quotient, it can be shown that $R_{\max}$ is given by the largest eigenvalue of $M^{\dagger}+M$:

[TABLE]

Because the matrix $M^{\dagger}+M$ is symmetric for $\beta=1$ (or Hermitian for $\beta=2$), its spectrum in the large $N$ limit is a translated Wigner’s semicircle $\rho_{M^{\dagger}+M}(\lambda)=\frac{1}{4\pi\sigma^{2}}\sqrt{8\sigma^{2}-(\lambda+2\mu)^{2}}$. Using this result we read out the rightmost edge and so the averaged maximal reactivity in the large $N$ limit reads

[TABLE]

where $\mu_{T}=\sqrt{2}\sigma$. We identify the stable transient regime for $\mu<\mu_{T}$ and stable non-transient regime for $\mu>\mu_{T}$ and presented it in Fig. 3.

Although we understand stable and unstable regimes quite well, it remains to inspect further the novel transient regime when $\mu\in(\mu_{S},\mu_{T})$. A natural question to ask is how abundant transient amplification is in the ensemble of trajectories. It is relevant since, by using a different criterion based instead on a reactivity averaged over the initial conditions $R_{\text{av}}$ given by Eq. (14), we find $\left<R_{\text{av}}\right>_{X}=-2\mu$. As $\mu>0$, $\left<R_{\text{av}}\right>_{X}$ is always negative and does not predict a transient behaviour.

.2.1 Density and abundance of transient trajectories

To inspect the question of abundance of transient trajectories, we define a probability density for the reactivity $g(r)=\delta(r-R(y_{0}))$ and consider both maximal and typical densities:

[TABLE]

Although the averaging over $X$ and $y_{0}$ is interchangeable $\overline{g_{\text{av}}}=\overline{g}_{\text{av}}$, the $\max$ operation and average is not $\overline{g_{\max}}\neq\overline{g}_{\max}$. Firstly we compute the density $\overline{g}$ averaged only over $X$:

[TABLE]

where an implicit dependence of $\overline{g}$ on $\beta$ is assumed. We use the delta function representation $\delta(x)=(2\pi)^{-1}\int dpe^{ipx}$, rewrite $ip(r-R)=ip(r+2\mu)-i\alpha\text{Tr}Y(X+X^{\dagger})$ where $\alpha=p(\text{Tr}Y^{2})^{-1}$ and set $Y=\ket{y_{0}}\bra{y_{0}}$. We compute the integral (19) by completing the square: $i\alpha\text{Tr}Y(X+X^{\dagger})+\frac{\beta N}{2\sigma^{2}}\text{Tr}X^{\dagger}X=\frac{\beta N}{2\sigma^{2}}\text{Tr}\left(X^{\dagger}+\frac{2i\sigma^{2}\alpha}{\beta N}Y\right)\left(X+\frac{2i\sigma^{2}\alpha}{\beta N}Y\right)+\frac{2\alpha^{2}\sigma^{2}}{\beta N}\text{Tr}Y^{2}$. The result is a quadratic Fourier integral:

[TABLE]

which no longer depends on the initial values $Y$ as $\text{Tr}Y^{2}=(\text{Tr}Y)^{2}$. The remaining integration gives a Gaussian distribution

[TABLE]

with mean $-2\mu$ and variance $\sigma_{\beta,R}^{2}=\frac{4\sigma^{2}}{\beta N}$. Eq. (21) is already the typical reactivity density $\overline{g_{\text{av}}}$ as it does not depend on the choice of initial conditions $y_{0}$, and so trivially $\overline{g}=\overline{g_{\text{av}}}$.

We turn to the extreme reactivity density which is related to the largest eigenvalue of $X^{\dagger}+X$:

[TABLE]

where an implicit dependence of $\overline{g_{\max}}$ on $\beta$ is assumed. In the literature on extreme value statistics Tracy and Widom (2009), one defines a cumulative distribution function $F_{N,\beta}(t)=\int[dH]e^{-\frac{\beta}{2}\text{Tr}H^{2}}\theta\left(t-\lambda_{\max}(H)\right)$ of the largest eigenvalue. By a simple rescaling, the extreme reactivity density therefore reads

[TABLE]

In this case, the order of operations is crucial – taking first the average over $X$ and then maximizing will reduce to the average scheme as $\overline{g}_{\max}=\overline{g}$. This is expected as the extreme scenario of any observable $O$ can be realized as an average scheme given by Eq. (14) but over a particular point-source density $\rho_{0}\sim\delta(y_{0}-v_{\max}(X))$ dependent on the maximal eigenvector $v_{\max}$ of $X$ itself. If so, the two averages no longer commute.

The abundances of transient trajectories are found as tail distributions of previously computed densities given in Eqs (21) and (22):

[TABLE]

which results in an implicit formula for the $\overline{N_{\max}}$ and an explicit one for $\overline{N_{\text{av}}}$:

[TABLE]

We compute the asymptotic forms of both abundances as $N\to\infty$. The abundance of a typical transient trajectory is asymptotically Gaussian:

[TABLE]

Although the abundance of an extreme transient trajectory is expressed in terms of the $F_{N,\beta}$ function which is not known in an explicit form, instead we cite two results valid in the $N\to\infty$ limit. To this end, we set $\mu=\mu_{T}-\Delta$ and inspect deviations around a typical value $\mu_{T}$ on different scales $\Delta$. The asymptotic formula of $\overline{N_{\max}}$ for large deviations $\Delta\sim\mathcal{O}(1)$ was found in Dean and Majumdar (2006, 2008) as:

[TABLE]

with $\phi(x)=\frac{1}{108}\left\{-x^{4}+36x^{2}+\sqrt{x^{2}+6}\left(x^{3}+15x\right)+\right.$ $\left.+27\left(\log(18)-2\log\left(\sqrt{x^{2}+6}-x\right)\right)\right\}$. We note the same formulas arose in a discussion of the symmetric May-Wigner model near its stability transition in Majumdar and Schehr (2014) which is however not equivalent to our case. For small perturbations $\Delta\sim\mathcal{O}(N^{-2/3})$ we cite another result:

[TABLE]

where $F_{\beta}$ is the Tracy-Widom distribution Tracy and Widom (1993, 1994). The abundance of extreme transient trajectories is therefore given as:

[TABLE]

We plot both $\overline{N_{\text{av}}}$ and $\overline{N_{\max}}$ given by Eqs (26) and (31) in Fig. 4 along with numerical results. The abundance of extreme transients $\overline{N_{\max}}$ is directly related to the transient/non-transient boundary at $\mu=\mu_{T}$ as it becomes a sharp theta function $\overline{N_{\max}}\sim\theta(\mu_{T}-\mu)$ in the (global) $N\to\infty$ limit. For large and intermediate values of $N$, the maximal abundance $\overline{N_{\max}}$ increases rapidly as we traverse the $\mu_{T}$ boundary and reaches unity upon entering the unstable regime near $\mu_{S}$. The nature of the abundance of typical transients $\overline{N_{\text{av}}}$ is different – it varies between $1/2$ for $\mu\to 0$ and [math] if $\mu\to\infty$ however reaches zero in the transient regime between $\mu_{S}$ and $\mu_{T}$ quite rapidly. Moreover, as the size of the matrix $N$ increases, $\overline{N_{\text{av}}}$ approaches zero for all values of $\mu$ according to Eq. (26) and does not result in a transition. Additionally, for any finite $N$ we find the average abundance being considerably smaller than the extreme one.

We recapitulate these two complementary viewpoints: 1) transient regime for intermediate parameters $\mu\in(\mu_{S},\mu_{T})$ is found as the abundance of extreme trajectories $\overline{N_{\max}}$ is close to unity. Moreover, it increases with the system size $N$ and reaches certainty for formally infinite systems. 2) On the other hand, according to the abundance of average trajectories $\overline{N_{\text{av}}}$, the number of typical transient trajectories is relatively small when the initial conditions are not especially tailored. In fact, their number on average decreases rapidly with the growth of the systems’ size as shown in Eq. (26).

Main conclusion is that although transient trajectories are (potentially) present in the whole transient regime $\mu\in(\mu_{S},\mu_{T})$ as shown by the behaviour of $\overline{N_{\max}}$, they are otherwise uncommon as dictated by $\overline{N_{\text{av}}}$.

.2.2 Generators of transient behaviour

Up to now we have considered the May-Wigner model of Eq. (3) with matrices drawn from Eq. (8) and found that when $\mu\in(\mu_{S},\mu_{T})$, transient trajectories are present although rare. Drawn by the interest of the transient behaviour itself, we ask a related question – although in May-Wigner model these trajectories are not found often, which features of a matrix we can tweak so that it produces generic transient trajectories. To put it differently, how to amplify the abundance $\overline{N_{\text{av}}}$ and simultaneously stay in the stable regime. To this end, we recall the definition of reactivity given in Eq. (11):

[TABLE]

where the presence or absence of transient behaviour was defined by the sign of $R$. If $X$ is drawn from Eq. (8), we have shown previously that although $\left<R_{\max}\right>_{X}$ can be positive (which sets the scale $\mu_{T}$), on average $\left<R\right>_{X}=-2\mu$ is always negative and thus no typical transient behaviour is present. To circumvent this we need to modify the matrix measure (8) accordingly.

Our aim is to render the reactivity (32) positive. For pure May-Wigner models defined by Eq. (8), it is always negative since both $\left<X\right>_{X}$ and $\left<X^{\dagger}\right>_{X}$ are zero. A simplest route of introducing a pdf with a non-zero mean $\left<X\right>^{\prime}_{X}=X_{0}\neq 0$ does not produce satisfactory result since then also the eigenvalue density of $X$ is modified resulting in an instability. The way out is to freeze the eigenvalues and tweak the remaining degrees of freedom. To achieve that, we introduce a Schur decomposition Horn and Johnson (1990):

[TABLE]

where the matrix $O$ is orthogonal ($\beta=1$) or unitary ($\beta=2$), $Z$ is a diagonal matrix with eigenvalues and $T$ is a strictly upper-triangular matrix encoding the non-orthogonality of the eigenbasis, hereafter we will refer to it as the non-orthogonality matrix.

For $\beta=1$, both $Z$ and $T$ have a block structure determined by the character of eigenvalues which are either purely real or form complex conjugate pairs. If we assume there are $k$ purely real eigenvalues and $\frac{N-k}{2}$ conjugate pairs of complex eigenvalues, the blocks of $Z$ and $T$ have dimensions $k\times k$ and $\frac{N-k}{2}\times\frac{N-k}{2}$ on the diagonal and $k\times\frac{N-k}{2}$ on the off-diagonal. The individual entries are either numbers or $2\times 2$ matrices on the diagonal and $2\times 1$ vectors on the off-diagonal.

For $\beta=2$, both $Z$ and $T$ have a simple structure - $Z$ is diagonal filled with eigenvalues and $T$ is strictly upper-triangular.

In the $\beta=1$ case, this decomposition induces a change of variables in the measure (8) computed in Lehmann and Sommers (1991); Edelman (1997) and for $\beta=2$ it is trivial. Because both 1) the Jacobian of the transformation (33) factorizes into $Z$ and $T$ dependent parts and 2) the Gaussian factor of Eq. (8)

[TABLE]

factorizes for both $\beta=1,2$, the eigenvalue matrix $Z$ and the non-orthogonality matrix $T$ fully decouple and can vary independently. Using Eq. (33), one finds that the averaged reactivity

[TABLE]

is likewise decoupled. The first term $\left<Z^{\dagger}+Z\right>_{X}$ is zero when $X$ is drawn from Eq. (8) as can be shown by a symmetry argument – for $\beta=1$ we reflect real eigenvalues $\lambda_{i}\to-\lambda_{i}$ along with the real parts of complex pairs $\text{Re}z_{i}\to-\text{Re}z_{i}$ and for $\beta=2$ we just set $z_{i}\to-z_{i}$. Although the second part $\left<T^{\dagger}+T\right>_{X}$ is also zero when averaged over Eq. (8), it does not need to be the case. In particular, we can fix $T$ by a constraint $\delta(T-T_{0})$ nad define a fixed $T_{0}$ May-Wigner model:

[TABLE]

where $P_{\beta}(X)$ was introduced in Eq. (8) and $c^{\prime}_{\beta}$ is the appropriate constant. When averaged over Eq. (36) denoted as $\left<\cdots\right>_{\tilde{P}_{\beta}}$, we find an average reactivity

[TABLE]

Due to independence of $Z$ and $T$, introducing a pdf in Eq. (36) does not change the spectrum of $X$. We define an external field/non-orthogonality parameter:

[TABLE]

so that for a fixed $T_{0}$ model given by Eq. (36), a generic transient behaviour is found when $\tau>2\mu$ and absent otherwise. The resulting phase diagram is shown in Fig. 5.

Now the parameter $\tau$ depends on the initial condition $y_{0}$. In particular, if the average $\left<\tau\right>_{y_{0}}$ is taken over the initial value vectors drawn from a symmetric density $p_{0}(y_{0})=p_{0}(-y_{0})$, it will vanish. A non-zero contribution is however produced if any asymmetry is present in $p_{0}$ or $y_{0}$ is held fixed. This was absent previously as the order parameter $\left<R\right>_{X}$ was completely decoupled from the initial condition $y_{0}$.

We interpret $T_{0}$ as a "field" conjugated to the order parameter $\left<R\right>_{X}$ akin to the the magnetic field and magnetization. Since the non-orthogonality matrix $T_{0}$ drives the transient phase transition, we consider the unperturbed system of $T_{0}=0$ in Eq. (36). It describes a normal matrix model, defined also by the condition $[X^{\dagger},X]=0$ and considered mostly when the matrix $X$ is complex Chau and Zaboronsky (1998); Wiegmann and Zabrodin (2003); Teodorescu et al. (2005). In this particular case, we compute an exact formula for the average norm:

[TABLE]

where we have used the formal solution given by Eq. (12) and $Y=\ket{y_{0}}\bra{y_{0}}$. Since $T_{0}=0$ we obtain $X=OZO^{\dagger}$ by the Schur decomposition and find

[TABLE]

where $P^{\prime}_{\beta}(Z)$ is the eigenvalue pdf for both values of $\beta=1,2$. The unitary/orthogonal integral is computed as $\int O_{ij}O^{\dagger}_{kl}[dO]_{\beta}=\frac{1}{N}\delta_{il}\delta_{jk}$ and the result reads:

[TABLE]

where $\tilde{\rho}_{X}(z)=N^{-1}\left<\sum_{i=1}^{N}\delta^{2}(z-z_{i})\right>_{P^{\prime}_{\beta}(Z)}$ is the spectral density of the normal matrix model.

In the large $N$ limit, the spectral density of a normal matrix $\tilde{\rho}_{X}$ also forms a circular law given in Eq. (16) with $\mu=0$. We plug it into Eq. (40) and find the average norm in the large $N$ limit as

[TABLE]

It is monotonically decreasing (as shown in the bottom inset of Fig. 5) and does not present transient behaviour in accordance with previous results.

.2.3 Non-orthogonality matrix and eigenvectors

The non-orthogonality matrix $T$ present in the Schur decomposition given by Eq. (33) is a crucial element in the development of transient behaviour. We will show in what sense the matrix $T$ is a measure of non-orthogonality and how it is related to other eigenvector related phenomena. In this section we restrict to the complex $\beta=2$ case. To this end, we diagonalize the matrix by a similarity transformation $S$:

[TABLE]

where $S^{-1}_{ij}$ is composed of left eigenvectors $\bra{L_{i}}_{j}$ and $S_{ji}$ of right eigenvectors $\ket{R_{i}}_{j}$:

[TABLE]

Left and right eigenvectors are bi-orthogonal $\braket{L_{i}}{R_{j}}=\delta_{ij}$ but not orthogonal in each space separately i.e. $\braket{L_{i}}{L_{j}}\neq\delta_{ij}$ and $\braket{R_{i}}{R_{j}}\neq\delta_{ij}$. In terms of the matrix $S$, these two relations are rewritten as $\braket{L_{i}}{L_{j}}=(S^{\dagger}S)_{ij}$ and $\braket{R_{i}}{R_{j}}=(S^{\dagger}S)^{-1}_{ij}$. A formal relation between $S$, $T$ and $Z$ is found by juxtaposing Eqs (33) and (42):

[TABLE]

If $T=0$ and eigenvalues are non-degenerate we find $O^{\dagger}S\sim 1$, $S$ becomes an unitary matrix with left and right eigenvectors rendered orthogonal. For non-zero $T$, one finds recurrence relations between the orthogonality matrices $S^{\dagger}S$, $Z$ and $T$ Chalker and Mehlig (1998). These orthogonality matrices are used in the definition of a 1-point eigenvector correlation function

[TABLE]

with weight $O_{ii}=(S^{\dagger}S)_{ii}(S^{\dagger}S)^{-1}_{ii}$. This object becomes a necessary ingredient in hydrodynamical description of complex dynamical matrices as established in Burda et al. (2015, 2014). The $O_{ii}$’s are also related to the eigenvalue condition number $\kappa(z_{i})$ as shown in Belinschi et al. (2017).

.3 Conclusions

In this paper we characterize transient phenomena in generic complex systems. Motivated by both physics and interdisciplinary applications, we argue that transient behaviour is complementary to the stability analysis and hints at non-linear features already on the linear level.

A seminal May-Wigner model is introduced as an example where we discuss transient dynamics. Based on the reactivity defined as an initial response of the system, we identify seminal unstable and stable regimes and find the latter to be additionally split into transient and non-transient regions. In the stable transient area, we compute the abundance (number) of transient trajectories in both extreme and typical choice of initial conditions. We conclude that although for certain special initial conditions trajectories do become transient, typically they do not show up. To amplify this typical transient behaviour, we introduce a modified May-Wigner model where only the non-orthogonality matrix $T_{0}$ is tweaked. This provides a model where a typical transient behaviour arises as a result of a non-zero value of $T_{0}$ with close relation to normal matrix models if $T_{0}$ vanishes. Lastly, the non-orthogonality matrix is discussed in relation to the eigenvector correlation function.

Transient dynamics have many faces – they are often described as a destabilizing mechanism Trefethen et al. (1993). In particular, this study points toward identifying crucial characteristics important especially in the context of early-warning signals of transitions in complex systems Scheffer et al. (2009). In the studies of neural networks, it is responsible for memory effects Ganguli et al. (2008). Paradoxically, in describing food-webs when the probed time span is relatively short wrt. the characteristic attenuation time of the transient, it can be reintrepreted as an effectively stable solution Hastings (2004).

This papers sheds light on the relevant characteristics of transient behaviour, enabling the tools needed to assess the severity of transient behaviour in the system at hand and their ultimate stability. In the spirit of recent work Fyodorov and Khoruzhenko (2016), in this work the phase diagram of May-Wigner models gets also refined to augment the stability questions. Additionally, transient features are robust when an additional structure is introduced into $X$ as presented in the Tang and Allesina (2014), pointing naturally into questions of universality. It is therefore an interesting point to study transient dynamics in general noise-plus-structure models of Grela and Guhr (2016) with special emphasis on the application to neuronal networks Ahmadian et al. (2015); Rajan and Abbott (2006).

Acknowledgements.

Author appreciates discussions with G. Schehr, S. Majumdar, M. A. Nowak and acknowledges the support of Grant DEC-2011/02/A/ST1/00119 of the National Centre of Science.