Spectral Early-Warning Signals for Sudden Changes in Time-Dependent Flow Patterns

Lagrangian coherent sets are known to crucially determine transport and mixing processes in non-autonomous flows. Prominent examples include vortices and jets in geophysical fluid flows. Coherent sets can be identified computationally by a probabilistic transfer-operator-based approach within a set-oriented numerical framework. Here, we study sudden changes in flow patterns that correspond to bifurcations of coherent sets. Significant changes in the spectral properties of a numerical transfer operator are heuristically related to critical events in the phase space of a time-dependent system. The transfer operator approach is applied to different example systems of increasing complexity. In particular, we study the 2002 splitting event of the Antarctic polar vortex.


Introduction
Characterizing and predicting the occurrence of sudden changes in the dynamics of complex systems-typified in this context by geophysical flow patterns-is crucial in understanding and controlling such critical events, which may have serious environmental consequences. Famous examples include the 2010 oil spill in the Gulf of Mexico and the sudden split of the Antarctic ozone hole in September 2002 (see Figure 1).  25 September, 2002. Velocity data from the ECMWF Interim data set (http://apps.ecmwf.int/datasets/). Such complex systems may be modeled as incompressible time-dependent dynamical systems where the sudden critical changes are characterized as critical transitions [1][2][3] of dominant slowly mixing patterns. The identification and characterization of such coherent flow patterns in the Lagrangian frame of reference has received an increasing interest in the last two decades. Many different identification methods have been developed and applied to fluid flow (see [4][5][6] for recent reviews). In this context, a coherent set [4,[7][8][9][10] is defined as a region in the fluid volume that only weakly mixes with its surroundings through a timedependent flow. Coherent sets were originally introduced based on transfer operators [7,8] and can be seen as non-autonomous, finite-time analogues of almost-invariant sets [11] (see [12] for a thorough discussion of the different constructions). A transfer operator is a linear Markovian operator that describes the evolution of densities of the underlying nonlinear dynamical system [13]. Within a set-oriented numerical framework [14], this infinite-dimensional operator can be approximated by the transition matrix of a finitestate Markov chain, and the sign structures of left and right singular vectors to singular values close to 1 are systematically used to yield initial and final time-robust coherent patterns [7,8,12]. We note that this is fundamentally different from the singular value decomposition methods often used in atmospheric sciences (see, e.g., [15] (chp. 6.3), [16]). In those frameworks, the linearized propagator of the underlying partial differential equation is investigated, whereas the respective transfer operator approach is based on a probabilistic description of the nonlinear and time-dependent Lagrangian particle dynamics in physical space.
In this work, we want to study how sudden changes (in the sense of happening on relatively short time scales) in time-dependent flow patterns relate to the behavior of the corresponding coherent sets and, thus, to the spectral properties of the numerical transfer operator. We first consider some very simple example systems before we study flow data from the Antarctic polar vortex splitting event in September 2002.
Early analytical results relating the transfer operator spectrum to classical bifurcations can be found in [17]. In [18], early-warning indicators for transitions between atmospheric flow regimes were defined based on the transfer operator of a dissipative autonomous atmospheric model. Bifurcations of almost-invariant and almost-cyclic sets in two-dimensional autonomous conservative systems and corresponding changes in the spectra of the transition matrices were observed in [19,20], but were not systematically studied. In [21,22], the bifurcation and splitting of almost-invariant sets, which are interpreted as dominant patterns in autonomous incompressible flows, were considered. There, almost-invariant sets were built around a global stationary state, whose bifurcation yielded the qualitative change of the underlying flow patterns. All of these previous works dealt with autonomous systems, and only very recently has a transfer-operator-based computational approach for studying changes in non-autonomous flow structures been proposed [23]. In that work, a numerical framework was introduced, which was mainly applied to (quasi-)periodically driven systems, but also to polar vortex flow data. Singular values and vectors, as well as their equivariance properties, were studied. In this paper, we focus on the non-autonomous transition phase between different flow regimes with the aim of obtaining a better understanding of the relation between phase space structures and spectra. We hope that these case studies can serve as a starting point for theoretical and practical advances in the field. This paper is organized as follows: In Section 2, we review the set-oriented framework for studying finite-time dynamics and coherent sets. Early-warning signals for sudden changes in coherent flow patterns are investigated in Section 3 by first considering a onedimensional non-autonomous pitchfork bifurcation in finite time, and then two explicit non-autonomous incompressible two-dimensional ordinary differential equations where flow patterns change in time in a controlled manner. In Section 4, we relate these findings to the behavior of coherent sets during the Antarctic polar vortex splitting scenario in September 2002.

Background on Non-Autonomous Dynamics and Coherent Sets
In this section, we describe a set-oriented framework for non-autonomous dynamical systems that allows one to study finite-time transport of macroscopic non-autonomous objects.

Non-Autonomous Dynamical Systems
We consider a non-autonomous ordinary differential equation: whereẋ denotes the derivative of x with respect to t. We assume that f : R × R d → R d satisfies conditions for global existence and uniqueness of solutions, and we denote the general solution or process of (1) by φ(t, t 0 , x 0 ), i.e., φ(·, t 0 , x 0 ) solves (1) for the initial condition x(t 0 ) = x 0 . The mapping φ fulfills the initial value property φ(t 0 , t 0 , x 0 ) = x 0 and the cocycle property [24,25] We are interested in the finite-time evolution of macroscopic objects and restrict the analysis of (1) to a finite time interval [t, t + τ] of length τ > 0 and a compact set X t ⊂ R d at a fixed initial time t ∈ R. We denote the image state space at final time t + τ by Y t+τ = φ(t + τ, t, X t ).

Coherent Sets
We are interested in identifying coherent sets, which are regions of the phase space that are minimally dispersive or maximally coherent over the finite time interval [t, t + τ] under consideration. Let µ t be a probability measure describing the initial mass distribution at the initial time t supported on X t . We assume that µ t is absolutely continuous with respect to the Lebesgue measure (volume measure), and thus possesses a density. The process φ maps the measure µ t to a measure ν t+τ at time t + τ via ν t+τ (A) = µ t (φ(t, t + τ, A)) for all measurable sets A ⊂ Y t+τ . The measure ν t+τ represents the mass distribution of objects of interest on Y t+τ at time t + τ. In the case of incompressible flows, both µ t and ν t+τ can be considered as (normalized) volume.
We note that any pair of measurable sets (A t , A t+τ ), where A t+τ = φ(t + τ, t, A t ), is perfectly coherent in the sense that its coherence ratio [7], given by Here, r measures the proportion of A t that is mapped to A t+τ under the action of the flow map φ(t + τ, t, ·). In the literature on nonautonomous dynamical systems, the sets A t and A t+τ with r(A t , A t+τ ) = 1 are denoted as fibers of an invariant non-autonomous set [25], and since for any given measurable set A t , a corresponding set A t+τ can always be found by evolving this set to time t + τ, the problem of finding coherent sets, that is, sets that maximize (2), is ill-posed in our setting. However, we are interested in coherent sets that are robust with respect to perturbations. While including diffusion regularizes the mathematical problem of finding optimal sets [8], it also makes it applicable to real-world systems that are naturally perturbed by noise. We note that in this perturbed setting, the above-defined coherence ratio will be close to 1, but not equal to 1 in general. We primarily seek coherent sets that have the following properties: (i) Maximizing the coherence ratio r as given in (2), (ii) Robustness under external perturbations, (iii) Closeness to a critical point or orbit that changes stability (this is crucial for setting up set-oriented analogues to finite-time bifurcations), (iv) Minimization of the surface-to-volume ratio (in conservative flows).
Such optimization problems can then be approximately solved by considering the Perron-Frobenius operator P t,τ : L 1 (M, ) → L 1 (M, ) associated with the flow map φ, where denotes the Lebesgue measure and M ⊂ R d . The transfer operator is defined by The interpretation is that if f is a density and f (x) is the density value in x at time t, then P t,τ f (x) describes the density value in φ(t + τ, t, x) at time t + τ induced by the flow map. P t,τ is a Markov operator, i.e., a linear operator whose spectrum is contained in the unit circle, which pushes forward densities under the action of the underlying nonlinear and time-dependent flow. In [8,12], it was shown that maximizing r in (2) can be described in the framework of optimizing an inner product involving a compact self-adjoint Markov operator Q obtained from a perturbation of P t,τ . In particular, it was shown in [8] that where λ 2 < 1 is the second largest eigenvalue of Q, and A c t = X t \ A t and A c t+τ = Y t+τ \ A t+τ are the complements of A t and A t+τ , respectively. That means that the existence of optimal coherent sets is directly tied to the existence of an eigenvalue λ 2 close to one. The supports of the positive and negative parts of the corresponding second eigenfunctions approximate (A t , A c t ). This also extends to families of k coherent pairs [26], where further eigenvalues close to 1 and the corresponding eigenfunctions have to be considered.
However, in order to avoid the technical functional analytic description underlying these ideas, we briefly review the concept of finite-time coherent sets in the finitary setting [7], where we work with finite-rank approximations of the infinite-dimensional operators.

Set-Oriented Framework
We now discretize the state spaces X t and Y t+τ . We subdivide X t and Y t+τ into partitions {B 1 , B 2 , . . . , B m } and {C 1 , C 2 , . . . , C n } containing a finite number sets, respectively, i.e., X t = m i=1 B i and Y t+τ = n i=1 C i . The partition elements (also called boxes) of X t satisfy where denotes the Lebesgue measure, and we assume the same for the partition elements of Y t+τ . We define the initial and final time-lumped finite state spaces [27] as respectively, and we define the finite-time transition matrix P t,τ [7] via where φ is the process induced by (1). Note that P t,τ is an m × n row stochastic matrix, and each entry P t,τ ij is the probability that a randomly chosen point in B i at time t will be mapped into the box C j at time t + τ. The matrix P t,τ is an Ulam finite rank approximation [28] of the Perron-Frobenius operator [7,13], generated by the underlying non-autonomous dynamical system (1), and naturally includes the required perturbation for studying coherent sets due to it being based on the finite partitions of X t and Y t+τ . In this paper, the numerical estimation of the entries in Equation (5) is given by [7] which is the proportion of the uniformly distributed initial test points (particles) {x k } N k=1 in box B i at time t that are mapped to box C j at time t + τ. We implemented this using the MATLAB-based software package GAIO [14]. The estimation of a transition matrix from particle trajectories is computationally expensive and depends on the number of boxes in the discretization, the number of test points, the flow time, the time step size, and other parameters of the chosen numerical solver for ordinary differential equations. However, the integration of particle trajectories can be done in parallel.

Singular Vectors and Coherent Partitions
We are interested in finding coherent sets using the sparse transition matrix P := P t,τ defined in (5). We consider an initial probability distribution p ∈ (0, 1) m defined via p i := µ t (B i ) for i ∈ {1, . . . , m}, so p is the probability row vector that approximates the initial distribution of mass µ t . The image probability vector q = pP represents the distribution of mass at final time t + τ. We assume that p > 0 and q > 0 (component-wise). We say that the pair (A t , A t+τ ) ∈ B m × C n forms a pair of coherent sets over the finite time interval [t, t + τ] if the following is satisfied: There exist index sets I ⊂ {1, 2, . . . , m} with A t = i∈I B i and J ⊂ {1, 2, . . . , n} with A t+τ = j∈J C j such that We form a transition matrix L from P via The matrix L is a normalization of P and, as P itself, depends both on the initial time t and the final time t + τ. It represents the transition matrix of the macroscopic dynamics with respect to the initial and final mass distributions p and q. L satisfies 1 m L = 1 n [8], where 1 m , 1 n are the vectors in R m and R n consisting of 1 in each component, respectively.
In [7,8], it was shown that (under some technical assumptions) the problem of finding optimal coherent sets can be approximated by considering the left eigenvectors u 2 ∈ R m of LL * and v 2 ∈ R n of L * L to the second largest eigenvalue λ 2 < 1. Here, L * = P . LL * and L * L are approximations of the compact self-adjoint operator Q briefly mentioned earlier and its dual. Note that these two eigenvalue problems can be turned into the task of finding leading singular values and corresponding left and right singular vectors of a sparse matrix (see [7] for the exact construction), which can be very efficiently computed by iterative schemes (e.g., svds in MATLAB).
The signed vector entries of u 2 and v 2 can be interpreted as relaxations of indicator functions of the sets A t and A t+τ and their complements. Thus, the vector u 2 defines fuzzy coherent sets on X t , whereas v 2 represents their image on Y t+τ . A partition of X t and Y t+τ into finite-time coherent pairs (A 1 t , A 1 t+τ ) and (A 2 t , A 2 t+τ ) can be approximated by considering the positive and negative level sets of u 2 and v 2 , respectively: and we can also pick (A 2 t , A 2 t+τ ) as the negative level sets of u 2 , v 2 , or, alternatively, define Thus, using singular vectors to the singular value σ 2 , we have constructed two complementary pairs of coherent sets.
In [7,12], a line search was used to find an optimal partition from u 2 , v 2 that maximizes (7). However, this approach is restricted to finding a partition into two sets: of a coherent set and its complement. In practice, there are often k > 2 singular values σ i , i ∈ {1, . . . , k}, close to 1 (ideally followed by a spectral gap), and the corresponding singular vectors highlight the location of coherent sets, which could be extracted using a hard-clustering approach, such as k-means. Alternatively, to preserve the eigenspace structure, one can project the singular vectors to a sparse basis [29], where the entries of each vector denote probabilities that the underlying box B i belongs to a specific coherent set.
The question of how many singular values and vectors are to be considered is a subtle one, as it is in other reduction techniques, such as principal component or Fourier analyses.
The singular values provide a ranking of the importance of the flow patterns (represented by the corresponding singular vectors) in terms of their coherence. As outlined above, the second singular vectors provide the most coherent partition into two sets, with further vectors typically picking up increasingly smaller structures. In some systems, one finds k leading singular values (i.e., close to one), followed by a clear spectral gap, indicating that the dominant coherent dynamical features are represented by the respective k singular vectors. When there is no clear gap in the spectrum, how many vectors are considered is the user's choice, and this can be considered as a limitation of the method.
We note that in the context of this paper, we do not postprocess the singular vectors, but consider them as representations of the underlying flow patterns. Finally, we note that the singular values and vectors and, thus, the flow structures of interest are robust with respect to perturbations by construction due to being obtained from a finite-rank approximation of the transfer operator that includes numerical diffusion in the order of magnitude of the box sizes. Robustness can also be shown analytically by using results from the perturbation theory of reversible Markov chains [30,31] and has also been studied numerically in example systems [32].

Changes in Spectrum and Flow Patterns
The aim of this work is to study changes in flow patterns under the action of a non-autonomous dynamical system in case studies and, in particular, to identify spectral signatures of splitting or significant changes in the shapes of coherent sets. In [8,12], several properties of the underlying transfer operator were discussed, including the existence of lower bounds in (4), as well as the regularity of the leading singular vectors and its dependence on the underlying small random perturbation. For analytical results on the geometric features of coherent sets, we refer to [9]. We state a few heuristic arguments that may help to relate changes in singular values and changes in coherent flow patterns.

1.
By construction, coherent sets are also characterized by minimizing diffusive transport, and this puts constraints on their shape [8,9]. Coherent sets that have a long boundary compared to volume would be less coherent than those of the same volume with a smaller boundary. Thus, when a coherent set is reshaped, e.g., from a ball to an ellipse of the same volume, we would expect a decrease in the corresponding singular value. This is due to the fact that the singular value directly enters the upper (and also lower) bounds on the coherence ratio.

2.
For the same reason, the size of a coherent set influences its coherence ratio, since this is a relative quantity. Large coherent sets are, in general, more coherent than smaller ones, as diffusive transport affects larger parts of the latter. Thus, when a coherent set grows in size, we would expect an increase in the corresponding singular value.

3.
Finally, when a coherent pattern becomes more dominant, e.g., when a coherent set splits, this might change the relative order of coherent patterns in terms of their coherent ratios. In this case, we would expect a crossing of the corresponding singular values, as also observed in [23]. For the related case of almost-invariant sets in autonomous dynamical systems, this has been studied in [19][20][21][22]. Computational approaches for the detection of singular value crossings are proposed in [23].
A clear practical limitation of the transfer operator studies is that, in general, one cannot expect a certain magnitude of the spectral changes when approaching and crossing a critical event, but rather, trends in the spectrum have to be considered. The detailed behavior depends on the system itself, the choice of the domain and flow time, and the level of (numerical) diffusion, and we note that this specific dependence is also an issue with early-warning signals (such as autocorrelation and variance) that are frequently used as statistical indicators of critical transitions [1][2][3].

Studying Bifurcations of Coherent Sets and Early-Warning Signals
In this section, we study different examples in order to develop a better understanding of the spectral footprints of bifurcations of coherent sets.

Spectral Signatures of a Finite-Time Bifurcation in One Dimension
We start the finite-time bifurcation analysis of coherent sets with the one-dimensional ordinary differential equatioṅ which is a non-autonomous version of a classical pitchfork bifurcation (see [33] (p. 146)).
The process of (10) is denoted by φ, and its finite-time dynamics are relatively simple to grasp. For all t ∈ R, the right-hand side vanishes for x = 0, which means that the differential equation is solved by the trivial solution. This solution is attractive as long as t < 0, and becomes repulsive for t > 0. This change in stability in time can also be described by finite-time notions of attraction and repulsion [34]. We note that for t > 0, the right hand side also vanishes for x = ±w(t), where w(t) = arctan(t), and, when considering time as bifurcation parameter, this corresponds to a supercritical pitchfork bifurcation. Due to the dissipative behavior, there are no coherent sets in our interpretation, as all trajectories will be attracted to the x-axis before the bifurcation and attracted to ±w(t) for t > 0. However, by including a set-valued perturbation and by considering the finite-time dynamics on short time intervals, one can artificially create a splitting of the phase space into two coherent sets when the unperturbed system crosses the bifurcation point. For detailed discussion of the setting and the computational details, we refer to Appendix A, and we only present the main results here. Using a set-oriented discretization, we compute forty transition matrices P t,τ , where t ∈ {−20, −19, . . . , 19} and τ = 1, from which we obtain the corresponding singular values σ t i , as well as the left and right singular vectors u t i and v t i . We refer to this as the sliding window approach. On the left-hand side of Figure 2  We see that both σ t 2 and σ t 3 are significantly smaller than 1 whenever t < 0, but both singular values rise sharply when approaching and passing the (classical) bifurcation point. This indicates that the perturbed system does not admit proper coherent sets for negative times, or, viewed differently, the whole state space is the only coherent set. Due to the emergence of the two attracting solutions after the (classical) bifurcation point, i.e., for positive times, we obtain two coherent sets, and this is clearly visible in the spectrum, with σ t 2 ≈ 1 and a considerable gap to σ t 3 for t > 0. We stress that the rise of the singular values already starts before the bifurcation, and can thus be considered as an early-warning signal. The corresponding changes in the singular vectors are very subtle and will be discussed in detail in Appendix A (using Figures A1 and A2).
Moreover, using the second left and right singular vectors u t 2 and v t 2 , we can numerically approximate the coherence ratio (2), where we used the sign structure of the vectors to define the coherent sets. The coherence ratio is plotted on the right-hand side of Figure 2 for the different time intervals. We note that for small values of t, the coherence values are much smaller than 1, and as a consequence, the second condition of Equation (7) is not fulfilled in this regime, which indicates that there are no coherent sets. On the other hand, when t approaches 0, the coherence ratio increases and reaches a value very close to 1, with the two coherent pairs ( The singular values and coherence ratios in Figure 2, as well as the changes in the right singular vectors when approaching the bifurcation point (see Figure A1g-i), are clear indicators of the finite-time bifurcation at t = 0, and can thus be simultaneously used as early-warning signals for the upcoming bifurcation.
The study of this relatively simple one-dimensional perturbed model demonstrates that variations in the singular values and corresponding singular vectors are relevant observables for understanding finite-time bifurcations, and they serve as early-warning signals in this particular example. In the remainder of the paper, we study several nonautonomous two-dimensional incompressible systems that undergo critical changes in coherent sets in their finite-time evolution. Based on the behavior of the singular values for different time windows [t, t + τ], we study early-warning signals of such critical changes.
The differential Equation (11) is asymptotically autonomous, since the time-dependent bifurcation parameter only varies in the interval (0, 1). As shown in Figure 3, the differential equation models a single rotating vortex for t ≤ 0 and a double rotating vortex for t ≥ 1, and the system is non-autonomous on the interval [0, 1], modeling a transition from a single to a double gyre pattern. As in the previous section, singular values are used to understand the finite-time transition of the vortex pattern, which is determined by the corresponding singular vectors. We note that this transition is different from a proper splitting, but nevertheless, the singular values undergo interesting changes during this transition. Figure 3. The velocity field of (11) at t < 0 (a), t = 0.3 (b), t = 0.4 (c), and t = 0.9 (d).
A time interval [t, t + τ] of length τ = 1 turns out to be useful to illustrate this by using the singular vectors and a sliding window approach with initial times t ≤ 0. The differential Equation (11) is integrated using a fourth-order Runge-Kutta scheme in the spatial domain X = [0, 1] × [0, 1]. This domain is subdivided into 2 12 equally sized boxes, and 100 uniformly distributed test points are initialized in each box. As the domain is invariant (i.e., X = X t = Y t+τ for all t, τ), the initial and final box coverings B m and C n coincide. We set up the transition matrix (8) for each [t, t + 1] and compute the singular value decomposition in order to extract the coherent patterns from the singular vectors and their corresponding singular values. We then explore the spectral changes accompanying the transition. More specifically, starting from a single pattern, i.e., when t + 1 < 0, we track its behavior as time evolves until t = 1. Figure 4 shows the changes in the second singular value σ t 2 of the transition matrices (8) for different time windows [t, t + 1]. We see three different phases: a constant phase, a decreasing phase, and an increasing phase. During the constant phase, the patterns of the left and right singular vectors u t 2 and v t 2 do not change as the time interval moves; they are both indistinguishable and look like the pattern shown in Figure 5a. The corresponding left and right singular vectors of the decreasing phase of the singular values in Figure 4 are shown in Figure 5. We clearly see that the left singular vectors u t 2 (Figure 5a-d) at initial times t and the respective right singular vectors v t 2 (Figure 5e-h) at final times t + 1 are no longer indistinguishable. Indeed, as the second gyre (red color) starts to rise from the right-hand side of the domain at times t + 1, we see that the initial gyre (the blue-colored pattern) starts to shrink (Figure 5e-h), while at times t, it does not (Figure 5a-d). This shrinking process of the initial vortex is captured by the singular values with a decreasing phase (Figure 4). The more the second gyre emerges and increases in size, the more the initial gyre undergoes a fold-like shrink. Finally, the rising singular values in Figure 4 correspond to the singular vectors shown in Figure 6.  The emerging regime of the second gyre no longer causes the shrinking of the first gyre. In fact, both gyres start to increase in size. The finite-time evolution of the system is then similar to the behavior after the bifurcation in the one-dimensional toy model (10). Indeed, the left singular vectors (Figure 6a-d) at time t exhibit two coherent sets, each of which is separately mapped to the coherent sets at times t + 1 (Figure 6e-h).
Although this case study does not represent a proper splitting, the trends of the singular values are still interesting, since they capture the topological changes leading from a single gyre to a double gyre. We proceed with the study of early-warning signals for a proper pattern splitting.

Early-Warning Signals for a Vortex Splitting Regime
We study a non-autonomous Duffing-type oscillator given by the differential equatioṅ An autonomous version of this model was studied in [21,22], where spectral indicators of the bifurcation of almost-invariant sets were investigated.
Here, we explore early-warning signals for a critical transition that occurs when coherent patterns split in a neighborhood of the origin. The latter is the trivial critical point of (12), which undergoes a pitchfork bifurcation when t crosses 0, as illustrated by the velocity fields in Figure 7. This bifurcation was used in [22] to study the splitting of dominant almost-invariant sets supported in a neighborhood of (0, 0). Figure 7. Examples of the velocity field of (12) for t < 0 (left) and for t > 0 (right).
The system is studied in the domain X = X t = [−2, 2] × [−2, 2], which is subdivided into 2 14 square boxes with 100 uniformly distributed test points in each box. This defines our initial covering B m . The time interval under consideration is [−15, 15]. We vary t from −15 to 5 in steps of 0.4 and compute 51 transition matrices for time intervals [t, t + τ], where τ = 10 and our model (12) are integrated with a fourth-order Runge-Kutta scheme. For the final covering C n of Y t+τ , we consider the image of X under the corresponding flow map on [t, t + τ] and discretize it with boxes of the same size. For all matrices, we compute the singular value decomposition in order to extract coherent patterns from the singular vectors. The first ten dominant singular values, plotted with respect to the final times t + τ, are presented on the left-hand side of Figure 8. We restrict our analysis to the first eight singular values and vectors in order to focus on the dominant flow patterns. This choice can also be justified by the spectral gap between σ t 8 and σ t 9 for negative final times t + τ. Some of the corresponding left and right eigenvectors are presented in Figure 9, where the upper panel shows the left singular vectors u t 2 , . . . , u t 8 (u t 1 , v t 1 are constant by construction) at initial times t indicated below the plots, and the lower panel shows the corresponding right singular vectors v t 2 , . . . , v t 8 at final times t + 10 given below the plots. Note that in order to focus on the pattern changes, we show the right singular vectors restricted to the original domain X. Y t+τ extends out of X and changes in time. Moreover, we used a uniform color scheme, and although vector entries may take larger or smaller values than in the interval [−2, 2], we restrict to this interval.
For t + τ < 0, we have nearly constant singular values, and the corresponding singular vector patterns also do not appear to change much, as seen in the first three columns of Figure 9. Moreover, both the second and third singular values remain constant until about t + τ = 10. Then, they start to decrease slightly (and later become constant again, which is not shown) as the system has entirely crossed the bifurcation point with t > 0 and evolves asymptotically into an autonomous system. The corresponding left and right singular vectors shown in the first two rows of the upper and lower panels in Figure 9 change slowly, but the highlighted sets in u t 3 and v t 3 stay topologically the same. Other interesting spectral events include crossings of singular values. When t + τ ≈ 0, the singular values σ t 5 and σ t 6 become very close and the curves cross back and forth (C), as can also be seen in the corresponding plots of the singular vectors, which look very similar in that regime, indicating rotational dynamics. v t 5 continues to highlight rotation, but outside of the central coherent set, whereas v t 6 highlights rotation within the bifurcating coherent set, which gets blocked when the set starts to split. The singular value curves cross once more (D) at the point where the coherent set has entirely split into two pieces, and the singular vectors u t 5 and v t 5 highlight that we have obtained two different coherent sets. Further crossings (E,F) indicate that the rotational pattern becomes less dominant.
We note that the singular value patterns correspond to real, optimally coherent flow structures with respect to the time interval under consideration. They are representative of the motion of bulks of particles. In Figure 10, we demonstrate that groups of Lagrangian particle trajectories experience the same changes as the flow patterns.
In summary, in this simple model, the splitting of a vortex-like coherent flow pattern is preceded by an elongation of the coherent set, and this leaves its spectral footprint in decreasing singular values as expected. This is accompanied by several crossings of singular values. In the next section, we find similar spectral indicators in a concrete real-world application of pattern splitting. Version November 18, 2020 submitted to Fluids 21 of 21

Spectral Study of the Sudden Antarctic Polar Vortex Split of September 2002 from Recorded Satellite Data
In September 2002, an unusual critical transition occurred in the Antarctic polar vortex region: The ozone hole suddenly split into two rotating blobs of air as a consequence of the stratospheric warming. This was the first event of this kind that was recorded in the meteorological history of the Southern Hemisphere, and many scientific studies [35][36][37] have focused on understanding the main causes of the splitting of the ozone layer. In [36], the vortex dynamics are associated with the evolution of chemical constituents, such as ozone O 3 and chlorine monoxide (ClO), before, during, and after the split. Furthermore, based on satellite velocity data, it is very well understood that the split occurred during the last week of September 2002, between 24 September and 27 September [37].
In the following, we study the sudden split by analyzing the finite-time global dynamics of the vortex from recorded velocity data. Using the available velocity data from the ECMWF Interim data set (http://apps.ecmwf.int/datasets/), we track the changes of the Antarctic polar vortex, as illustrated earlier in Figure 1. Its different geometric shapes during the splitting event of September 2002 were studied in [38] by means of finite-time Lyapunov exponents.
The global ECMWF data are given at a temporal resolution of 6 h and a spatial resolution of a 0.5 • in longitude and latitude degrees, respectively. We focus on the stratosphere over the Southern Hemisphere and consider the velocity data from 1 September 2002 to 17 October 2002 on a 600 K isentropic surface. The restriction to two-dimensional dynamics can be justified by the observation that isentropic surfaces are nearly invariant over short times of about two weeks [39], and thus, particle motion is predominantly confined to a two-dimensional surface. Thus, the equation of motion of the atmospheric dynamics in that region is interpreted as a two-dimensional non-autonomous ordinary different equation: determined on discrete points (t, x, y). Using linear interpolation in space and time, we can integrate solutions again by a fourth-order Runge-Kutta scheme and obtain an approximate flow map. For the set-oriented analysis of the global dynamics, we consider the region south of 45 • S, discretized by m = 8920 square boxes {B 1 , . . . , B m } of side length 93.75 km. We cover the image domain with a corresponding collection of boxes {C 1 , . . . , C n } (with varying n), where the C i have the same size as the initial boxes. We choose n = 49 sample points uniformly distributed in each box and compute the transition matrices as in (8). The polar vortex at different times is then considered as a coherent set approximated by the singular vectors of (8), and the singular values are used to identify possible early-warning signals.
Following previous studies [7,23], we use a flow time of τ = 14 days, in accordance with the validity of the two-dimensional restriction to the isentropic surface. We compute the transition matrices on [t, t + τ]. Our first initial time t is 0:00 on 1 September 2002, and we increase t in increments of 6 h. Figure 11 shows the leading ten singular values plotted with respect to the final times between mid-September and the end of October 2002. One clearly sees that between 20 and 30 September, the spectrum fluctuates very strongly. This is exactly the period where the Duffing-type bifurcation of the polar vortex happens, starting with a weakening of the vortex on 21 September and the splitting process from 24-26 September [37]. One can observe many crossings of singular values in that regime, but these are very hard to analyze because they often correspond with structural changes outside the vortex. The polar vortex is surrounded by the polar jet stream, and there are other smaller vortices that enter and leave our domain of interest.
As the second singular vectors are known to represent the polar vortex flow pattern during the split [23], we focus on the behavior of the corresponding singular value σ t 2 for final times between 18 September and 4 October 2002 (Figure 11 (inset)). One observes separate particular behaviors: nearly constant singular values from 16 September to 20 September and decreasing singular values from 23 September to 25 September, followed by a rising and then approximately constant regime.  In the previous examples, we observed that nearly constant singular values yield corresponding singular vectors whose support or level sets do not change in size, corresponding to a stable regime. The Antarctic polar vortex appears to be stable from 16 September to 20 September, and then starts to weaken. From 23 September to 24 September, the singular values in Figure 11 begin to decrease considerably. The corresponding left and right singular vectors are shown in Figure 13.
The shape changes of the left and right singular vectors' patterns are very similar to those in the high-order non-autonomous Duffing oscillator in Figure 9, including the decrease in magnitude of the corresponding singular value. Thus, a spectral early-warning signal of the split occurs from 23 September to 24 September, where we start to notice a breaking up of the rotating vortex into two blobs of air, as shown in Figure 13.
The post-splitting scenario begins on 27 September. From this date onwards, the northwestern blob of the ozone hole starts to disappear, while the southwestern blob begins to increase again in size to recover the initial shape before splitting. This can be classified as a reformation phase. The leading singular vectors continue to capture the polar vortex as shown in Figure 14.
Finally, beyond the second dominant singular vectors and corresponding singular values, it may be crucial to explore pattern structures of more singular vectors. We restrict our study to the polar vortex region itself, thereby excluding the polar jet stream and other coherent structures that are otherwise highlighted in subdominant singular vectors. In Figure 15, the second, third, and fourth right singular vectors' patterns exhibit, as expected, a finer partition of the polar vortex. While the second singular vectors v t 2 in Figure 15a-c continue to highlight the elongation and splitting of the vortex, the third and the fourth vectors, v t 3 and v t 4 , may help to explain the splitting scenario of the polar vortex from a local point of view with a net symmetric separation of the finite-time dynamics of the local flow into two parts (blue and yellow colors), where each part lies in one blob of air (see Figure 15f-h). We note that the discussion of the singular values and vectors at final times t + τ is effectively based on using past velocity data only, i.e., up to τ = 14 days to the date under consideration. Forecasts could be included in the computation accordingly (e.g., by using 10 days of past data and a four-day forecast of the velocity field). In addition, finer analysis with respect to the crossing events in the singular value curves will be done in the future. The study [23] also highlighted the structural flow pattern changes during the polar vortex split, but the corresponding behavior of the singular values was not analyzed in much detail. Moreover, it would be challenging to study the full three-dimensional system, which would give much more insight into the Lagrangian aspects of the vortex split.

Conclusions
A set-oriented dynamical system approach was used to study early-warning signals of splitting of coherent sets in case studies. This study is mainly implemented around singular value decompositions of stochastic transition matrices, which are generated from finitetime evolutions of time-dependent systems. Thus, dominant singular vectors and their sign patterns yield numerical approximations of slowly mixing regions of Lagrangian particle dynamics, while corresponding singular values exhibit the potential to measure their shape changes. On this basis, trends in the singular values and sign patterns of singular vectors are mutually used as observables to anticipate sudden changes that are imminent in the non-autonomous dynamics of vortex-like patterns. This technique was used to identify early indicators of the sudden split of the Antarctic polar vortex in September 2002 from recorded satellite velocity data, and this real-world critical transition may be classified as a rare event. Though it was demonstrated on two-dimensional incompressible example systems, the probabilistic approach can, in principle, be applied to a wider class of systems, including three-dimensional flows.
The proposed methodology has some limitations: Several parameters, both in the underlying systems and the numerical approximations, influence the magnitude of the spectral changes that can be observed, ruling out general statements about significance and reliable forecasting of critical events. However, early-warning signals (such as autocorrelation and variance) that are frequently used as statistical indicators of critical transitions have similar issues [1][2][3]. The approach presented here analyzes mass transport in physical space, and is thus restricted to low-dimensional systems. It also does not take into account other relevant system information, such as pressure or temperature, which are highly relevant in the observation of atmospheric flows and numerical weather prediction. We note that despite our low-dimensional context in this paper, transfer operator techniques have been applied successfully in high-dimensional systems, including molecular dynamics (see, e.g., [40]), but in such contexts, a reduction to effective lower-dimensional state spaces is necessary. In addition, in [18], early-warning indicators of transitions between atmospheric flow regimes were studied using almost-invariant sets, which are time-independent versions of coherent sets. These almost-invariant sets were computed as numerical approximations of particular meta-stable regimes (using a reduced transfer operator, as proposed in [41]) detected in an atmospheric barotropic model via principal component analysis techniques (see also [42]).
The aim of our work was to apply dynamical systems and ergodic theoretical concepts to study changes in Lagrangian flow patterns in example systems, which might inspire further theoretical and practical advances. From a mathematical perspective, future studies may address a combination of geometric methods and the above probabilistic approach to find lower bounds of the dominant singular values in order to abstractly control the changes in the patterns. In view of the recently introduced data-based approaches for coherent structure identification [43][44][45][46][47] that make use of spectral graph clustering, it will be interesting to study the extent to which similar indicators of flow pattern changes, as discussed in this paper, are applicable. These may be the first steps towards using Lagrangian probabilistic methodologies in real operational settings. As, in principle, forecasts can be included in the computation of flow patterns, it would be challenging to study how uncertainties in the velocity forecasts influence the resulting flow patterns.
The changes in the global dynamics are also visible in the structure of the transition matrices L in Figure A1a-c and Figure A2a-c, where the nonzero entries are plotted in blue color. Here, the horizontal axes correspond to the reachable states at time t + 1, while the vertical axes correspond to the initial state at time t. We note that for t 0, the final state space Y t+1 is only a small subset of the initial state space X t as a consequence of the contraction induced by the presence of the global attractor {0} of the original system, but it grows as the system approaches and passes through the bifurcation point at t = 0.
The contraction of all initial points towards the trivial attractor and the subsequent set-valued perturbation results in an overall mixing, so that there are no proper coherent sets as long as t < 0. In particular, for a t 0 that is sufficiently small, there exists a maximal neighborhood R t of 0 that is contained in the set-valued image Φ t+τ t (x) of every point x ∈ [−2, 2]: R t = x∈[−2,2] Φ t+τ t (x). This is due to the strong contraction towards 0 and the (relatively) large choice of ε. Thus, a lot of mixing takes place in this region. As t approaches 0, the attraction towards the origin gets weaker and weaker, and as a consequence, R t shrinks and becomes empty close to the bifurcation point (for our choice of ε = 0.4, this happens at t ≈ −1.67).
This is also visible in the leading left and right singular vectors for times t < 0 (see Figure A1d-i). While the left singular vectors are unaffected by approaching the bifurcation point, the right singular vectors clearly display R t as a shrinking zero plateau in a neighborhood of the origin for both the second and the third singular vectors. For final times t + 1 ≥ 0, such a region R t is empty, which is also strikingly visible in the corresponding singular vectors in Figure A2d-i.   [1,2] in (c,f,i).