Resonances in a Chaotic Attractor Crisis of the Lorenz Flow

TL;DR

This paper investigates the behavior of Ruelle-Pollicott resonances in the Lorenz flow during a boundary crisis, revealing that stable resonances approach the imaginary axis while unstable ones do not, challenging existing early warning methods.

Contribution

It introduces a novel analysis of Ruelle-Pollicott resonances in the Lorenz system during a boundary crisis, highlighting their different behaviors and implications for crisis prediction.

Findings

Stable resonances approach the imaginary axis during the crisis.

Unstable resonances do not indicate proximity to the crisis.

Traditional early warning indicators may be ineffective near boundary crises.

Abstract

Local bifurcations of stationary points and limit cycles have successfully been characterized in terms of the critical exponents of these solutions. Lyapunov exponents and their associated covariant Lyapunov vectors have been proposed as tools for supporting the understanding of critical transitions in chaotic dynamical systems. However, it is in general not clear how the statistical properties of dynamical systems change across a boundary crisis during which a chaotic attractor collides with a saddle. This behavior is investigated here for a boundary crisis in the Lorenz flow, for which neither the Lyapunov exponents nor the covariant Lyapunov vectors provide a criterion for the crisis. Instead, the convergence of the time evolution of probability densities to the invariant measure, governed by the semigroup of transfer operators, is expected to slow down at the approach of the crisis. Such convergence is described by the eigenvalues of the generator of this semigroup, which can be divided into two families, referred to as the stable and unstable Ruelle-Pollicott resonances, respectively. The former describes the convergence of densities to the attractor (or escape from a repeller) and is estimated from many short time series sampling the phase space. The latter is responsible for the decay of correlations, or mixing, and can be estimated from a long times series, invoking ergodicity. It is found numerically for the Lorenz flow that the stable resonances do approach the imaginary axis during the crisis, as is indicative of the loss of global stability of the attractor. On the other hand, the unstable resonances, and a fortiori the decay of correlations, do not ag the proximity of the crisis, thus questioning the usual design of early warning indicators of boundary crises of chaotic attractors and the applicability of response theory close to such crises.

Equations

$$\{\begin{array}{ll}\dot{x}& =\sigma (y-x)\\ \dot{y}& =x(\rho -z)-y\\ \dot{z}& =xy-\beta z\end{array}(x,y,z)\in {\mathbb{R}}^3,$$

$$DF(x,y,z)=\left(\begin{array}{ccc}-\sigma & \sigma & 0\\ \rho -z& -1& -x\\ y& x& -\beta \end{array}\right),(x,y,z)\in {\mathbb{R}}^3.$$

$$\mathrm{div}F=\mathrm{trace}DF=-1-\sigma -\beta ,$$

$$f=\sum _{i=1}^n{f}_i\frac{{\chi }_i}{\eta ({\chi }_i)}\mathrm{such}\mathrm{that}{f}_i={⟨f,{\chi }_i⟩}_{\eta }\mathrm{and}{\int }_{X}fd\eta =\sum _{i=1}^n{f}_i.$$

$$({\mathcal{P}}_t^{\eta }f{)}_j={⟨{\mathcal{P}}_t^{\eta }f,{\chi }_j⟩}_{\eta }$$

$$[{\mathbf{P}}_t^{\eta }{]}_{ij}$$

$${\lambda }_k=\frac{1}{\tau }\log {\zeta }_k(\tau ).$$

$$[{\mathbf{P}}_{\tau }^m{]}_{ij}=\frac{\#\{x_0^{(i)}\in B_i,x_{\tau }^{(i)}\in B_j\}}{\#\{x_0^{(i)}\in B_i\}},$$

$$[{\mathbf{P}}_{\tau }^{\mu }{]}_{ij}=\frac{\frac{1}{T_{\mathrm{samp}}}{\sum }_{s=0}^{T_{\mathrm{samp}}}{\chi }_i(x_{t_s}){\chi }_j(x_{t_s+\tau })}{\frac{1}{T_{\mathrm{samp}}}{\sum }_0^{T_{\mathrm{samp}}}{\chi }_i(x_{t_s})}=\frac{\#\{x_{t_s}\in B_i,x_{t_s+\tau }\in B_j\}}{\#\{x_{t_s}\in B_i\}}.$$

$$\dot{x}(t)$$

$${\mathcal{U}}_t^{\eta }:g\to g\circ {\Phi }_t,t\ge 0,$$

$${\int }_{X}g(x){\mathcal{P}}_t^{\eta }f(x)\eta (dx)={\int }_X{\mathcal{U}}_t^{\eta }g(x)f(x)\eta (dx)$$

$${\int }_A{\mathcal{P}}_t^{\eta }f(x)\eta (dx)={\int }_{{\Phi }_t^{-1}A}f(x)\eta (dx),t\ge 0.$$

$$\mu ({\Phi }_t^{-1}A)=\mu (A),\mathrm{for}\mathrm{any}A\in \mathcal{B}.$$

$${\int }_{X}g({\Phi }_{t}x)\mu (dx)={\int }_{X}g(x)\mu (dx):={⟨g⟩}_{\mu },$$

$${⟨g⟩}_{\mu }=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^{T}g({\Phi }_{t}x)dt:=\bar{g},\mathrm{for}\mu -\mathrm{almost}\mathrm{every}x.$$

$$C_{f,g}(t):={\int }_{X}f(x)g({\Phi }_{t}x)\mu (dx)-{⟨f⟩}_{\mu }{⟨g⟩}_{\mu },t\ge 0,$$

$$\underset{t\to \infty }{\lim }\mu (A\cap {\Phi }_t^{-1}B)=\mu (A)\mu (B),\mathrm{for}\mathrm{any}A,B\in \mathcal{B},$$

$$C_{f,g}(t)=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_0^T(f(x)-\bar{f})(g({\Phi }_{t}x)-\bar{g})dt,$$

$${\mathcal{U}}_t^X{1}_X=1_{{\varphi }_t^{-1}X\cap X}=1_X,t\ge 0.$$

$$W_s(x)=\{y\in X:d({\Phi }_t(x),{\Phi }_t(y))\to 0,t\to \infty \}.$$

$${\psi }^X(y)={\psi }^{\Lambda }(x),\mathrm{whenever}y\in W_s(x)$$

$${\psi }^X(y_n)={\psi }^{\Lambda }(x_n)$$

$$\underset{n\to \infty }{\lim }{\psi }^X(y_n)=\underset{n\to \infty }{\lim }{\psi }^{\Lambda }(x_n)={\psi }^{\Lambda }(x)={\psi }^{\Lambda }(y),$$

$${\mathcal{U}}_t^X{\psi }^X(y)={\psi }^X({\Phi }_{t}y)$$

$$={\psi }^{\Lambda }(x)\mathrm{whenever}y\in W_s({\Phi }_t^{-1}x)$$

$$={\psi }^{\Lambda }({\Phi }_{t}z)\mathrm{whenever}y\in W_s(z)$$

$$=\zeta {\psi }^{\Lambda }(z)\mathrm{whenever}y\in W_s(z)=\zeta {\psi }^X(y).$$

Full text

This behavior is investigated here for a boundary crisis in the Lorenz flow, for which neither the Lyapunov exponents nor the covariant Lyapunov vectors provide a criterion for the crisis. Instead, the convergence of the time evolution of probability densities to the invariant measure, governed by the semigroup of transfer operators, is expected to slow down at the approach of the crisis. Such convergence is described by the eigenvalues of the generator of this semigroup, which can be divided into two families, referred to as the stable and unstable Ruelle-Pollicott resonances, respectively. The former describes the convergence of densities to the attractor (or escape from a repeller) and is estimated from many short time series sampling the phase space. The latter is responsible for the decay of correlations, or mixing, and can be estimated from a long times series, invoking ergodicity.

It is found numerically for the Lorenz flow that the stable resonances do approach the imaginary axis during the crisis, as is indicative of the loss of global stability of the attractor. On the other hand, the unstable resonances, and a fortiori the decay of correlations, do not flag the proximity of the crisis, thus questioning the usual design of early warning indicators of boundary crises of chaotic attractors and the applicability of response theory close to such crises.

1 Introduction

It is a problem of fundamental relevance in mathematical, natural, and applied sciences to understand under which conditions a system may undergo abrupt changes under perturbation and, if so, predict when these changes will occur. Much of our understanding of such transitions comes from the bifurcation theory of autonomous dynamical systems Guckenheimer1983 ; Ruelle1989a ; Kuznetsov1998 , with extensions to nonautonomous Kloeden2011 and random Arnold2003 dynamical systems. In particular, local bifurcations, taking place for example when a stationary point or a limit cycle loses stability, are characterized by the critical exponents of these invariant sets. They yield a local measure of the relaxation rate of trajectories to these sets. As the latter become less stable, these exponents approach zero, resulting in the slowing down of the convergence of trajectories to the attractor. This critical slowing down has allowed to design early-warning signals of critical transitions by monitoring the rate of decay of correlations Held2004a , peaks in power spectra Kleinen2003a , or of recovery from perturbations Nes2007 . See Scheffer2009 , for a review, and Lenton2011a for applications to climate science.

Most physical systems of interest are, however, chaotic, in the sense that they support an invariant measure with positive Lyapunov exponents Eckmann1985 ; Young2002 . Such sets can be involved in global bifurcations during which they may lose their attracting character or undergo topological transformations (Ruelle1989a, , Chap. 17). A particularly important type of bifurcation is the boundary crisis Grebogi1983 , during which a chaotic invariant set ceases to be attracting due to, loosely speaking, its collision with a saddle 111A saddle is defined here as a possibly chaotic invariant set, which is attracting neither in forward nor in backward time (see Ashwin1996 for a more rigorous definition)., the so-called edge state. Most trajectories then undergo a transient before they converge to another attractor, if any. Such crises have been identified in the Hénon map and the Lorenz system Grebogi1983 , and are also found in high-dimensional turbulent flows Skufca2006 ; Schneider2007 ; Eckhardt2008 and climate models Bodai2015 ; Lucarini2017 .

Much less is known regarding changes in the statistical properties of these systems, and, in particular, if critical slowing down may be observed. While for chaotic attractors the notion of critical exponents can be generalized to that of Lyapunov exponents characterizing the stability of chaotic trajectories Oseledets1968 ; Eckmann1985 , the latter do not in general allow to infer to what extent a chaotic set is attracting. The Lyapunov spectrum is indeed calculated from a linearization about orbits on the invariant set and thus only provides local information. On the other hand, the size of the basin of attraction leads to the notion of global stability of an attractor (see e.g. LaSalle1976 ). For example, during the boundary crisis of an attractor, the size of the basin of attraction shrinks as the criticality is neared and vanishes when the set becomes unstable. It has been proposed in Faranda2014e to approach the problem studying large fluctuations of the systems using extreme value theory and relating the change in the qualitative properties of the extremes to the approach to the critical transition (see also Faranda2014b ).

From a statistical physics point of view, the divergence of nearby trajectories characterized by the positive Lyapunov exponents is a manifestation of chaos at a microscopic level. Macroscopically, the decay with time of correlations associated with the loss of memory to initial conditions of ensembles (as they mix) is a clear manifestation of chaos. The evolution in time of the correlation function between any appropriate observables is fully determined by the semigroup of transfer operators $\mathcal{P}_{t}^{\mu},t\geq 0$ Lasota1994 , governing the evolution of densities with respect to an invariant measure $\mu$ (e.g. supported by an attractor).

It is a classical result from ergodic theory Halmos1956 ; Arnold1968 that correlations between observables vanish for time lags going to infinity only if there are no eigenvalues of the transfer operator in the unit disk other than the eigenvalue 1. A more difficult problem, which is still a matter of investigation Baladi2001 ; Young2013 , is to characterize the rate of mixing. This rate depends on the position in the complex plane of the poles of the correlation spectrum of a pair of observables, which correspond to the Ruelle-Pollicott resonances Pollicott1985 ; Ruelle1986 . The latter are given by the eigenvalues with nonzero real part of the generator of the semigroup of transfer operators acting on anisotropic Banach spaces adapted to the dynamics of contraction and expansion of chaotic systems Liverani1995a ; Blank2001a ; Gouezel2006 ; Butterley2007 ; Faure2014 ; Baladi2017 . Note that, while both the presence of positive Lyapunov exponents and mixing are a manifestation of chaos, their relationship is nontrivial Collet2004 ; Alves2004 ; Pires2011 ; Slipantschuk2013 . In the case of uniformly hyperbolic systems, correlations are expected to undergo an initial fast decay associated with large Lyapunov exponents, but have an asymptotic decay bounded above by the smallest positive Lyapunov exponent Collet2004 . On the other hand, nonuniformly hyperbolic systems may have arbitrarily slow mixing rates.

While transfer operators acting on densities with respect to the invariant measure allow to study the ergodic and mixing properties, global information about the invariant set supporting this measure (e.g. an attractor) should be studied from transfer operators $\mathcal{P}_{t}^{m},t\geq 0$, acting on densities with respect to the Lebesgue measure (i.e. in phase space). This has recently led to new developments in the theory and applications of the stability of dynamical systems Vaidya2008 ; Mauroy2016 . In this case, the spectrum of the generator of the transfer semigroup, does not only capture the rate of mixing, but also the rate of convergence (escape) of densities to (from) an invariant measure supported by an attractor (repeller). It follows that, as an attractor becomes less attracting at the approach of a crisis, densities are expected to take more time to converge to the invariant measure, resulting in the slowing down of the decay of correlations, i.e. in the mixing. This slowing down should thus be associated with a decrease of the spectral gap in the spectrum of the transfer semigroup. This characterization of the stability of chaotic attractors has been used in Tantet2015a to give numerical evidences that, in a high-dimensional climate model undergoing a boundary crisis, the spectral gap in some eigenvalues of the transfer operators in heavily projected spaces indeed shrinks, explaining the slowing down of the decay of correlations observed from simulations.

It is, however, still unclear whether this change in the spectral properties of the transfer semigroup is a generic property of dynamical systems undergoing an attractor crisis. Moreover, the question emerges on whether such changes can be detected from time series on the attractor alone - which can be investigated looking at the properties of $P_{t}^{\mu}$ - or if perturbations of the system away from the attractor are needed, so that the eigenvalues of $P_{t}^{m}$ are the key objects of interest. We refer to the resonances that belong to the spectrum of $P_{t}^{m}$ but not of $P_{t}^{\mu}$ as the stable resonances, which are associated to the nearing of the orbits towards the attractor, while we use the expression unstable resonances for the eigenvalues of $P_{t}^{\mu}$.

This issue is directly related to the generalization of the Fluctuation-Dissipation Theorem (FDT) to nonequilibrium systems by Ruelle2009 (see also Cessac2007 ). In this regard, the problem is that in the case of deterministic systems the natural fluctuations explore only the unstable manifold of the attractor, whereas generally an external perturbation will impact both the unstable and stable directions. A specific example in the context of a geophysical system is considered in Gritsun2017 . On the other hand, for systems with an invariant measure which is absolutely continuous with respect to the Lebesgue measure, such as stochastically perturbed dynamical systems leading to hypoelliptic diffusions Hairer2010 , both transfer semigroups can be identified and the FDT is expected to hold. Note that a theory of the linear response to noise of the statistics of the system has been developed by Lucarini2012c , while the changes in the spectrum of the transfer semigroup due to noise have been studied by Gaspard2002a and by Tantet2016 , for the particular case of the Hopf bifurcation.

In this study, a chaotic boundary crisis in the Lorenz flow is analyzed in terms of eigenvalues of the transfer semigroup. Although this appears are a case study, we think that it can be of more general interest by revealing possible statistical properties of more general chaotic dynamical systems. The Lorenz flow is system of three ordinary differential equations derived in the seminal paper Lorenz1963a via a spectral truncation of the fluid equations for Rayleigh-Bénard convection. As one of the first examples of low-dimensional flows with robust chaos and because it exhibits a wide range of dynamical phenomena, it has been extensively studied (see Sparrow1982 for a review). Most of the dynamical properties of Lorenz-like attractors have been obtained for the geometric Lorenz system introduced in Guckenheimer1979 . However, these properties were later found to persist for general singular-hyperbolic attractors, yielding a paradigm for robust low-dimensional chaotic systems Araujo . In particular, the standard Lorenz attractor was proved to be a singular-hyperbolic chaotic attractor in Tucker1999 . As a consequence, the Lorenz attractor supports a unique physical measure with non-zero Lyapunov exponents (Araujo, , Chap. 7). Numerical evidences suggest that the standard Lorenz flow obeys linear Reick2002a and nonlinear response theory Lucarini2009b . However, whether the Lorenz flow is mixing and has exponential decay of correlations is still an open problem (Araujo, , Chap. 10.2). In addition, the Lorenz flow is known Grebogi1983 to undergo a boundary crisis for some parameter values, while its low-dimensionality renders numerical applications tractable. For these reasons, the Lorenz flow constitutes an interesting test bed to better understand the relationship between the eigenvalues of the transfer semigroup and the dynamics of the boundary crisis.

The analysis presented here is numerical and descriptive in nature, but should allow to build intuition for further studies. The Lorenz flow and the bifurcations of interest for this study are briefly introduced in section 2. The results presented in this section are mainly a reproduction of the ones presented in Sparrow1982 . However, we also discuss the fact that that neither the Lyapunov exponents nor the alignment of the covariant Lyapunov vectors, can give an indication of the approach of the crisis. This provides evidence that, contrary to what was found for other types of crises Rollins1984 ; Pompe1988 ; Mehra1996 ; Beims2016 , these important dynamical quantities are in general not useful to flag vicinity to a boundary crises. In section 3, two different methods are presented to approximate the transfer semigroups $P_{t}^{m}$ and $P_{t}^{\mu}$ to yield the stable and unstable resonances, respectively. The main results are presented in section 4, for different parameter values about the crisis, allowing to monitor the changes in the resonances during the crisis and to analyze whether these changes can be observed from the dynamics on the attractor alone. A summary and discussion is given on the implications of these results regarding the possibility of designing early-warning indicators of chaotic attractor crises, the effect of the addition of noisy perturbations as well as the viability of response theory close to the crisis. In Appendix A, we recapitulate some concepts of ergodic theory of dynamical systems that we deem relevant for the interpretation of our results ; The expert reader might want to skip this material. Appendix B gives evidence that the numerical results of section 4 are robust to the resolution, sampling and transition time.

2 Attractor crisis in the Lorenz system

Here, we summarize some general properties of the Lorenz flow, as well as the bifurcations that will be important for the rest of this study.

2.1 Dissipativity and boundedness

The Lorenz flow $\Phi_{t}$, with time $t$ in $\mathbb{R}$, is generated by the following set of Ordinary Differential Equations (ODEs),

[TABLE]

where $\rho,\sigma$ and $\beta$ are positive parameters and the dot indicates differentiation with respect to time. In this study, the parameters $\sigma$ and $\beta$ are set to classical values of $10$ and $8/3$, respectively. On the other hand, the Rayleigh number $\rho$ will be varied from [math], for which all trajectories converge to the stationary point at the origin, to the classical value of $28$, at which Lorenz obtained the celebrated ”butterfly” attractor.

The ODE (2.1) is invariant with respect to the change of variable $(x,y,z)\to(-x,-y,z)$, so that to each solution corresponds another symmetrically related one. The Jacobian of the vector field $F$ of the right-hand side of (2.1) is given by,

[TABLE]

The contraction rate of volumes under $F$ is given by its divergence,

[TABLE]

It is constant and, for the values of $\sigma$ and $\beta$ considered here, negative so that volumes contract uniformly in phase space. The vector field $F$ is therefore said to be dissipative.

From the quadratic nature of the equations (2.1), it was shown in Lorenz1963a (see also Lorenz1979b ) that solutions initiated in a specific ellipsoid remain there forever. In particular, for $\sigma=10$ and $\beta=8/3$, trajectories remain in the ball $R_{o}$ bounded by the sphere $S_{o}$ of radius $\rho+\sigma$ and center $(0,0,\sigma+\rho)$. This probability together with that of positive volume contraction ensures that trajectories initialized inside $R_{o}$ remain in $R_{o}$ and must converge to a non-wandering set of Lebesgue measure 0. The sphere $S_{o}$ is used in sections 3 and 4 to bound the domain on which the transfer operators are approximated.

2.2 Route to chaos in the Lorenz flow

In this section, we briefly describe the series of local and global bifurcations leading the chaotic attractor crisis of interest (e.g. Guckenheimer1983 ; Kuznetsov1998 for references on bifurcation theory). These results are summarized in the bifurcation diagram in figure 1. They are not new and should rather be considered as a reproduction of those discussed in Sparrow1982 . However, they will be useful to interpret later results discussed in section 4.

We start by increasing the control parameter $\rho$ from zero. For $\rho=0$, the stationary point $O$ at the origin is stable. Solving the algebraic equation $F(x)=0$ analytically reveals that two symmetric stationary points $p^{+}$ and $p^{-}$ other than $O$ exist for $\rho>\rho_{P}=1$. The eigenvalues of the Jacobian $DF$ evaluated at each point reveals that $O$ loses stability at $\rho_{P}$ while the points $p^{\pm}$ are stable. A pitchfork bifurcation thus occurs at $\rho_{P}$, as is illustrated in the left of the diagram in figure 1, on which the branch of the stationary point $O$ is represented in black, while that for $p^{+}$ is represented in red.

It is then numerically found that the Jacobian evaluated at $p^{\pm}$ quickly acquires a complex conjugate pair of eigenvalues. The latter eventually cross the imaginary axis at $\rho_{\mathrm{Hopf}}\approx 24.74$. This loss of stability of the two stationary points is associated with two symmetric subcritical Hopf bifurcations. That is, two symmetric unstable periodic orbits, denoted $\Gamma^{+}$ and $\Gamma^{-}$, merge with the stationary points $p^{+}$ and $p^{-}$, respectively. The branch of the orbit $\Gamma^{+}$ is represented in green in the diagram in figure 1.

Starting from $\rho=\rho_{\mathrm{Hopf}}$, we now decrease $\rho$ to continue the unstable periodic orbits $\Gamma^{\pm}$ as they emerge from the stationary points $p^{\pm}$. The two orbits $\Gamma^{\pm}$ are continued numerically via pseudo-arclength continuation (see e.g. (Kuznetsov1998, , Chap. 10)), yielding a numerical integration of the orbit, its period $T$, as well as the fundamental matrix $D\Phi_{T}$ at time $T$. The latter allows to calculate the Floquet exponents characterizing the stability of the limit cycle Guckenheimer1983 as the complex logarithm of the eigenvalues of $D\Phi_{T}$ divided by $T$ Hartman1964 . The results of this continuation are represented in figure 2. On the upper left panel, one can see that the period of the limit cycles is increasing rapidly as $\rho$ is decreased below the Hopf bifurcation. The real parts of the Floquet exponents are represented in the lower left panel. The Floquet exponent zero is associated with the direction of the flow, while one of the exponents is positive and the last is negative, indicating that the periodic orbits $\Gamma^{\pm}$ are unstable and of saddle type. As expected from normal form theory Guckenheimer1983 ; Kuznetsov1998 , the instability of the periodic orbits grows away from the Hopf bifurcation, i.e. as $\rho$ is decreased from $\rho_{\mathrm{Hopf}}$ and their radius increases. This can be seen on the right panel of the same figure 2, where the periodic orbits for different values of $\rho$ are represented in the $(x+y,z)$ plane.

The increase in the period of the periodic orbits $\Gamma^{\pm}$ and the decrease of their distance to the stationary point $O$ at the origin give numerical evidence that the periodic orbits eventually collide with $O$. It is thus expected that, at a particular value $\rho=\rho_{\mathrm{Homo}}$, the union of each periodic orbit with $O$ yields a homoclinic orbit on which $O$ is an asymptotic point for both forward and backward trajectories. The numerical method used here does not permit the continuation of the periodic orbits all the way to the homoclinic bifurcation. However, direct numerical integration of a trajectory starting in the neighborhood of the stationary point allows to get an approximation of the homoclinic orbit. The latter is obtained at $\rho=13.93\approx\rho_{\mathrm{Homo}}$, in agreement with Sparrow1982 . It is plotted by a thin dashed black line on the right panel of figure 2 and is also represented by the green ellipse in the diagram in figure 1. This homoclinic bifurcation has been studied in the geometric Lorenz system Guckenheimer1979 ; Guckenheimer1983 where a nontrivial invariant set $\Lambda$, with the structure of a singular horseshoe (Araujo, , Chap. 3), is known to emerge from its unfolding, i.e. for $\rho>\rho_{\mathrm{Homo}}$. An infinity of unstable periodic orbits are embedded in this horseshoe, including the two simple periodic orbits $\Gamma^{\pm}$. This set is thus responsible for transient chaos Kaplan1979a , but is not yet attracting for $\rho$ close to $\rho_{\mathrm{Homo}}$.

2.3 Attractor crisis

The bifurcation of interest in this study occurs as $\rho$ is increased from $\rho_{\mathrm{Homo}}$, when the chaotic invariant set $\Lambda$ becomes an attractor 222This attractor must be only partially hyperbolic, since it contains the singularity at the origin, which prevents the stable and unstable manifolds to be continuous.. The study of the one-dimensional Lorenz map Lorenz1963a ; Sparrow1982 ; Guckenheimer1983 , gives evidence that this crisis occurs for $\rho=\rho_{A}\approx 24.06$. That is, for $\rho_{\mathrm{Homo}}<\rho<\rho_{A}$, the chaotic set $\Lambda$ is of saddle type, while for $\rho>\rho_{A}$ and $\rho$ not too large, $\Lambda$ is an attractor. The branch of the chaotic set $\Lambda$ is represented in cyan in the diagram in figure 1, with $\Lambda$ represented by a lemniscate. Importantly, for $\rho_{\mathrm{Homo}}<\rho<\rho_{\mathrm{Hopf}}$, the strange attractor $\Lambda$ coexists with the two unstable periodic orbits $\Gamma^{\pm}$ as well as with two other attractors, the stable stationary points $p^{\pm}$. Following Sparrow1982 , an informal phenomenological description of this crisis is supported here by numerical integration.

Figure 3.(a) represents the minimal Euclidean distance between an aperiodic orbit of the chaotic set $\Lambda$ and the periodic orbits $\Gamma^{\pm}$, versus $\rho$. One can see that, close to the crisis, this numerical calculation of the distance is close to zero. This is more obvious from figure 3.(b), in which a long integration of an aperiodic orbit is represented in blue together with the periodic orbits $\Gamma^{\pm}$, in black, for $\rho\approx\rho_{A}$. In fact, below the crisis, for $\rho<\rho_{A}$, the periodic orbits $\Gamma^{\pm}$ are embedded in the chaotic saddle $\Lambda$ Sparrow1982 . On the other hand, this distance increases as $\rho$ is increased from $\rho_{A}$ due to the shrinkage of $\Gamma^{\pm}$ as the Hopf bifurcation is approached. The chaotic set $\Lambda$ thus appears to become on attractor only once the periodic orbits $\Gamma^{\pm}$ have left it. Since all the orbits in $\Lambda$ for $\rho<\rho_{A}$ cannot be mapped to those of $\Lambda$ for $\rho>\rho_{A}$, the flow on this set is not structurally stable at the boundary crisis; there is a genuine global bifurcation at $\rho=\rho_{A}$.

This bifurcation is better understood from the geometry of the stable manifold of the stationary point $O$. Numerical integrations of one branch of the unstable manifold of $O$ are represented in blue in figure 4, for different values of $\rho$ about the crisis. The periodic orbit $\Gamma^{-}$ is also represented by a thick black line, together with orbits belonging to the stable manifold of $\Gamma^{-}$, represented by thin black dashed lines. These numerical results suggest that:

Between the homoclinic bifurcation and the attractor crisis (for $\rho_{\mathrm{Homo}}<\rho<\rho_{A}$, Fig. 4.(a)), the unstable manifold of $O$ is connected to the stable manifold of the stable stationary point $p^{-}$ emanating from the pitchfork bifurcation, i.e. orbits repelled by $O$ converge to $p^{-}$, 2. 2.

At the crisis ($\rho=\rho_{A}$, Fig. 4.(b)), as a consequence of the shrinkage of the periodic orbit $\Gamma^{-}$ with increasing $\rho$, , the trajectories along the unstable manifold of $O$ no longer convergences to $p^{-}$, but instead connect to the stable manifold of the periodic orbit $\Gamma^{-}$, 3. 3.

After the crisis, ($\rho>\rho_{A}$ Fig. 4.(c)), the small size of $\Gamma^{-}$ prevents the unstable manifold of $O$ to dive into the stable manifold of $\Gamma^{-}$. These trajectories having to converge to an attractor (see Sect.  section 2.1), must wander along the chaotic set $\Lambda$, which is then attracting.

2.4 Lyapunov exponents and covariant Lyapunov vectors

Local bifurcations of stationary points or limit cycles are fully characterized by the crossing of the imaginary axis by the characteristic exponents of these invariant sets, i.e. the eigenvalues of the Jacobian for the stationary points, the Floquet exponents for the limit cycles Guckenheimer1983 ; Kuznetsov1998 . In these cases a bifurcation point corresponds to a loss of hyperbolicity due to the crossing of the imaginary axis by one or several critical exponents.

The notion of hyperbolicity can be generalized to nontrivial sets. A very strong form of hyperbolicity is that of uniform hyperbolicity for which the tangent space can be split into a one-dimensional space in the direction of the flow and contracting and expanding spaces, the stable and unstable spaces, respectively, with uniform decay and growth bounds, respectively (Katok1996, , Part 4).

A weaker form of hyperbolicity can be introduced making reference to the Lyapunov exponents. Lyapunov exponents describing the growth rate of perturbations applied to a trajectory provide a generalization of the characteristic exponents to chaotic invariant sets Oseledets1968 ; Eckmann1985 . Nonuniform hyperbolic systems do not obey uniform hyperbolicity but admit an invariant measure with nonzero Lyapunov exponents almost everywhere, except for the exponent associated with the direction of the flow Barreira2002 . In this case, the covariant vectors in the direction of which the growth of perturbations is given by the respective Lyapunov exponents, namely the Covariant Lyapunov Vectors (CLVs, e.g. Kuptsov2012 ). The CLVs associated with positive (neutral, negative) Lyapunov exponents are tangent to the unstable (center, stable) manifold along an orbit and span the unstable (center, stable) tangent space. In general, the angles between the manifolds can be arbitrarily small, as opposed to the uniform hyperbolic case.

The set $\Lambda$ in the Lorenz flow is chaotic, so that it has one positive Lyapunov exponent $\lambda_{+}$, one zero exponent $\lambda_{0}$ corresponding to the direction of the flow and one negative exponent $\lambda_{-}$, which is larger in absolute value than $\lambda_{+}$, in order for volumes to contract. Moreover, the stationary point at the origin belongs to $\Lambda$. These results were first proved for the geometric Lorenz attractor Guckenheimer1979 and later for the Lorenz flow for the classical parameter values $\rho=28,\sigma=10$ and $\beta=8/3$ Tucker1999 . The Lorenz attractor is thus at most partially hyperbolic Pesin2004 . In fact, this type of partially hyperbolic flow with a singularity in the attractor happens to be the prototype of robust chaotic flow in three dimensions Araujo and has been called singular hyperbolic. In the case of the Lorenz flow, the tangent space splits into a stable space, with contraction, and a center-unstable space where, loosely speaking, weak contraction is allowed in addition to expansion. This implies that the angle between the center-unstable space and the stable space is uniformly bounded away from zero Pesin2004 , and so are the angles between the stable CLV and each of the unstable and central CLVs. However, due to the stationary point at the origin, where the center manifold does not exist, the angle between the unstable and the central CLVs may be arbitrarily small as this singularity is approached.

Decreasing $\rho$ towards the attractor crisis at $\rho_{A}$, one may wonder whether the positive Lyapunov exponent becomes negative. However, since for $\rho<\rho_{A}$ a chaotic invariant set persists, the positive Lyapunov exponent must remain positive. Indeed, the negative value of the average volume contraction rate implies that the Lyapunov exponents different from zero must be either both negative or, if one of them is positive, the other must be negative and larger in absolute value. However, if both non-zero Lyapunov exponents of the invariant set $\Lambda$ were negative for $\rho<\rho_{A}$, this set would not be chaotic, leading to a contradiction. This heuristic argument is supported by numerical estimations 333Following the algorithm described in Kuptsov2012 , also yielding the CLVs. of the Lyapunov exponents of $\Lambda$ for varying values of $\rho$, represented in figure 5. One can see on this figure that the Lyapunov exponents vary only slightly close to the crisis and that none of them cross the imaginary axis.

In addition, since a break of hyperbolicity may be due to a homoclinic tangency where the stable and unstable manifolds are tangent, the angle between the CLVs has been suggested in Ginelli2007 and Yang2009 as an indicator of the degree of hyperbolicity of the dynamics. Moreover, it has been conjectured in Beims2016 that the alignment of the CLVs could provide a criterion to predict crises, although the latter are understood there as chaotic bursts (in other words, an extreme fluctuation in an otherwise weakly chaotic trajectory) rather than as an attractor crisis. However, due to the singularity in the Lorenz attractor, the angle between the unstable and the central CLVs can be arbitrarily small even away from the crisis, so that the loss of hyperbolicity should not provide a precursor of the crisis. This is confirmed from the computations of the CLVs for varying $\rho$, as well as for the angle between the unstable and the stable CLVs, and between the center and the stable CLVs (not shown here).

3 Approximation of the Ruelle-Pollicott resonances

The Ruelle-Pollicott resonances are given by the eigenvalues of the generator of the semigroup of transfer operators, governing the evolution of densities. While analytical results on the properties of this semigroup for chaotic systems are difficult to obtain, the relationship between the ergodic properties of dynamical systems and the spectrum of the semigroup of transfer operators is well known. A summary of this is given in appendix A. Here, we explain how these operators can be approximated from transition matrices estimated from time series. While the method is well known, we make clear the difference between two variants allowing either for the approximation of both the stable and the unstable resonances from the transfer operators $\mathcal{P}_{t}^{m}$ with respect to the Lebesgue measure, or only for the approximation of the unstable resonances from the transfer operators $\mathcal{P}_{t}^{\mu}$ with respect to the invariant measure supported by a chaotic attractor.

Calculating the Ruelle-Pollicott resonances is a difficult task. Explicit formulas have only been obtained for low-dimensional systems from a trace formula or from the decomposition of the eigenvectors on basis functions Hasegawa1992b ; Gaspard1992 ; Gaspard1995 ; Gaspard2001a . For high-dimensional systems, to our knowledge, only qualitative results limited to uniformly hyperbolic systems have been obtained regarding the distribution of the resonances in the complex plane Gouezel2006 ; Butterley2007 ; Faure2014 .

Numerical methods are thus required to approximate the Ruelle-Pollicott resonances and the associated eigenvectors. For still relatively low-dimensional systems, but with chaotic dynamics, one approach is to calculate the spectrum of transition matrices resulting from the projection of the transfer operators on a finite family of basis functions. The transition probabilities can then be estimated from time series. This Galerkin truncation with estimation is referred to as Ulam’s method ulam1964collection ; Dellnitz1999 in the literature. The method is not limited to the use of characteristic functions and has also been referred to as the Extended Dynamical Mode Decomposition Klus2015a ; Williams2015 . In particular, Ulam’s method has been applied to the Lorenz flow in Lucarini2016 to calculate its linear response to forcing from transition matrices approximating the transfer operators. See also Shutte1999 , Chekroun2014 and Chekroun2016 for generalizations to high-dimensional systems and Tantet2015 for an application to atmospheric regimes detection.

Following this approach, the transfer operators $\mathcal{P}_{t}^{\eta}$, defined by (A.3) with respect to some probability measure $\eta$, are projected on a truncated family $G=\{\chi_{1},...,\chi_{n}\}$ of orthogonal basis functions. That is, $\left<\chi_{i},\chi_{j}\right>_{\eta}=0$ when $i\neq j$, where the scalar product is defined as $\left<f,g\right>_{\eta}=\int_{X}f(x)g(x)\eta(dx)$ and $\eta(f)=\int_{X}f(x)\eta(dx)$. Any vector of components $\mathbf{f}=(f_{1},...,f_{n})$ in $\mathbb{R}^{n}$, defines an observable

[TABLE]

For any such $f$, the component of $\mathcal{P}_{t}^{\eta}f$ on the basis function $\chi_{j}$ is then given by

[TABLE]

where we have defined the transition matrix $\mathbf{P}^{\eta}_{t}$ with elements the normalized correlations

[TABLE]

In this study, we will only consider families of characteristic functions $G=\{\mathbf{1}_{B_{1}},...,\mathbf{1}_{B_{n}}\}$ on a grid of disjoint boxes $\{B_{1},...,B_{n}\}$ such that $\cup_{i=1}^{n}B_{i}\subseteq\mathcal{B}(R_{o})$, where $\mathcal{B}(R_{o})$ is the Borel $\sigma-$algebra of the ball $R_{o}$ defined in section 2 444In the uniform discretization of the ball $B_{o}$ in spherical coordinates, in section 3.3, different grid boxes will have different Lebesgue volumes. The component $f_{i}$ of an observable $f$ being the integral of $f$ over the box $B_{i}$, this component will tend to be larger for boxes of larger volume. This does not affect the eigenvalues, but clearly affects the eigenvectors, where large values of the components might be associated to large boxes, ceteris paribus.. In this case, the correlations are simply given by the transition probabilities

[TABLE]

Discrete approximations of the eigenvalues, eigenvectors and adjoint eigenvectors of $\mathcal{P}_{\tau}$ can then be obtained from those of a transition matrix $\mathbf{P}^{\eta}_{\tau}$ at some time $t=\tau$. Assuming that the family of transition matrices $\mathbf{P}^{\eta}_{t},t\geq 0,$ preserves the semigroup property of $\mathcal{P}^{\eta}_{t},t\geq 0,$ the spectral mapping theorem (Engel2001, , Chap. I.3, IV) ensures that the eigenfunctions of the generator of the semigroup are given by those of $\mathbf{P}^{\eta}_{\tau}$ and that the generator eigenvalues $\lambda_{k},k\geq 0,$ can be calculated from the eigenvalues $\zeta_{k}(\tau)$ of $\mathbf{P}^{\eta}_{\tau}$ according to

[TABLE]

The eigenvalues of the generator, as poles of the resolvent of the semigroup, correspond to the Ruelle-Pollicott resonances. If the semigroup property is preserved, they should not depend on the transition time $\tau$ and are thus more amenable to analysis 555The dependence of the $\lambda_{k}$’s on $\tau$ thus constitute an important test for the quality of the approximations. This is discussed in appendix B.. Before that, however, it is necessary to estimate the transition probabilities in $\mathbf{P}^{\eta}_{\tau}$.

3.1 Estimation of $\mathcal{P}^{m}_{t}$ from short time series sampling the phase space

In the original version of Ulam’s method (see e.g. Dellnitz1999 ), the transition probabilities are estimated from an ensemble of $N_{s}$ short time series $x^{(i)}_{t},0\leq t\leq\tau$, $1\leq i\leq N_{s}$, integrated numerically, with initial states sampling a given volume in phase space ($R_{o}$, for the Lorenz flow). The transition probabilities are then estimated by

[TABLE]

Because the initial states uniformly sample the Lebesgue measure $m$ (restricted to $R_{o}$), the transition matrix $\mathbf{P}^{m}_{\tau}$ gives an approximation of the transfer operator $\mathcal{P}^{m}_{\tau}$, defined in (A.3) as the dual of the Koopman operator $\mathcal{U}^{m}_{\tau}$ with respect to the scalar product induced by $m$.

3.2 Estimation of $\mathcal{P}_{t}^{\mu}$ from a long time series on the attractor

While the semigroup of transfer operator $\mathcal{P}_{t}^{m},t\geq 0$ encapsulates all the information necessary to propagate densities in phase space, we now take the point of view of an observer of a system in a statistical steady-state. The latter has then only access to time series of trajectories on the attractor and sampling the physical measure $\mu$ (the statistical steady-state) associated to it (see remarks 1 and 2). In particular, according to the formula (A.9), this is all one needs to calculate the correlation function (A.7).

It turns out that this information is encapsulated in the semigroup $\mathcal{P}_{t}^{\mu},t\geq 0$, defined in (A.3) as the dual of the Koopman semigroup $\mathcal{U}^{\mu}_{t}$ with respect to the scalar product induced by $\mu$. It follows that the transition matrix $\mathbf{P}^{\mu}_{\tau}$ approximating the transfer operator $\mathcal{P}_{\tau}^{\mu}$ can be estimated from a single long time series $x_{t_{s}},0\leq t_{s}\leq T_{\mathrm{samp}}$, of length $T_{\mathrm{samp}}$ and converged to the attractor $\Lambda$, according to

[TABLE]

The second equality in (3.4) is valid only for the case of characteristic basis functions so that the transition probabilities are given by the number of times the time series transits from one box to another in a time $\tau$ normalized by the number of trajectories in the initial box. In practice, the time series is sampled at a finite rate for a finite time, so that the time means yield the Maximum Likelihood Estimator (MLE, see e.g Billingsley1961 ) of the transition probabilities.

In this study, the invariant measure $\mu$ of interest is that supported by the chaotic set $\Lambda$, so that it may only be physical when $\Lambda$ is attracting, i.e. after the crisis, for $\rho>\rho_{A}$.

3.3 Experiment design

In order to discuss the changes in the evolution of statistics during the crisis, in the next section 4, the spectrum of both semigroups $\mathcal{P}_{t}^{m}$ and $\mathcal{P}_{t}^{\mu}$ are approximated via the methods presented in the previous sections 3.1 and 3.2. Before that, let us summarize the choice of the numerical parameters used to obtain these results. This choice is based on careful robustness tests given in appendix B and can be divided into five items:

Numerical integration: we use a Runge-Kutta sheme of order four with a time step of $10^{-4}$ time units. 2. 2.

Grid: the ball $R_{o}$ (see Sect. 2) is discretized into $3.2\cdot 10^{7}$ boxes by dividing each spherical coordinate $(r,\theta,\phi)$ in 400, 200 and 400 intervals of the same length, respectively. 3. 3.

Transition time: the generator eigenvalues were calculated from transition matrices estimated for a transition time $\tau=0.05$ time units. 4. 4.

Number and length of trajectories: $6.4\cdot 10^{9}$ trajectories of length the transition time $\tau$ and uniformly sampling $R_{o}$ were integrated in order to estimate $\mathbf{P}^{m}_{\tau}$, following section (3.1). To estimate $\mathbf{P}^{\mu}_{\tau}$, following section 3.2, 24 trajectories of length $T_{\mathrm{samp}}=1\cdot 10^{5}$ where integrated, with a spinup of $10^{4}$ time units removed. 5. 5.

Eigenproblem solver: to solve the eigenproblem for the transition matrices, the block Krylov Schur algorithm implemented in the Anasazi package Baker2009 of the Trilinos library Heroux2003 was used. The generator eigenvalues are then calculated according to (3.2) from those of the transition matrix.

4 Results

In this main section, the generator eigenvalues, calculated according to (3.2) and following both methodologies presented in section 3, are analyzed for varying control parameter values $\rho$ around the attractor crisis. The difference between the behavior of the eigenvalues approximating the stable and the unstable resonances will be emphasized and the consequences regarding critical slowing down discussed.

4.1 Stable and unstable resonances from $\mathbf{P}^{m}_{\tau}$

The approximations of both the stable and unstable resonances from the eigenvalues of the transition matrix $\mathbf{P}^{m}_{\tau}$ (Sect. 3.1), are represented in figure 6, for values of $\rho$ ranging from 22 (between the homoclinic bifurcation and the boundary crisis) and 25 (after the boundary crisis). Overall, a nontrivial arrangement of eigenvalues is found, as expected from the presence of the chaotic set $\Lambda$.

In order to focus on the generator eigenvalues closest to the imaginary axis, expected to be affected by the crisis, the real parts of the leading eigenvalues are represented in figure 7 versus $\rho$. One can see that for $\rho$ smaller than $\rho_{A}\sim 24.06$, there are two eigenvalues zero. As $\rho$ is increased to $\rho_{A}$, a third eigenvalue approaches the imaginary axis. As $\rho$ is increased further that $\rho_{A}$, there is an exchange of eigenvalues, so that only one eigenvalue remains zero and two eigenvalues get further from the imaginary axis. There are also other eigenvalues approaching and then escaping from the imaginary axis, although with a gap remaining finite.

Zero eigenvalues are associated with physical measures corresponding to attractors. For $\rho<\rho_{A}$ there are two attractors, namely the two stationary points, so that it does not come as a surprise that two zero eigenvalues are found. For $\rho_{A}<\rho<\rho_{\mathrm{Hopf}}$, however, there are three coexisting attractors, namely the two stationary points and the strange attractor $\Lambda$, so that one would expect three zero eigenvalues instead of only one. It thus appears that the discretization does not allow to resolve the boundary of the different basin of attractions. One can interpret the errors induced by the discretization as numerical diffusion Froyland2011a hindering the distinction of the three different basins of attraction. The three attractors are, at our level of description, merged as a result of diffusion.

A more detailed understanding is obtained from the eigenvectors, represented in figure 8, associated with the three leading generator eigenvalues (from top to bottom), calculated from the transition matrices $\mathbf{P}^{m}_{\tau}$ for $\rho=22$ (left) and $\rho=24.5$ (right). The first eigenvectors (top) give an approximation of the physical measure of the system by a Lebesgue-density which is stationary, as far as the numerical approximations are concerned. One can see that, for $\rho<\rho_{A}$ (left), this stationary density has two modes over the two stable stationary points $p^{\pm}$. There is, however, some spread around these points, due to the numerical approximation. The second eigenvector is also nonzero only about $p^{\pm}$, but with opposite values over each stationary point 666Note that the secondary eigenvectors being orthogonal to the leading adjoint eigenvectors, which is constant (see Appendix A.4), they integrate to zero with respect to the Lebesgue measure.. It thus accounts for the fact that transitions between these points are very rare (and only possible due to the discretization). The third eigenvector, however, is positive about the stationary points but negative about the chaotic set $\Lambda$. For $\rho$ close, but smaller than $\rho_{A}$, the set $\Lambda$ is invariant but weakly unstable, so that the measure $\mu$ supported by $\Lambda$ is not yet physical. Transitions to $p^{\pm}$ from points around $\Lambda$, but not exactly on it, are thus possible although very slow, since the crisis is near. This explains the dichotomy between $p^{\pm}$ and $\Lambda$ visible in the third eigenvector.

For $\rho_{A}<\rho<\rho_{\mathrm{Hopf}}$ (right panels) the second and third eigenvectors are similar to those for $\rho<\rho_{A}$. However, the stationary density (top right panel) now spreads about $\Lambda$, in agreement with the fact that $\Lambda$ is now an attractor. For $\rho<\rho_{A}$, one can thus interpret the real part of the third eigenvalue as the escape rate from the chaotic set $\Lambda$ to each stationary point, while for $\rho>\rho_{\mathrm{Hopf}}$, the second and third eigenvalues can be interpreted as the escape rate from the stationary points to $\Lambda$. However, while this interpretation is valid for $\rho>\rho_{\mathrm{Hopf}}$, for the value of $\rho$ between $\rho_{A}$ and $\rho_{\mathrm{Hopf}}$ represented in the left panels, the stationary points $p^{\pm}$ are still attracting. One would thus expect three eigenvalues to be zero (see appendix A.4) instead of one (Fig. 7). It thus appears that the discretization is still too coarse to get around numerical diffusion Froyland2011a preventing the resolution of the three distinct basins of attraction. This difficulty is likely to arise from the small size of the basins of attraction of $p^{\pm}$ close to the Hopf bifurcation and from the possibly convoluted geometry of their boundaries.

As a conclusion, the first key result of this study is the clear indication of the crisis in the shrinkage of the spectral gap in the approximation of the generator eigenvalues for $\mathcal{P}_{t}^{m},t\geq 0$. However, as discussed in section 3.2 (see also appendix A.3), for this change in the spectrum to be detectable from observations on the attractor alone, the eigenvalues of the generator of $\mathcal{P}_{t}^{\mu},t\geq 0$ with respect to the physical measure $\mu$, for $\rho>\rho_{A}$, should also be affected. This is tested in the next section 4.2.

4.2 Unstable resonances from $\mathbf{P}^{\mu}_{\tau}$ and decay of correlations

To better understand which of the changes in the evolution of statistics are available from observations on the attractor, the unstable resonances alone are calculated from the generator eigenvalues of $\mathbf{P}^{\mu}_{\tau}$ (Sect. 3.2). The resulting generator eigenvalues are represented in the left panels of figure 9, for values of $\rho$ ranging from about $\rho_{A}$ to $\rho=28$. As recalled in appendix A.3, the spectrum of the transfer semigroup $\mathcal{P}_{t}^{\mu},t\geq 0,$ governs the decay of correlations between any pair of observables in appropriate functional spaces. As an example, the right panels of Fig. 9 represent sample estimates 777The sample correlation function is estimated from a long time series initialized in the basin of attraction of $\Lambda$ as a discrete approximation of the time average (A.9) (see e.g. VonStorch1999b ). of the correlation function $C_{z,z}(t)$ for the observable $z:(x,y,z)\to z$. This is but one example of correlation function, which will be sufficient for the present discussion. To facilitate the comparison with the results of the previous section 4.1, the corresponding generator eigenvalues from $\mathbf{P}^{m}_{\tau}$ have also been represented by black crosses in the top left panel in figure 9, for $\rho=24.1$.

As discussed in appendix A.4, we expect the eigenvalues of $\mathbf{P}_{\tau}^{\mu}$ to constitute a subset of the eigenvalues of $\mathbf{P}_{\tau}^{m}$, within numerical errors. This seems to be indeed the case as the generator eigenvalues calculated from $\mathbf{P}^{\mu}_{\tau}$ (blue dots in the top left panel) roughly correspond to some of the generator eigenvalues calculated from $\mathbf{P}^{m}_{\tau}$ (black crosses). However, the complex plane is less densely populated by generator eigenvalues of $\mathbf{P}^{\mu}_{\tau}$ than of $\mathbf{P}^{m}_{\tau}$. This is not surprising, since the transfer operators $\mathcal{P}_{t}^{\mu}$ only gives access to the unstable resonances and not to the stable resonances describing the relaxation of densities transversally to $\Lambda$.

The second key result of this study is that the generator eigenvalues of $\mathbf{P}^{m}_{\tau}$, which where found in section 4.1 to approach/escape from the imaginary axis during the crisis, do not belong to the set of generator eigenvalues of $\mathbf{P}^{\mu}_{\tau}$. In fact, the spectral gap between the leading generator eigenvalues of $\mathbf{P}^{\mu}_{\tau}$, in figure 9, is not significantly affected by the attractor crisis. As a result, neither is the rate of decay of the correlation function $C_{z,z}(t)$, whose sample estimates are represented in the right panels of the same figure 9. This applies to the correlation function between any pair of observables (not shown here). In conclusion, for this attractor crisis in the Lorenz flow, critical slowing down is not observable from time series converged to $\Lambda$. This can be interpreted by the fact that the evolution of densities is affected by the weaker stability of the attractor close to the crisis, but that this change in stability is not felt along the attractor. The stability of the attractor should here be understood as the rate of convergence of densities to it, as measured by the spectral gap between the zero and the secondary generator eigenvalues of $\mathcal{P}_{t}^{m}$, as opposed to the rate of divergence of nearby trajectories measured by the Lyapunov exponents. In conclusion, only the stable resonances are affected by the crisis.

5 Summary and discussion

‘ The global stability of the chaotic attractor in the Lorenz flow as it undergoes a boundary crisis is investigated numerically in terms of Ruelle-Pollicott resonances. The results are now summarized and their implications regarding early-warning of chaotic attractor crises and response theory discussed.

For the boundary crisis of the Lorenz flow, neither the Lyapunov exponents nor the covariant Lyapunov vectors provide a precursor. For that purpose, global information appears to be necessary. While the Lyapunov spectrum describes the divergence of nearby trajectories, the semigroup of transfer operators governs the convergence/escape of densities to/from an invariant measure. The rate at which this occurs is characterized by the eigenvalues of the generator of the semigroup with negative real part, the Ruelle-Pollicott resonances. These resonances are divided into a stable and an unstable family describing the evolution of densities about and on the support of the invariant measure, respectively. The gaps between the stable (unstable) resonances and the imaginary axis thus provides an indicator of the global stability (mixing rate) of a possibly chaotic invariant set.

To monitor the changes in the spectrum during the crisis, a discretization of the transfer operators is estimated either from many short time series sampling the phase space or from a few long trajectories converged to the attractor. The difference between the two approaches is essential, as the first one yields an approximation of transfer operators acting on functions of the full phase space, thus yielding information on global stability, while the second one gives a restriction of the transfer operators to the attractor supporting a physical measure. The second approach thus only allows to approximate the unstable resonances associated with mixing dynamics on the attractor, while the first also gives an approximation of the stable resonances corresponding to contraction towards the attractor, thus providing global information on the properties of the system outside the attractor.

The main result of this study is twofold. First, as the crisis is approached, some stable resonances approach the imaginary axis, as indicative of the weaker stability of the attractor. After the crisis, these stable resonances get further from the imaginary axis, their distance to the imaginary axis giving a measure of the escape rate from the chaotic saddle.

Second, only the stable resonances are affected by the crisis. This implies that, in the boundary crisis considered here, the approach of the crisis cannot be inferred from time series evolving along the attractor, i.e. from long observations of the unperturbed system. This explains why the correlation functions between different observables are not showing any sign of slowing down of their decay before the crisis.

As a consequence, early-warning indicators Scheffer2009 based on correlation functions Held2004a or power spectra Kleinen2003a are in general unable to give a precursor of attractor crises showing the same behavior as the one studied here. On the other hand, indicators based on the recovery from perturbations Nes2007 may be useful. The perturbation should then be transverse to the attractor, in the sense that it should project on the stable manifold of the attractor. In fact, in the case of a stationary point, the recovery rates as given by the eigenvalues of the Jacobian of this point Nes2007 coincide with the leading Ruelle-Pollicott resonances Gaspard1995 ; Lan2013 . A conclusion of this study is that, in the chaotic case, the recovery rates are not given by the Lyapunov exponents, but rather by the stable resonances.

In fact, that some resonances correspond to motions on the attractor, which can be observed from long time series, while other correspond to motions transverse to the attractor, is the central difference between the response theory for systems in thermodynamic equilibrium Kubo1957a and Ruelle’s response theory for dissipative systems Ruelle2009 . Indeed, invariant measures of equilibrium systems have a density with respect to the Lebesgue measure, allowing for the expression of the linear term in the response to forcing in terms of correlation functions. This is the statement of the celebrated flucutation-dissipation theorem. On the other hand, for dissipative systems such as the Lorenz flow Lucarini2009b , the invariant measure supported by the attractor is singular and an additional term accounting for perturbations along the stable manifold of the attractor must be taken into account. In any case, the susceptibility of the system to forcing depends on the presence of resonances close to the imaginary axis Ruelle2009 ; Cessac2007 . Eventually, when resonances touch the imaginary axis, response theory breaks down, which is what is observed in this study during the attractor crisis, where response theory cannot be expected to work.

The story can, however, be very different in the presence of noise. Noisy perturbations may indeed push the system away from the attractor, in the direction of its stable manifold, allowing for the system to explore the phase space. This argument can be made rigorous by verifying the Lie bracket Hörmander condition Hormander1968a ; Hairer2011 to ensure, from the hypoellipticity of the generator, that, thanks to the noise, the transfer semigroup has a smoothing effect. If this condition is verified, the invariant measure may admit a density, allowing for the fluctuation-dissipation theorem to hold Hairer2010 . In this case, the transfer semigroup in phase space and the one restricted to the support of the invariant measure $\mu$ can be identified with the same semigroup on $L^{2}_{\mu}$. As a consequence, recovery rates and decorrelation rates can both be identified with the real part of the leading generator eigenvalues Chekroun2016 . Moreover, when the unperturbed deterministic system is hyperbolic, singular perturbation theory may be applied to relate the eigenvalues of the perturbed system to those of the unperturbed system Gaspard2002a ; Tantet2016 .

A final issue remains regarding the genericity of the results obtained here for the Lorenz flow. The crisis scenario described here might be of general relevance for the study of three-dimensional flows, since the robust chaotic sets in three dimensions are singular hyperbolic. It is not a priori clear whether the basic difference between the behavior of stable vs. unstable resonances near the crisis persists in higher dimension. In particular, indications were given in Tantet2015a of the shrinkage of the spectral gap of transfer operators and of the slowing down of the decay of correlations during a boundary crisis in a deterministic and autonomous general circulation climate model, featuring exclusively internally generated fluctuations. Nonetheless, those results are not conclusive enough, because we looked at the transfer operator in a severely projected space and we could not follow the statistical properties mentioned above up to the exact value of the parameter determining the crisis.

A next step in our analysis is to try to understand what determines whether critical transitions due to boundary crises in a deterministic chaotic system are flagged by the properties of the unstable resonances or whether looking at the stable resonances is instead required. These two scenarios pertain to cases where the attractor contains all the needed information, or, instead, a neighborhood of the attractor needs to be considered. In the latter case, a natural way to probe the system is by adding a small stochastic perturbation.

A possible way ahead is to look at whether a role in this sense is played by different scenarios pertaining to the geometrical properties of the basin boundary, and in particular at looking at the details of how the edge state, i.e. the unstable periodic orbits in the case studied here, and the attractor collide at the bifurcation. It is important to note that the edge state separating two co-existing attractors is the gate for the noise-induced transitions between the attractors. Therefore, near the crisis, adding noise does not only allow for probing a small region near the attractor, but also leads to sampling the properties of edge state.

Appendix A Transfer operators, resonances and decay of correlations

Let us shortly review the spectral theory of transfer operators relevant for this study. A more detailed exposition can be found in Tantet2015a and references therein.

A.1 Semigroup of transfer operators and ensemble

We consider a dynamical system on a phase space $X$ in $\mathbb{R}^{d}$ and with global invertible flow $\Phi_{t}:x_{0}\to x(t),t\geq 0$ on $X$ generated by a sufficiently smooth autonomous vector field $F:X\to X$, i.e.

[TABLE]

In other words, the transformation $\Phi_{t}$ associates to any initial condition $x_{0}$ in $X$ the corresponding solution to (A.1) at time $t$.

When the system is chaotic, it is fruitful to follow the evolution of probability densities and observables rather than that of solutions of (A.1). For that purpose, we endow the phase space $X$ with its Borel $\sigma$-algebra $\mathcal{B}$ and some probability measure $\eta$ (to be specified) on $(X,\mathcal{B})$. The flow $\Phi_{t},t\geq 0$ 888 The flow will be assumed to be nonsingular for the Lebesgue measure, so that it maps sets of null measure into sets of null measure. induces a family of linear operators, the Koopman operators,

[TABLE]

acting on observables $g$ in the space bounded measurable functions $L^{q}_{\eta}(X),1\leq q<\infty$. On the other hand, there exists a family of linear operators $\mathcal{P}^{\eta}_{t},t\geq 0$, the transfer operators or Perron-Frobenius operators on $L^{p}_{\eta}(X)$, with $1/p+1/q=1$, and such that the duality relation

[TABLE]

holds for any observable $f$ in $L^{p}_{\eta}(X)$ and $g$ in $L^{q}_{\eta}(X)$. It follows that $\mathcal{U}^{\eta}_{t}$ is the adjoint of $\mathcal{P}^{\eta}_{t}$. Moreover, taking $g$ in (A.3) as the indicator $\mathbf{1}_{A}$ of some set $A\in\mathcal{B}$, one has that 999For a probability density $f$, i.e. $f\geq 0$ and $\int_{X}f(x)\eta(dx)=1$, (A.4) expresses the fact that the probability to find a member $x$ sampled from an initial ensemble $f$ in a set $A$ after some time $t$ is none-other than the probability of this member to be initially in the preimage of this set by the flow.

[TABLE]

The families of linear operators $\mathcal{P}^{\eta}_{t}$ and $\mathcal{U}^{\eta}_{t},t\geq 0$ inherit from the semigroup property of the flow, thus imposing a strong constrain on the operators for different times.

The relationship between the nonlinear flow $\Phi_{t}$ and these semigroups provides the connection between the ergodic properties of dynamical systems and the functional analysis of linear operators Yosida1980 ; Engel2001 ; Davies2007 , as was first recognized in Koopman1931 and Neumann1932 .

A.2 Ergodicity, mixing and correlations

A key concept in ergodic theory is that of an invariant measure for the flow $\Phi_{t}$, that is, a probability measure $\mu$ such that

[TABLE]

In other words, a measure is invariant if, according to this measure, the probability of a state to be in some set does not change as this state is propagated by the flow. It follows then that the average with respect to the invariant measure $\mu$ of any integrable observable $g$ is also invariant with time, i.e.

[TABLE]

so that the invariant measure gives a statistical steady-state.

A flow $\Phi_{t}$ with an invariant measure $\mu$ has interesting statistical properties when $\mu$ is ergodic, that is, when the sets $A$ which are invariant, i.e. $\Phi_{t}^{-1}A=A$, are either of measure 0 or 1 101010In particular, the ergodicity of the invariant measure implies that each set of positive measure is visited infinitely often by orbits starting from $\mu$-almost every point (Eisner2015, , Lemma 6.15).. Then, by the celebrated individual ergodic theorem of Birkhoff (Lasota1994, , Chap. 7.3), the average of any $\mu$-integrable observable $g$ is such that

[TABLE]

Thus, when $\mu$ is ergodic, the time mean is independent of the initial state $x$ except for a set of null measure.

Remark 1

There may exist many ergodic measures. This is actually the case for the Lorenz flow for $\rho_{\mathrm{Homo}}<\rho<\rho_{A}$, where, as seen in section 2, three attractors coexist. Each attractor then supports at least one invariant measure (Katok1996, , Chap. 4). In this case, the equality (A.6) between ensemble averages and time averages may hold only for initial states belonging to the attractor supporting the measure, while time averages for two initial conditions in different basins of attraction will not coincide in general. More useful for experiments is then the eventual physical property of the measure. The latter ensures that the equality (A.6) between ensemble averages and time averages holds not only for initial states in a set of positive measure $\mu$, but also for states in any set of positive Lebesgue measure $m$ in the basin of attraction of a given attractor Young2002 .

A particularly important quantity in ergodic theory is the correlation function,

[TABLE]

between any observables $f,g\in L^{2}_{\mu}(X)$. It gives a measure of the relationship between the two observables as one evolves with time. The decay of correlation functions with time is a macroscopic manifestation of chaos. Indeed, this decay is equivalent to the mixing property (see e.g. (Lasota1994, , Chap. 4) and (Katok1996, , Chap. 4)),

[TABLE]

of the invariant measure $\mu$, which is stronger than ergodicity. In other words, the probability for some state to be in any set $B$ after some time $t$ is independent of the probability of this state to be in any set $A$ initially. Any information about the initial state of an ensemble is thus gradually forgotten.

Remark 2

For a physical measure $\mu$ supported by an attractor $\Lambda$, the correlation function can be estimated by the time mean

[TABLE]

from a single time initialized in the basin of attraction of $\Lambda$. For this reason, we will see in section 3.2 that complete information on $\mathcal{P}^{\mu}_{t},t\geq 0$ can be obtained from from observations on the corresponding attractor alone. Unfortunately, for the forced and dissipative systems considered here, volumes contract on average (Gallavotti2014a, , Chap. 2.8) so that the invariant measure $\mu$ is supported by an attractor with zero Lebesgue measure $m$. It follows that, as opposed to the semigroup $\mathcal{P}^{m}_{t},t\geq 0$, no information on the dynamics away from the attractor is carried by $\mathcal{P}^{\mu}_{t}$.

A.3 Spectral theory of ergodicity and mixing

One can see from the definition (A.7) that the correlation function is fully determined by the family of transfer (Koopman) operators $\mathcal{P}_{t}^{\mu}$ ($\mathcal{U}_{t}^{\mu}$) on $L^{2}_{\mu}(X)$, i.e. with respect to $\mu$. In fact, classical results relate the ergodic properties of measure-preserving dynamical systems to the behaviour of the semigroups Halmos1956 ; Arnold1968 ; Lasota1994 . Note first, that the invariance of the measure $\mu$ together with the invertibility of the flow ensure that the semigroups are isometries and constitute a unitary group on $L^{2}_{\mu}(X)$. As a consequence, the spectrum of the operators is contained in the unit circle $|z|=1,z\in\mathbb{C}$ and $\mathcal{P}_{t}^{\mu}$ and $\mathcal{U}_{t}^{\mu}$ have the same eigenvalues and eigenfunctions, for any real $t$. Moreover, the following results hold.

Theorem A.1

If the measure-preserving dynamical system $(X,\mathcal{B},(\Phi_{t})_{t\in\mathbb{R}},\mu)$ is

1. (i)

ergodic, then $1$ is a simple eigenvalue of $\mathcal{P}_{t}^{\mu},t\in\mathbb{R}$, and conversely. 2. (ii)

mixing, then the only one eigenvalue is one.

The first item of this theorem can be understood from the fact that ergodic systems have only one stationary density with respect to the invariant measure. The second is due to the fact that eigenfunctions associated with eigenvalues on the unit circle do not decay and thus prevent mixing for general observables.

The rate of decay of correlations, or mixing, characterizes the (weak) convergence of ensembles to a statistical steady-state and has been the subject of intense research these past decades. Once again, the latter can be studied from the spectral properties of the semigroups Liverani1995a . However, the eigenvalues responsible for such decay, the Ruelle-Pollicott resonances Pollicott1985 ; Ruelle1986 , lie inside the unit disk and are thus not accessible from operators on regular functional spaces. Anisotropic Banach spaces of distributions capturing the dynamics of contraction and expansion of the system and on which the semigroups are contracting should instead be considered Liverani1995a ; Blank2001a ; Gouezel2006 ; Butterley2007 ; Faure2014 ; Baladi2017 111111Contraction by semigroups can be interpreted in physical systems as determining entropy production. Recognizing this has been essential to understand how reversible microscopic evolutions can lead to irreversible macroscopic properties of systems out of thermal equilibrium (see the pioneering work Misra1979 and Gaspard1998 ; Garbaczewski2002 for reviews)..

A certain degree of hyperbolicity is necessary to adapt these spaces to the dynamics, and the robustness to perturbations of the spectrum or the convergence of numerical algorithms is not guarantied in general (see Keller1998 ; Gouezel2006 and Baladi1999 ; Froyland2007a , respectively, for results in this direction). In the present study, it is found, however, that the discretization scheme used in section 3 yields reasonably good approximations of the resonances.

The spectral properties of the semigroups with respect to an invariant measure $\mu$ thus allows for the study of the ergodic properties of this measure. On the other hand, for the study of the global stability of this set, operators acting on a neighborhood of positive Lebesgue measure of a chaotic set should be considered.

A.4 Spectral theory of global stability

In the introduction, we have advocated on heuristic grounds that the global stability of invariant sets could be studied from the evolution of densities and observables. In this case, however, the transfer (Koopman) operators $\mathcal{P}_{t}^{m}$ ($\mathcal{U}_{t}^{m}$) with respect to the Lebesgue measure $m$ should be considered. Indeed, this should allow one not only to study the mixing dynamics of an invariant measure, but also the contraction to/escape from such a measure. This has allows for new developments in the theory of the global stability of stationary points and periodic orbits from the behavior of $\mathcal{P}_{t}^{m}$ Vaidya2008 and $\mathcal{U}_{t}^{m}$ Mauroy2016 . This approach should also be amenable to chaotic invariant sets. However, the same difficulties as mentioned in the previous section A.3, regarding the functional analytic framework appropriate to capture the resonances inside the unit disk, are encountered. In order to better interpret the numerical results obtained in section 4, let us however give a few comments regarding the eigenfunctions of the Koopman operators $\mathcal{U}^{X}_{t}$ and $\mathcal{U}^{\Lambda}_{t}$ on the space $C(X)$ of continuous functions on the compact metric space $(X,d)$ (e.g. $X=\mathbb{R}^{n}$ with $d$ the distance induced by the Euclidean norm) and its restriction to the support $\Lambda$ of some invariant measure $\mu$, respectively. Since $(X,d)$ is a normal space and $\Lambda$ is a closed subset of $X$, the Tietze-Urysohn Extension Theorem (e.g. Lang1993 ) ensures that for any continuous function $g^{\Lambda}$ on $\Lambda$, there exists a continuous function $g^{X}$ on $X$ whose restriction to $\Lambda$ is $g^{\Lambda}$. The restriction $\mathcal{U}^{\Lambda}_{t}$ of $\mathcal{U}^{X}_{t}$ on $C(\Lambda)$ is thus defined such that $\mathcal{U}^{\Lambda}_{t}g^{\Lambda}=\mathcal{U}^{X}_{t}g^{X}$.

A.4.1 Multiple attractors

We first discuss the coexistence of multiple attractors with disjoint basins of attraction, (such as in the Lorenz flow for $\rho_{\mathrm{Homo}}<\rho<\rho_{A}$), in terms of multiplicity of the eigenvalue $1$. When the phase space $X$ is forward invariant, i.e. such that $X\subset\Phi_{t}^{-1}X$ then one has that

[TABLE]

Thus, the forward invariance of $X$ implies that the function $\mathbf{1}_{X}$, constant on $X$, is an eigenfunction of the Koopman semigroup associated with the eigenvalue $1$. This is for example the case when the flow has a single globally asymptotically stable attractor on $X$. If instead, there is a single repeller and $\Phi_{t}^{-1}X\subset X$, then $\mathcal{U}^{X}_{t}\mathbf{1}_{X}=\mathbf{1}_{\phi_{t}^{-1}X\cap X}=\mathbf{1}_{\phi_{t}^{-1}X}$ so that $\mathbf{1}_{X}$ is no longer an eigenfunction associated with $1$.

When several attractors $\Lambda_{1},...,\Lambda_{l}$ coexist, each basin of attraction $B(\Lambda_{i})$ is forward invariant and $\mathcal{U}^{X}_{t}\mathbf{1}_{B(\Lambda_{i})}=\mathbf{1}_{B(\Lambda_{i})},t\geq 0$. Thus, to each attractor $\Lambda_{i}$ corresponds a Koopman eigenfunction $\mathbf{1}_{B(\Lambda_{i})}$ associated with a generator eigenvalue zero. However, only $l$ of the $l+1$ functions $\mathbf{1}_{X},\mathbf{1}_{B(\Lambda_{1})},...,\mathbf{1}_{B(\Lambda_{l})}$ are independent, so that there are only $l$ eigenvalues 1. For example, when two attractors coexist, a possibility is to have as independent eigenfunctions $\mathbf{1}_{X}$ and $\mathbf{1}_{B(\Lambda_{1})}-\mathbf{1}_{B(\Lambda_{2})}$. Adding suitably defined (e.g. Gaussian) noise leads to a unique invariant measure; reducing the intensity of the noise to zero leads in the limit to selecting a special invariant measure constructed as linear combination of all the invariant measure of the deterministic case. Note, however, that the numerical approximations of the basins of attraction may be difficult when the geometry of the boundary is convoluted. The fact that the numerical discretization of the phase space leads effectively to introducing some noise in the system might explain why also in the deterministic case one can miss the presence of various coexisting and independent invariant measures.

A.4.2 Correspondence between the eigenvalues of $\mathcal{U}_{t}^{\Lambda}$ and $\mathcal{U}_{t}^{X}$

Let us discuss the correspondence between the eigenvalues of $\mathcal{U}_{t}^{X}$ and $\mathcal{U}_{t}^{\Lambda}$. Following Cessac2007 , we refer to the eigenvalues of $\mathcal{U}_{t}^{X}$ as the unstable resonances and those of the $\mathcal{U}_{t}^{\Lambda}$ which do not correspond to the latter as the stable resonances.

What we show below is strictly applicable only for eigenvalues on the unit circle, as discussed in Remark 3. Nonetheless, we believe that it may be practically relevant in the region near the unit circle, and, using more advanced mathematical tools, could be extended for the unit disk, away form the essential spectrum.

We consider the particular case where the invariant measure $\mu$ is supported by a uniformly hyperbolic (compact invariant) set $\Lambda\subset X$ for the continuously differentiable flow $\Phi_{t},t\in\mathbb{R}$ on $(X,d)$. For convenience, we will assume that $X$ is the basin of attraction of $\Lambda$. The global stable manifold $W_{s}$ through the point $x\in X$ can then be characterized topologically by

[TABLE]

It follows directly that $W_{s}$ is invariant, i.e. $\Phi_{t}W_{s}(x)\subset W_{s}(x)$ and $\cup_{x\in\Lambda}W_{s}(x)=X$, since $X$ is the basin of attraction of $\Lambda$. Moreover, as a consequence of the stable manifold theorem (Katok1996, , Chap. 6), $x\mapsto W_{s}(x)$ is continuous.

Proposition A.2

Let $\Phi_{t},t\in\mathbb{R}$ be a continuously differentiable flow on the compact space $X$. Let $\Lambda$ be a uniformly hyperbolic compact invariant set and denote by $\mathcal{U}_{t}^{\Lambda}:C(\Lambda)\to C(\Lambda),t\geq 0,$ the restriction to $C(\Lambda)$ of the Koopman operator $\mathcal{U}^{X}_{t}:C(X)\to C(X),g\mapsto g\circ\Phi_{t}$. Assume that $\psi^{\Lambda}$ in $C(\Lambda)$ is an eigenfunction of $\mathcal{U}_{t}^{\Lambda}$ associated with the eigenvalue $\zeta\in\mathbb{C}$ for some $t>0$, i.e. $\mathcal{U}_{t}^{\Lambda}\psi^{\Lambda}=\zeta\psi^{\Lambda}$. Then the function $\psi^{X}$ such that

[TABLE]

is in $C(X)$ and is an eigenfunction for $\mathcal{U}_{t}^{X}$ associated with the eigenvalue $\zeta$.

In other words, $\psi^{X}$ takes on a leaf of the global stable manifold a constant value given by that of $\psi^{\Lambda}$ on this leaf.

Proof

Let us first verify that $\psi^{X}$ is indeed in $C(X)$. Note first that $\psi^{X}$ is defined on $X$, since $W_{s}=X$. The continuity of $\psi^{X}$ follows from that of $\psi^{\Lambda}$ and $x\mapsto W_{s}(x)$. To see this, let $\{y_{n}\}$ be a sequence in $X$ converging to $y$. For each $y_{n}$ there is an $x_{n}$ in $\Lambda$ such that $y_{n}\in W_{s}(x_{n})$. Thus

[TABLE]

From the continuity of $x\mapsto W_{s}(x)$, the limit $x$ of $x_{n}$ exists and is such that $y\in W_{s}(x)$ and $\psi^{X}(y)=\psi^{\Lambda}(x)$. From the continuity of $\psi^{\Lambda}$, it follows that

[TABLE]

so that $\psi^{X}$ is continuous on the metric space $(X,d)$.

That $\psi^{X}$ is an eigenfunction follows directly from the invariance of the global stable manifold:

[TABLE]

Thus, to each eigenfunction $\mathcal{U}_{t}^{\Lambda}$ on $C(\Lambda)$ corresponds an eigenfunction of $\mathcal{U}_{t}^{X}$ on $C(X)$ associated with the same eigenvalue, i.e. the spectrum of $\mathcal{U}_{t}^{\Lambda}$ is a subset of the spectrum of $\mathcal{U}_{t}^{X}$.

Remark 3

There is, however, a major caveat to the applicability of this results. Indeed, it cannot be applied to eigenfunctions associated with eigenvalues inside the unit disk, since $\mathcal{U}_{t}^{\Lambda}$ has a roughening effect due to the contraction on the unstable manifold manifold of $\Lambda$, backward in time. For that purpose, spaces a distributions should be considered for $\psi^{\Lambda}$.

Remark 4

A similar result holds for the nonuniformly hyperbolic case Barreira2002 . However, the stable foliation is then only measurable so that eigenfunctions in spaces of measurable functions should be considered.

Appendix B Robustness of the resonances

In this appendix, we explain and give support to the choice of parameters used to obtain the results of section 4 and given in section 3.3.

B.1 Numerical integration

In this study, the time step of integration is particularly important for the position of the attractor crisis. It appears that for larger time steps, the attractor crisis occurs for values of $\rho$ smaller than $\rho_{A}$. A time step of $10^{-4}$ allows for the crisis to occur within $1\%$ of $\rho_{A}$, while keeping the numerical integration tractable.

B.2 Discretization

The choice of the grid used to approximate the transfer operators by transition matrices is the most critical step of the numerical application in this study. First, in order for probabilities to be conserved, the grid should cover a bounded set within which any trajectory remains. This condition if fulfilled by a grid covering the ball $R_{o}$ (see Sect. 2). A discretization of $R_{o}$ is then easily implemented when working in spherical coordinates $(r,\theta,\phi)$.

Second, due to the fine-grained geometry of the eigenvectors associated with chaotic sets, the generator eigenvalues are slow to converge with the grid resolution 121212There is no general result regarding the convergence with the grid resolution of eigenvalues of transition matrices to the Ruelle-Pollicott resonances. This is, however, the case for uniformly hyperbolic systems Gouezel2006 ; Butterley2007 for which Ulam’s method may converge (see Froyland2007a , for the case of hyperbolic maps).. In this study, the focus is, however, on the generator eigenvalues close to the imaginary axis, which are expected to be more robust to perturbations that eigenvalues further from the imaginary axis. The results of a test of convergence with respect to the grid resolution $n_{d}$-by-$n_{d}/2$-by-$n_{d}$ is given in figure 10. One can see that for $n_{d}\geq 400$, the real part of the first nonzero eigenvalue is close to convergence. To get an idea of the value to which this real part would converge for higher resolutions, the dashed orange line represents a least square fit of an exponential $ae^{bn_{d}},$ with $a,b\in\mathbb{R}$, to it. The quality of this fit and the fact that the fitted curve is an exponential converging to zero suggests that this eigenvalue, for $\rho\approx 24$, would converge to the imaginary axis if the resolution were to be further increased. Eventually, a grid resolution of 400-by-200-by-400 is chosen, allowing for the real part of the first nonzero generator eigenvalue to remain within $2\%$ of the corresponding value obtained for $n_{d}=500$. The same grid is used to estimate $\mathbf{P}^{\mu}_{\tau}$, for which similar numerical convergence is also observed (not shown here).

B.3 Transition time

The transition time $\tau$ for which the transition matrices $P_{\tau}^{m}$ and $P_{\tau}^{\mu}$ are estimated is a key parameter. In theory, the spectral mapping formula (3.2) allows to calculate the generator eigenvalues $\lambda_{k}$ from the eigenvalues $\zeta_{k}(\tau),k\geq 0$ for the transfer operator $\mathcal{P}_{\tau}^{\eta}$, for any transition time $\tau$. A first issue is, however, that, in taking the complex logarithm divided by $\tau$ to get the $\lambda_{k}$ from the $\zeta_{k}(\tau)$, the imaginary part of the $\lambda_{k}$ is only known modulo $2\pi/\tau$. An arbitrary choice of the principal part would then be valid only for true $\lambda_{k}$ such that $|\Im(\lambda_{k})|\leq\pi/\tau$, so that this window shortens as $\tau$ is increased (see (Crommelin2009a, , Sect. 2.4)).

Second, a compromise should be found in order to estimate correctly as many eigenvalues as possible, while only approximating those for which the eigenvectors can be resolved for a given grid. Indeed, the longer the transition time $\tau$, the smaller $|\zeta_{k}(\tau)|=e^{\Re(\lambda_{k})\tau}$ for $\lambda_{k}$ with a small real part. Thus, in order to be able to estimate the eigenvalue $\zeta_{k}(\tau)$ numerically, $\tau$ should be sufficiently small for $e^{\Re(\lambda_{k})\tau}$ to be larger than a threshold under which numerical errors become important. This threshold depends on several factors such as the sampling, the nonnormality of the transfer operators and roundoff errors Crommelin2009a .

On the other hand, it is not always a good strategy to take $\tau$ short to approximate as many $\lambda_{k}$ far from the imaginary axis as possible. Indeed, eigenvalues further from the imaginary tend to be associated with eigenvectors with more and more changes of sign. For the latter to be resolved, the grid resolution should be higher and higher. However, if for a fixed grid, eigenvalues $\zeta_{k}(\tau)$ associated with eigenvectors which cannot be appropriately resolved have not decayed, their imprecise approximation will also have an impact on eigenvalues closer to the imaginary axis. As a rule of thumb, the lower the grid resolution, the larger should $\tau$ be, so as to approximate only the eigenvalues for which the eigenvectors can be properly resolved at this resolution. For a grid of 400-by-200-by-400, a transition time $\tau$ of 0.05 time units is found to give of a good compromise.

B.4 Number and length of trajectories

Ulam’s method relies on the estimation of transition probabilities from time series. The quality of these estimations depend on the sampling. To test the sampling, one strategy is to estimate confidence intervals (see e.g. (Chekroun2014, , SI) and Tantet2015 ). Another approach, followed here, is to directly test the convergence of the generator eigenvalues with respect to the number of samples $N_{s}$. This convergence of the transition probabilities and the eigenvalues with $N_{s}$ is known to occur at a rate of $\mathcal{O}(N_{s}^{1/2})$ Billingsley1961 ; Crommelin2009a . When estimating $\mathbf{P}_{\tau}^{m}$ from many short time series, the number of samples is given by the number of trajectories. In this case, the robustness of the eigenvalues to $N_{s}$ is shown in the left panel in figure 11. One can see that a number $6.4\cdot 10^{9}$ of trajectories is more than enough for an estimation on the chosen grid of 400-by-200-400. When estimating $\mathbf{P}_{\tau}^{\mu}$ from a few long time series, the number of samples is given by the number of trajectories by their lengths $T_{\mathrm{samp}}$ divided by their sampling rate. Here, a sampling rate of $100$ samples per time unit is used and 24 long trajectories are used (in order to distribute each on a computer thread). The robustness of the eigenvalues to $T_{\mathrm{samp}}$ is shown in the right panel in figure 11. One can see that a length $T_{\mathrm{samp}}$ of $1\cdot 10^{5}$ time units is more than enough for an estimation on the chosen grid of 400-by-200-400.

B.5 Numerical eigenvalue problem

We have seen that the fine-grained geometry of the eigenvectors of the Lorenz flow require a large grid resolution in order to achieve numerical convergence, if only for the leading eigenvalues. While, the estimation of the transition matrices from time series is, an embarrassingly parallel problem, which can be easily distributed on several nodes of a calculator, solving the eigenproblem for such large transition matrices is more challenging, both in terms of computations and memory. On the other hand, the sparse structure of the transition matrices allows to use iterative eigensolvers. Here, we have chosen the block Krylov Schur algorithm implemented in the Anasazi package Baker2009 of the Trilinos library. This algorithm is a common and robust choice for such problems and the Anasazi implementation makes it straightforward to distribute it on several nodes of a calculator.

Acknowledgements.

AT would like to thank Sebastian Schubert for very helpful discussions on Lyapunov exponents and vectors and for providing his numerical code to estimate them. AT would also like to thank Andrey Gritsun for valuable discussions on the role of periodic orbits in chaotic systems as well as their numerical continuation. AT and VL acknowledge support from the DFG Project MERCI.

VL would like to thank C. Liverani, T. Kuna, and V. Baladi for various inspiring conversations. VL acknowledges the support of the Horizon2020 Project CRESCENDO and of the DFG SFB/Transregio project TRR181.

HD likes to acknowledge the support of the Netherlands Center for Earth System Science (NESSC) funded by the Netherlands Foundation for Scientific Research (NWO). Part of this work was carried out on the Dutch national e-infrastructure (Cartesius) with the support of SURF Cooperative under the project SH284.