Alignment of Lyapunov Vectors: A Quantitative Criterion to Predict Catastrophes?

We argue that the alignment of Lyapunov vectors provides a quantitative criterion to predict catastrophes, i.e. the imminence of large-amplitude events in chaotic time-series of observables generated by sets of ordinary differential equations. Explicit predictions are reported for a Rössler oscillator and for a semiconductor laser with optoelectronic feedback.

deviation. And there are also papers where one looks to the deviation of the statistics of the pulse amplitude as a parameter varies. For a review, see the collection of review papers presented by Akhmediev et al. 3 . Here, however, we are not concerned with extreme events but, instead, with large-amplitude pulses commonly observed in experiments and analysed using small sets of ordinary differential equations.

Results
As described in Methods, the angles between LVs are calculated from the expression 20 Our illustrative examples below involve chaotic dynamics and three Lyapunov exponents λ  ordered such that λ 1 > 0, λ 2 = 0 and λ 3 < 0. One has the following situations for the angles between LVs: (i) θ 13 is the angle between the unstable and stable manifolds, (ii) θ 23 is the angle between the stable manifold and the flow direction, (iii) θ 12 is the angle between the unstable manifold and the flow direction.
Catastrophes in Rössler's oscillator. While studying flows in three dimensions, Rössler was interested in reinjection mechanism where a trajectory, after having a slow motion in one part of a manifold in the phase space, jumps through large excursions to another branch of the manifold. A prototype set of equations with the minimum ingredients for producing sequences of such large excursions is 28 : Here, (x, y, z) are real variables evolving continuously as a function of time, (   x y z , , ) are the corresponding velocities, and (a, b, c) are real parameters controlling the oscillator. These equations give rise to oscillations of x and y which are amplified for a > 0 and result in a spiraling-out motion. Such outbound oscillations are coupled to the z variable in a nonlinear way which induces the reinjection back to the center region, with the cycle of spiraling-out and reinjection repeating indefinitely. Figure 1 displays a typical chaotic trajectory, obtained for a = 0.38, b = 0.3, and c = 4.82 and initial condition near the unstable fixed point x = y = z = 0. Colors represent the values of the aforementioned angle θ 23 . Parallel to the xy plane there is spiraling-out motion, represented in green (θ 23 ~ π/2), while the reinjection loops are shown in black and red, associated with θ 23 approaching either 0 or π, respectively. The attractor in Fig. 1 was obtained solving Eq. (2) numerically with the standard fourth-order Runge-Kutta algorithm and fixed time-step h = 0.05. The first 10 7 steps were discarded as transient, with LVs computed subsequently for 8 × 10 5 steps (including forward and backward time motion). The chaotic trajectory shown is of screw type 28 with irregular oscillations in amplitude and in reinjection times.
The time-series of the representative trajectory shown in Fig. 1 consists of large pulses along the z-axis (associated with a homoclinic orbit 30 ) and are separated by irregular time intervals during which the trajectory has chaotic oscillations basically confined to the (x, y) plane. The large peaks in the z-axis start just before the trajectory gets reinjected back near the origin (i.e. near to the unstable fixed point x = y = z = 0). Our purpose is to show that using the angle between the LVs it is possible to predict when large peaks in the z-axis occur. From Fig. 1 one sees clearly that θ π  23 well before the large peaks emerge out of the (x, y) plane. Figure 2 illustrates the relation between the LVs and the emergence of large spikes. The black curve shows the temporal evolution of z, the red line displays θ 23 , and the blue line depicts the variation of the derivative θ = θ  d dt 23 23 . For easier visualization, we plot θ −  23 . In this range, z presents two large peaks (with z > 10, marked P) and two intermediate peaks (with z ~ 5, marked NP). The blue curve shows that for all four peaks the quantity θ  23 has a local minimum slightly before the peaks. The beginning of two such minima is indicated by vertical green lines. The key here is to observe that before large peaks in z, fast variations of θ  23 occur when, in the red curve, θ 23 rapidly changes either from π/2 → 0 (leftmost peak) or from π/2 → π (rightmost peak).
The above observations suggest the interplay of two conditions prior to the occurrence of spikes. First, an alert condition occurs when the presence of a minimum in θ  23 signals to an imminent spike but without indication of its amplitude. A pair of alert conditions are indicated by the vertical green lines in Fig. 2 (using θ >  1 23 ). Second, an effective catastrophe condition is detected by inspecting how close θ 23 approaches either 0 or π, defined by two suitable thresholds θ 23 (min) and θ 23 (max) . This catastrophe condition corresponds to LV alignment condition, i.e. when the stable manifold aligns along the direction of flow thereby providing information concerning the intensity of the imminent spike. The detection of a pair of catastrophes is indicated by vertical magenta lines in Fig. 2  ). Accordingly, peaks labeled P (predicted catastrophes) satisfy both alert and catastrophe conditions, while peaks labeled NP (non-predicted catastrophes) satisfy only the alert condition. Of course, the discrimination between P and NP can be controlled by tuning the aforementioned thresholds.
To quantify the relation between the LV alignment condition and the intensity of the spikes, Fig. 3 displays local maxima z > 1 as a function of the θ 23 maxima (close to π) and minima (close to 0) for Rössler's oscillator (Fig. 2). The two vertical red lines mark the thresholds θ 23 (min) and θ 23 (max) . To obtain Fig. 3 we evolved the trajectory for a long time checking when maxima (minima) of θ 23 were larger (smaller) than 1.9 (1.2). After this, we computed the maximum (minimum) of θ 23 until a peak in z appeared. In other words, we determined how close θ 23 approaches 0 (or π) near peaks with z > 1. From Fig. 3 one sees that all peaks with large z values are associated with angles θ 23 close to either 0 or π.
To assess the performance of the alignment conditions over extended intervals of time the trajectory was integrated much longer, for t = 5 × 10 4 time steps. The LV alignment was checked automatically using a pair of threshold values θ = . 0 1 23 (min) and θ = .
3 04 23 (max) . In this way, a total of 3745 peaks were found to obey the alert condition θ >  1 23 , i.e. peaks which are potentially catastrophes. Of this total, 2372 where indeed found to be P (predicted) catastrophes, namely z max > 9.11, the value corresponding to the crossing in Fig. 3 of the red lines with the min and max curves. The remaining 1373 spikes were smaller NP (non predicted) events. Spikes in the classes P and NP were separated by imposing numerically the alarm and LV alignment conditions. Thus, all detected spikes were correctly accounted for. In addition, on average we determined alert times to precede spikes by ~2.2 time units, and catastrophe times to precede by ~0.85 time units. Of course, much earlier prediction times are obtained using  the catastrophe prediction time increases to 1.51, almost doubles, on the average. In this case 2944 peaks of type P were detected. We remark that the large interval of integration chosen above was used simply to check that the method was able to detect the imminence of large pulses to get a feeling for the distribution of large events in the model considered. Large integration intervals are by no means needed to predict catastrophes.

Catastrophes in semiconductor laser with optoelectronic feedback. Trajectories for Rössler's equa-
tions are known to be closely related to the presence of Shilnikov homoclinic orbits in the system [28][29][30] . To check a possible influence of homoclinicity on the prediction times, in the next example we consider a situation where Shilnikov's theorem does not apply 31,32 , namely we consider a semiconductor laser with optoelectronic feedback. The aim here is to show that the alignment of LVs is able to predict large peaks also in absence of homoclinicity. In dimensionless units the laser is governed by the equations 31,32 Here, the variables x, y, w are proportional to the photon density, the carrier density, and the laser intensity, respectively, while    x y w , , are corresponding rate of change. The feedback is controlled by f = α(w + x)/ (1 + s(w + x)), γ is the ratio between the population relaxation rate and the photon detuning, ε is the high-pass frequency in the feedback loop, δ 0 is the solitary laser threshold, s is a saturation coefficient of the amplifier and α is a coefficient proportional to the photodetector responsivity. In our simulations we fix α = 1, γ = 1 × 10 −3 , ε = 2 × 10 −5 , δ 0 = 1.017 and s = 11, values for which chaotic dynamics is observed experimentally 32 .
As before, Eqs (3)-(5) were solved numerically with the fourth-order Runge-Kutta algorithm and fixed time-step h = 0.05. The first 10 7 time-steps were discarded as transient and the LVs were computed for 9.5 × 10 5 times-steps (including forward and backward time motion). Figure 4 shows a typical trajectory with colors representing the values of θ 23 for every point of the trajectory. We started from (x 0 , y 0 , w 0 ) = (1, 1.5, 0.5) and observed a small spiraling-out motion in the (x, y) coordinates, then an abrupt increase of w, and then a larger spiraling-in motion leading to the reinjection process. It is easy to realize that when w increases, θ 23 changes colors and assumes values approaching 0 or π. This can be more easily observed in the time evolution of w shown in Fig. 5 (in units of 10 4 time intervals). Figure 5(a) displays the scaled variable w′ = w × 100 and Fig. 5(b), for the same times, presents the θ 12 (red line) and θ 23 (blue line). Comparing Fig. 5(a) and (b) we observe that preceding all large peaks in w, one finds θ 12 → 0 and, simultaneously, θ 23 → π. Note that in contrast to the Rössler's system, in this case both stable and unstable manifolds align along the direction of flow. Thus, once again, the LVs from the stable and unstable manifolds tend to align along the flow direction before large-amplitude events. In other words, large peaks of w occur shortly after the alignment of stable and unstable manifold along the flow direction. As seen in Fig. 5(b), in this case both the alert and catastrophe times (i.e. prediction times) can be close to each other. On average, prediction of large peaks occurs at about 1 time unit (10 4 time-steps).

Discussion
The alignment of LVs along the flow direction provides a straightforward and apparently reliable means of predicting catastrophic events in chaotic dynamical systems. This was corroborated for two familiar dynamical systems, namely a Rössler oscillator, and a semiconductor laser with electronic feedback. We are not aware of any other quantitative criterion able to predict large-amplitude events.
To predict large-amplitude events, one should evolve a trajectory and monitor θ  23 . Whenever its magnitude becomes larger than a suitable system-dependent threshold, a peak should be expected. To estimate the peak intensity, evolve the trajectory further and check how close θ 23 approaches 0 or π. The closeness of the approach defines the peak intensity and, therefore, can anticipate a catastrophe. In realistic situations where noise is present, the validity of the alignment conditions depends on the survival of the spikes under the noise influence.
An open challenge is to check the effectiveness of LVs alignment to predict catastrophes for higher-dimensional systems, a considerably more complicated framework for which there are yet no efficient methods to compute angles between invariant manifolds. Nevertheless, from a theoretical point of view we anticipate no problems in predicting catastrophes regardless of the underlying dimensionality of the system. This, of course, remains to be explored.

Methods
Here, we describe briefly the essential details for the numerical implementation of the procedure to compute the angle between the two invariant subspaces of a 3D dissipative dynamical system. We are interested in the forward evolution in the tangent space of an initially orthonormal basis = G a a a ( , (3) . Such evolution is governed by the equation where J n = J n−1 J n−2 … J 0 is the product of the Jacobian matrix of the map evaluated for every orbital point. It is implicit here that at the initial time n = 0 we have already computed the orbit for a sufficiently long time (forward transient time) so that the orthonormal basis of the tangent space converged already to the asymptotic Gram-Schmidt (GS) vectors 17 . In order to avoid divergences the matrix ∼ G n is renormalized for each step n. As usual, this can be done by the QR decomposition = .  and contains the information obtained in the GS orthonormalization procedure of ∼ G n . Here, = G g g g ( , , ) n n n n are the GS vectors after the orthonormalization. Since g g , n n (1) (2) and g n (3) are orthonormal by construction, they can only provide information about the local rates of expansion (contraction) of the vectors, stored in the diagonal elements of R n , from which the standard Lyapunov exponents can be obtained. However, to calculate the angle between the invariant subspaces in the tangent space, we need the LVs determined from the relations  LVs have normalized length so that the columns of → C j must be normalized to 1. The initial condition for C n , before starting the backward evolution, can be a generic nonsingular upper triangular matrix, which is the GS basis.