Koopman early warning signals for bifurcation and rate-induced tipping

Introduction

Critical transitions, or tipping points, are abrupt shifts in system dynamics that can occur when a stable regime loses resilience or can no longer track changing external conditions. Such transitions have been studied across complex systems, including in ecology, biology, and climate science [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. Because these transitions can be severe, irreversible, and strongly amplified by stochastic fluctuations, considerable effort has been devoted to the development of early warning signals (EWS) capable of anticipating their onset [13, 14, 15, 16].

Most classical EWS are rooted in bifurcation-induced tipping (B-tipping), where a slowly varying parameter drives the system towards loss of local stability. Near such transitions, perturbations recover increasingly slowly, a phenomenon known as critical slowing down (CSD) [17]. In stochastic systems, CSD typically manifests as increasing lag-1 autocorrelation (AR1) and variance [18], making these quantities among some of the most widely used indicators. Spectral theory based on combining Koopman framework [19] with response theory also provides a mathematical explanation for this behaviour [20, 15]. A distinct mechanism arises in rate-induced tipping (R-tipping) [21, 22, 23], where transitions occur not because stability is lost, but because external conditions change too rapidly for the system to track the evolving attractor. In this case, the instantaneous dynamics may remain locally stable throughout the transition pathway, so that classical indicators based on CSD may fail to provide reliable warning signals [24].

In realistic stochastic and externally forced systems, different tipping mechanisms may coexist and can be difficult to distinguish across multiple interacting timescales [25, 26, 27]. This is especially relevant in climate dynamics, where internal variability, nonstationary forcing, and multiscale feedbacks jointly shape the onset of critical transitions. Existing EWS are typically tailored to either stability loss (B-tipping) or tracking failure (R-tipping) separately, raising the need for frameworks that remain informative across both regimes.

Here we address this challenge using Koopman operator theory, which represents finite-dimensional nonlinear dynamics through the linear evolution of infinite-dimensional observables [28, 29, 19, 30]. Previous works have connected B-tipping and CSD to Koopman spectral structure, including spectral gap closure (Figure 1a), the amplification of stochastic fluctuations (Figure 1b) [15, 31, 32, 33, 20, 34], and the resulting deterioration of Koopman representation quantified through residual Koopman mode decomposition (ResKMD) near critical transition [34]. In this work, we extend this idea to the rate-induced setting using a control-aware Koopman formulation [35], addressing both tipping mechanisms under a unified framework.

We evaluate our framework across a hierarchy of systems, ranging from prototypical models to reduced-order and fully coupled simulations of the Atlantic Meridional Overturning Circulation (AMOC), a widely discussed climate tipping element. We further show that learned embeddings through deep learning (DL), such as through multi-layer perceptron (MLP)-based autoencoders, outperform prescribed observables, such as random Fourier (RFF) [36] and time-delay features (TDF [37]), in terms of better separation of tipping and non-tipping trajectories, especially in high-dimensional settings. Beyond early warning detection, our framework also provides interpretable spectral diagnostics, including decay timescales, spectral gap structure, and changes in memory preceding critical transition. Together, these results demonstrate how operator-theoretic methods can unify prediction and mechanistic interpretation of different tipping phenomena in stochastic systems.

Generalized Koopman framework for early warning of tipping points

The central quantity in our framework is the residual obtained through ResKMD [38, 39, 40]. This residual measures how accurately the observed dynamics can be represented using a finite set of Koopman modes. Far from tipping, the dynamics are often well described by a compact low-dimensional spectral representation dominated by a small number of Koopman modes. Near tipping, we argue that this representation progressively deteriorates: in the bifurcation setting because slow modes accumulate near the stability boundary, and in the rate-induced setting because the system can no longer track the moving attracting branch under changing external conditions. The Koopman residual therefore acts as a dynamical measure of loss of spectral representability (see Figure 1).

Bifurcation-induced setting. In the autonomous setting relevant to B-tipping, we consider a stochastic dynamical system with state $\omega_t \in\mathbb{R}^{N}$ and observables $\psi(\omega_t) \in\mathbb{C}^{M}$ evolved through the first-moment Koopman operator $\mathcal{K}^{(1)}:\mathbb{C}^{M}\to\mathbb{C}^{M}$. We approximate the dynamics using a retained Koopman subspace $\mathcal{S}_{m}^{b}$. In the bifurcation setting, the Koopman residual separates naturally into two contributions: a truncation component measuring imperfect finite-dimensional Koopman representation, and a stochastic component measuring variability in the observable evolution

Here $\mu$ denotes a reference measure on state space and $\Phi_{\varepsilon}$ the one-step stochastic map. Near B-tipping, both contributions are theoretically shown to increase as slow Koopman modes approach the stability boundary and stochastic fluctuations become amplified [34]. In this sense, ResKMD generalizes classical CSD indicators such as AR1 and variance within a spectral operator-theoretic framework [41] (see discussion in Supplementary Information S1.1, S1.2).

Rate-induced setting. The R-tipping setting requires a different interpretation. Here the key issue is not necessarily loss of local stability, but failure of the system to track a moving attracting branch under time-dependent forcing. We therefore augment the Koopman representation with a control variables describing the forcing history [38, 39]. Let $u_t \in\mathbb{R}^{U}$ denote the control variable evolving according to $u_{t+1}=g(u_t)$, let $\Psi(\omega,u) \in\mathbb{C}^{M}$ be an observable dictionary on the joint state-control space, and let $\mathcal{K}_{c}^{(1)}:\mathbb{C}^{M}\to\mathbb{C}^{M}$ denote the corresponding first-moment Koopman operator. The resulting control-aware residual retains the same structure as in the bifurcation setting, separating truncation and stochastic contribution, but now in the joint state-control space

Proposition 1: Control-aware Koopman residual decomposition

The proof is provided in Supplementary Information S2.2.

In order to make the rate-induced mechanism explicit, we define the tracking error relative to the moving attracting branch (see fig:1c)

which measures the deviation of the evolving system state from the instantaneous attracting branch associated with the control variable $u_t$. In the rate-induced setting, tipping corresponds to growth of this tracking error despite the dynamics remaining locally stable. Linearizing around $\omega^{\star}(u_t)$ yields (derivation in Supplementary Information S2.1)

The first term describes local restoration toward the attracting branch, while the second term $d_t$ denotes the branch drift induced by the changing control variable, and $\epsilon_t$ collects stochastic and higher-order contributions. In this picture, R-tipping occurs when branch drift is not sufficiently compensated by local restoration.

We assume that the tracking error lies in, or is well approximated by, the span of the chosen observables. Concretely, there exists a bounded linear map $C:\mathbb{C}^{M}\to\mathbb{R}^{N}$ relating observables to the tracking error $\bar{\omega}_t$. Motivated by the tracking error dynamics above, our main theoretical result connects increased control-aware residuals to unresolved branch drift.

Theorem 1: Residual lower bound from unresolved branch drift Assume that the local restoring dynamics are captured by the retained Koopman subspace, then the control-aware residual satisfies:

where $Q_{m d}$ and $Q_{m\epsilon}$ denote the components of the branch drift and high-order terms unresolved by the retained Koopman representation. Proof is provided in Supplementary Information S2.4.

Theorem 1 shows that the control-aware residual increases when branch drift, $d_t$, contains components unresolved by the retained Koopman representation. Thus, unlike B-tipping, where residual growth reflects spectral slowing near instability, increased residuals in R-tipping reflect loss of tracking under time-dependent forcing. Nonetheless, unresolved branch drift is not the only possible source of increased residuals. Additional unresolved dynamics through e.g., local restoring dynamics not fully captured by the retained Koopman subspace, and stochastic fluctuations may also contribute. Theorem 1 instead identifies unresolved branch drift as a minimal mechanism through which rate-induced tracking failure contributes to residual growth.

In summary, both of these residual quantities in (1) and (2) will be used as our Koopman-based EWS in B- and R-tipping, respectively. See Methods and Supplementary Information for complete derivation and formalism.

Koopman EWS skills on prototypical systems

We first benchmark our proposed Koopman early warning indicators on four prototypical stochastic systems spanning both B- and R-tipping before applying the framework to AMOC simulation data. Aside from traditional indicators such as AR1 and variance, we also plot the Koopman spectral gap, defined as $\gamma:= -\operatorname{Re}(\lambda_{L,1})$, where $\lambda_{L,1}$ is the leading non-zero Koopman generator eigenvalue, i.e., the non-zero eigenvalue with real part closest to zero. Near CSD, $\gamma$ is expected to decrease toward zero [31]. For B-tipping, we consider May’s harvesting model and the Rosenzweig–MacArthur consumer–resource model, which exhibit fold and Hopf bifurcations, respectively [42]. For R-tipping, we consider a forced saddle-node normal form and a Bautin oscillator, where transitions arise from failure to track a moving attracting branch under time-dependent forcing [24]. Additional system details, parameter choices, and implementation procedures are provided in Methods and Supplementary Information.

Figure 2 (see also Supplementary Section S3 for additional results) shows representative results for the Rosenzweig–MacArthur model (B-tipping) and the Bautin oscillator (R-tipping). In the bifurcation setting, Koopman residuals increase more clearly for tipping trajectories than for non-tipping trajectories. This is consistent with the theoretical mechanism: as the system approaches the bifurcation, slow Koopman modes move toward the stability boundary, finite-dimensional spectral representations deteriorate, and stochastic fluctuations become amplified along slowly decaying directions. Learned observables obtained from MLP embeddings further improve the separation between tipping and non-tipping trajectories relative to prescribed dictionaries such as random Fourier features (RFF), consistent with the quantitative AUROC scores shown in Figure 8.

The rate-induced setting reveals a qualitatively different picture. Classical CSD indicators show inconsistent separation between tipping and non-tipping trajectories. For instance, AR1 and variance increase in both tipping and non-tipping cases, while the spectral gap approaches zero near the critical transition in both groups. The control-aware Koopman residual, in contrast, increases substantially more for tipping trajectories, consistent with the unresolved branch drift mechanism described earlier. This supports the interpretation that, in R-tipping, the dominant precursor is loss of tracking rather than CSD, and is also consistent with the quantitative AUROC scores shown in Figure 8.

Koopman EWS skills on empirical observations

We further evaluate our Koopman EWS on a diverse collection of real-world observations spanning human and natural systems. These datasets include cyanobacteria recovery dynamics in a microcosm experiment [43], phonation onset of vocal cord [44], ATP collapse under hypoxia in plant cells [45], paleoclimate transitions associated with the end of greenhouse Earth [46], and electricity blackout [47]. Together, these systems probe a broad range of temporal scales, noise characteristics, and underlying transition mechanisms beyond idealized dynamical models. Additional dataset descriptions and details are provided in the Methods section.

Figure 3 (see also Figure S4 for additional results) shows that the Koopman-based indicators remain informative across these heterogeneous empirical settings. In systems associated with classical CSD, such as the microcosm and paleoclimate examples, the Koopman residual exhibits increasing trends consistent with growing dynamical memory and spectral slowing prior to transition. In more complex or noisy settings, including the blackout and phonation examples, the residual-based indicators remain sensitive to changes in the underlying dynamical structure even when conventional indicators appear less consistent. Overall, these results suggest that the Koopman framework can extract informative dynamical signatures from empirical observations across a wide range of domains and observational conditions.

Tipping of the Atlantic Meridional Overturning Circulation

The Atlantic Meridional Overturning Circulation (AMOC) is a major component of the climate system that transports heat and salinity through the Atlantic Ocean. Its weakening or collapse has long been discussed as a potential climate tipping phenomenon because it can reorganize large-scale ocean circulation and strongly affect regional and global climate [4, 10]. We therefore use AMOC as a physically relevant testbed for our Koopman framework in both the reduced box-model setting and coupled global circulation model (GCM).

For the reduced setting (see Figure 4 for sample realizations), we adopt the two-dimensional box-model from [23, 48], obtained from a five-box ocean model by treating Southern Ocean and bottom water salinities as slow variables and diagnosing Indo-Pacific salinity from salt conservation. The resulting state variables are the North Atlantic and tropical Atlantic salinities $(S_N, S_T)$, from which the AMOC strength $Q$ is diagnosed (see Methods (48)) and supplied as the scalar observable to all early warning estimators. The salinity dynamics are piecewise-defined according to the sign of $Q$ (see Methods (49)–(52)) and driven by surface freshwater hosing $H(r,t)$. Departing from the original deterministic formulation, we add white noise with amplitude $\sigma_N = \sigma_T = 1.0$ to both salinity equations. We configure the model for both tipping classes: in the B-tipping setup, $H$ is ramped linearly past the deterministic fold at $H^\ast= 0.4236$ [48], with non-tipping runs capped just below threshold; in the R-tipping setup, tipping and non-tipping trajectories follow an identical hyperbolic-secant forcing profile (see Methods (53)) differing only in the forcing rate $r$.

As illustrated in Figure 5, the results from AMOC box-model mirror the patterns seen in the prototypical systems. In the B-tipping setting, AR1, variance and ResKMD (RFF and MLP) all exhibit increasing trend prior to transition, consistent with CSD phenomena. In the R-tipping setting, by contrast, classical indicators fail to provide useful early warning signal as the non-tipping trajectories tend to be falsely characterized as tipping (e.g., increasing variance, diminishing spectral gap). The control-aware ResKMD, however, provides a better tipping signal as the indicator abruptly rises near the tipping threshold at a much faster rate than their CSD-based counterpart. This is also evident in our quantitative AUROC measures in Figure 8 where control-aware ResKMD outperforms across the board.

We next consider the 10-member NASA GISS-E2.1-G ensemble studied in [49] under the SSP2-4.5 scenario. The ensemble exhibits branching behavior in AMOC strength at $48^{\circ}\mathrm{N}$: eight members remain on a relatively strong overturning branch, whereas two evolve toward a markedly weaker state.

fig:6a illustrates this separation using depth profiles of the Atlantic overturning streamfunction, with non-tipping realizations retaining a stronger overturning cell and tipping realizations showing pronounced weakening, consistent with the AMOC strength trajectories in fig:6b. As a physically motivated control variable, we use annual Denmark Strait sea-ice mass flux, which [49] linked to the weakening branch through enhanced freshwater release (fig:6c). For this study, we treat the weak-branch members as tipping trajectories and the remainder as non-tipping trajectories ($26^\circ$N results are reported in Supplementary Information S5.1). Before supplying this control variable to the control-aware ResKMD, we smooth it using a 10-year running mean. Sensitivity to alternative control variables, such as maximum mixed layer depth (MLD) across different ocean basins, is reported in Supplementary Information S5.2.

fig:7a shows that the control-aware Koopman indicators provide the clearest separation between tipping and non-tipping AMOC trajectories in the coupled NASA GISS ModelE ensemble at $48^\circ$N. Classical CSD-based metrics, especially AR1, remain highly variable across ensemble members and do not show a consistent increase prior to trajectory divergence. Variance exhibits a smoother trend, but still provides weaker separation. By contrast, control-aware ResKMD, particularly with MLP observables, increases more strongly and earlier for the tipping members, with a marked rise around year 2150, while remaining comparatively small for the non-tipping members. In the MLP case, the tipping trajectories eventually decrease and converge after the transition, consistent with relaxation toward a different attracting state. These results suggest that, in the coupled setting, the relevant precursor is not fully captured by indicators based on the AMOC time series alone. Instead, it is more clearly expressed as a mismatch between AMOC evolution and the supplied control variable in the augmented Koopman space. This interpretation is consistent with the idea proposed in [50] that the strong trajectory separation observed in [49] may be associated with a close encounter with an edge state, which acts as a dynamical saddle and reduces predictability; see also [51].

Beyond early warning detection, our Koopman framework provides access to interpretable spectral diagnostics. Here we focus on the $48^\circ$N results using the TDF observables setting (see Figure S5 for similar results in $28^\circ$N). We first examine the distribution of the discrete-time Koopman eigenvalues $\lambda_{K,j} \in\mathbb{C}$ (see fig:7b). The tipping ensemble exhibits a broader and more outward radial distribution, with eigenvalues lying closer to the unit circle. Since the decay of discrete-time Koopman modes is governed by $\lvert\lambda_{K,j}\rvert$, this indicates slower decay, weaker damping, and thus reduced restoring stability of the stronger AMOC branch in the tipping ensemble.

To obtain physically interpretable decay rates and oscillation frequencies, we map to the continuous-time generator spectrum through $\lambda_{K,j} = e^{\lambda_{L,j}\Delta t} = e^{(\sigma_j+i\omega_j)\Delta t}$, where $\sigma_j = \operatorname{Re}(\lambda_{L,j})$ is the exponential decay rate, $\omega_j = \operatorname{Im}(\lambda_{L,j})$ is the angular frequency, and $\Delta t$ has annual resolution. From this we compute the e-folding time $T_{e,j} = 1/|\sigma_j|$. We report the dominant $T_e$ associated with the slowest-decaying non-constant mode. As illustrated in Figure 7, the non-tipping ensemble tends to exhibit larger $T_e$, suggesting a more persistent slow mode associated with an eventual recovery toward the stronger AMOC branch rather than collapse.

Finally, we examine spectral gap behavior by repeating the EDMD fit in a rolling window of length 20% of the available record and retaining, within each window, the three modes with smallest $|\operatorname{Re}(\lambda_{L,j})|$. Figure 7 shows that the tipping ensemble tends to exhibit slow modes that move more quickly toward the stability boundary $\operatorname{Re}(\lambda_{L,j}) = 0$ and remain closer to zero than in the non-tipping ensemble. This indicates spectral densification and a reduced spectral gap in the tipping case, consistent with a weakening of the isolated point-spectrum description. It is also consistent with the larger residuals, which indicate that the dynamics are less well captured by a small number of isolated Koopman modes prior to AMOC weakening.

Discussion

In this work, we develop a Koopman-based framework for EWS in stochastic systems undergoing either B-tipping or R-tipping. In the B-tipping setting, we built on recent results connecting ResKMD to CSD and showed how the residual provides a principled indicator that combines truncation error and stochastic fluctuation within a single spectral quantity. With a control setup we generalise this framework to R-tipping scenarios, yielding a control-aware ResKMD that captures loss of tracking relative to a moving attracting branch. Unlike the B-tipping case, where residual growth is tied to eigenvalues approaching the unit circle, the R-tipping mechanism is instead governed by unresolved drift associated with the time-varying control variables.

Across a hierarchy of examples, from low-dimensional prototypical systems to reduced box-model and fully coupled AMOC simulations, we found that the proposed Koopman-based indicators recover the expected warning behavior in B-tipping problems and provide useful discrimination in R-tipping regimes where classical EWS such as AR1 and variance can be weak or inconsistent. In the AMOC applications, the framework also yielded useful spectral diagnostics beyond detection itself, including $T_e$, spectral gaps, and the distribution of slow Koopman modes. These quantities provide additional physical interpretation of resilience loss and timescale separation, and highlight a central advantage of the operator-theoretic viewpoint: the same EDMD/ResDMD machinery used to construct the EWS also gives access to interpretable spectral structure. In real-world application where only a single realization is available, our ResKMD indicator can provide robust uncertainty estimate by performing multiple runs with different observable configuration.

Some considerations point to a number of promising active research directions. The practical performance of ResKMD and its control-aware extension depends on the choice of observables, truncation level, and the supplied control or driver variable. As seen in the AMOC experiments, some control choices are substantially more informative than others, and identifying the most dynamically relevant forcing channel is itself a nontrivial modeling problem. Furthermore, in the AMOC R-tipping box-model experiments, although the control-aware ResKMD indicator exhibits a substantially faster increase in the tipping trajectories, this increase alone may be difficult to interpret in isolation when non-tipping reference trajectories are unavailable. This highlights an important limitation for practical deployment in real-world systems when R-tipping dominates and counterfactual non-tipping realizations are typically inaccessible. Future work should therefore investigate more robust calibration strategies, uncertainty quantification, and methods for identifying absolute rather than purely comparative signatures of rate-induced loss of tracking.

While the B-tipping case admits a clearer theoretical divergence result, the R-tipping case is more problem dependent: increased residuals arise through unresolved branch drift rather than through a universal spectral route (see our discussion in Remark 1). Our results also rely on finite-sample EDMD- and ResDMD-type approximations, whose statistical error, dictionary dependence, and robustness in high-dimensional settings deserve further study. On the theory side, sharper guarantees for control-aware ResKMD under partial observability, model misspecification, and nonstationary noise would further strengthen the framework’s foundations, as would connecting the present framework more explicitly to Koopman generators, continuous spectrum, and ResDMD with control in more general settings. Understanding how residual growth relates to tracking failure in nonautonomous stochastic systems is a particularly interesting open question. Just like B-tipping can be seen as a manifestation of the closure of the spectral gap and of the divergence of response operators in the autonomous setting [20, 15], it is tempting to investigate whether recent advances on linear response for nonautonomous systems [52, 53] can shed light on R-tipping and be related to the results presented here. On the application side, an important next step is to test these EWS more broadly across other climate tipping elements and ESM simulations, where multiple coupled control variables, spatial heterogeneity, and partial observations [54, 55] make tipping detection especially challenging.

Overall, to our knowledge our results present the first framework of an unified EWS analysis across B-tipping and R-tipping in stochastic systems using operator-theoretic formalism. Rather than treating these two tipping mechanisms with separate heuristic indicators, ResKMD offers a common language in which both loss of stability and loss of tracking can be studied through spectral approximation error. This study provides a significant step toward EWS methods that remain flexible and physically interpretable in the stochastic, high-dimensional, and strongly driven systems that arise in AMOC and other complex applications.

Methods

Stochastic dynamics

Consider a discrete stochastic dynamical system with state $\omega_t \in\Omega\subseteq\mathbb{R}^N$:

Here, $\Phi:\Omega\to\Omega$ is the deterministic flow map and $\epsilon_t$ denotes stochastic forcing. Unless otherwise stated, we assume that $\epsilon_t$ has zero conditional mean and finite second-moment, so that the relevant expectations and covariances are well defined. In the synthetic experiments, we use independent Gaussian noise with covariance $\Sigma_{\epsilon}$. Expectations below are taken with respect to the probability law induced by the noise sequence:

Bifurcation-induced tipping

To define bifurcation-induced tipping (B-tipping), consider a parameterized family of deterministic dynamical systems $\omega_{t+1} = \Phi_{\beta}(\omega_t)$, where $\beta\in\mathbb{R}$ is a slowly varying parameter. Suppose that, for each $\beta< \beta^{\star}$, the system admits a stable fixed point, $\omega^{\star}(\beta)$, satisfying $\Phi_{\beta}(\omega^{\star}(\beta)) = \omega^{\star}(\beta)$. The set of equilibria $\omega^{\star}(\beta)$ defines an attracting equilibrium branch parameterized by $\beta$. To characterize local stability along this branch, we linearize the dynamics around $\omega^\star(\beta)$. The corresponding Jacobian is

Let $\{\lambda_{\mathcal{J},i}(\beta)\}_{i=1}^{N}$ denote the eigenvalues of $\mathcal{J}(\beta)$ ordered so that $|\lambda_{\mathcal{J},1}(\beta)| \ge\cdots\ge|\lambda_{\mathcal{J},N}(\beta)|$. We define the spectral radius of the Jacobian as $\rho(\mathcal{J}(\beta)) := \max_i |\lambda_{\mathcal{J},i}(\beta)|$. Local asymptotic stability requires $\rho(\mathcal{J}(\beta)) < 1$, so the spectrum lies strictly inside the unit circle. As the system evolves along the equilibrium branch toward the critical parameter value $\beta\to\beta^\star$, the leading spectral modes approach the unit circle, signaling weakening restoring stability and the onset of CSD.

For a fixed parameter value $\beta$, linearizing the perturbation $\bar{\omega}_t := \omega_t - \omega^\star(\beta)$ gives

where $\epsilon_t$ collects stochastic forcing and neglected higher-order terms. Let $\ell_{\mathcal{J},1}(\beta)$ denote the left eigenvector associated with the dominant eigenvalue $\lambda_{\mathcal{J},1}(\beta)$, and define

with $\eta_t := \ell_{\mathcal{J},1}(\beta)^\top\epsilon_t$ denoting the projected stochastic innovation. This yields the familiar spectral interpretation of CSD. Assuming stationarity, $\operatorname{Var}(z_{t+1}) = \operatorname{Var}(z_t)$, and uncorrelated innovations, $\operatorname{Cov}(z_t,\eta_t) = 0$, we obtain AR1 (derivation in Supplementary Information S1.1)

and variance (derivation in Supplementary Information S1.2)

These relations are exact for one-dimensional linear systems and provide a useful leading mode approximation when the local dynamics are dominated by a single slow mode. If several modes decay on comparable timescales, AR1 and variance should instead be interpreted as effective scalar summaries rather than exact spectral diagnostics. Thus, as $\rho(\mathcal{J}(\beta)) \to1$, perturbations decay more slowly, and both AR1 and variance increase.

Rate-induced tipping

Rate-induced tipping (R-tipping) refers to the loss of tracking of a moving attracting branch caused by sufficiently rapid variation in a time-dependent control variable, even when the corresponding frozen dynamics remain locally stable [21, 22]. We consider the stochastic control-dependent system

where $u_t \in\mathbb{R}^{U}$ is a time-dependent control variable. For each fixed value of $u$, we associate the frozen system $\omega_{t+1} = \Phi(\omega_t,u)$, obtained by holding the control constant. Assume that, for each $u$ in the regime of interest, this frozen system admits an equilibrium $\omega^\star(u)$ satisfying $\Phi(\omega^\star(u),u) = \omega^\star(u)$, so that $u \mapsto\omega^\star(u)$ defines a frozen equilibrium branch. We further assume that this branch is locally stable for each fixed $u$, with

Thus, if the control variable is frozen at any single value, nearby trajectories would relax locally toward the corresponding equilibrium $\omega^\star(u)$. The key issue in the non-autonomous setting is whether the evolving trajectory $\omega_t$ remains close to the moving equilibrium branch $\omega^\star(u_t)$ as the control variable changes in time. To quantify this, we define the tracking error, $\bar{\omega}_t := \omega_t - \omega^\star(u_t)$. Linearizing around the moving branch $\omega^\star(u_t)$ yields (derivation in Supplementary Information S2.1)

The term $d_t := \omega^*(u_t)-\omega^*(u_{t+1})$ is the control-induced drift of the frozen equilibrium branch, that is, the displacement of the instantaneous attracting state caused by the change in the control variable from $u_t$ to $u_{t+1}$. (16) therefore shows that R-tipping arises through competition between two effects: local restoration toward the current frozen equilibrium, represented by $\mathcal{J}(u_t)\bar{\omega}_t$, and control-driven drift, represented by $d_t$. In particular, tipping can occur even when $\rho(\mathcal{J}(u)) < 1$ for every fixed $u$, provided the control changes rapidly enough that the trajectory can no longer track the moving branch.

Stochastic Koopman operators

We now define the Koopman-theoretic framework used throughout the paper. Koopman theory represents nonlinear state evolution through linear operators acting on observables $\psi: \Omega\to\mathbb{C}$ [28, 29]. We let $\mu$ be a finite positive measure on the state space $\Omega$. Consider the space $L^2(\Omega,\mu)$ equipped with the standard inner product

The first-moment stochastic Koopman operator $\mathcal{K}^{(1)}$ acts on observables of $L^2(\Omega,\mu)$ by expectation

As a linear operator on $L^2(\Omega,\mu)$, $\mathcal{K}^{(1)}$ can be studied through its spectrum, whose structure depends on the underlying dynamics. When Koopman eigenfunctions exist, the corresponding eigenpairs $(\lambda,\phi)$ satisfy

To connect Koopman spectral structure with stochastic variability and classical CSD indicators such as variance, we additionally define a second-moment Koopman operator

The corresponding one-step variance is

In the data driven setting, we use observable dictionaries built from random Fourier features (RFF [36]), time-delay features (TDF [37]), and deep embeddings, such as using multi-layer perceptron (MLP)-based autoencoders [56].

Residual Koopman framework for B-tipping

We now formalize the residual Koopman framework used for the autonomous, bifurcation-induced setting. Throughout this subsection, the sub/superscript $b$ denotes the bifurcation setting. Let $\psi= [\psi_1,\ldots,\psi_M]^{\top}$ be a vector-valued observable dictionary of dimension $M$. When non-constant Koopman eigenfunctions are available, let $\{(\lambda_n^b,\phi_n^b)\}_{n\geq1}$ denote the retained eigenvalue–eigenfunction pairs, ordered by decreasing magnitude. We exclude the constant invariant eigenfunction with eigenvalue one, since it does not represent decay or early warning behavior. For a truncation level $m \in\mathbb{N}^+$, we define:

Following the standard Koopman mode decomposition [28], the truncated Koopman representation of the vector observable $\psi$ on $\mathcal{S}_m^b$ is:

where $v_n^b \in\mathbb{C}^M$ is the Koopman mode associated with eigenfunction $\phi_n^b$. We then define the residual Koopman mode decomposition (ResKMD) by

This residual separates into a deterministic approximation component and a stochastic fluctuation component. Specifically

Here $\operatorname{Cov}$ denotes the covariance of the next-step observable $\psi(\Phi_\varepsilon(\omega))$ induced by stochasticity at fixed current state $\omega$. The first term in (25) measures the error incurred by representing the expected observable evolution using only the retained Koopman subspace $S_m^b$, while the second term measures the stochastic spread of the lifted one-step evolution.

To interpret the truncation term, consider the simplifying case in which the relevant Koopman eigenfunctions form an orthonormal basis for the observable subspace of interest. Then the observable can be expanded as $\psi=\sum_{n\geq1}v_n^b\phi_n^b$, and the truncation error becomes

With the orthonormal-basis and eigenvalue-ordering assumptions above, the truncation error is bounded below by the first omitted Koopman mode and above by a worst-case tail contribution controlled by the largest omitted eigenvalue [34]

Moreover, for each non-constant Koopman eigenfunction, the stochastic one-step evolution can be written as

Here $\eta_{n,t}$ is the one-step stochastic innovation associated with the Koopman eigenfunction $\phi_n^b$. Under the same stationarity and innovation-independence assumptions used in Equations 12 and 13, this gives

From (27), the contribution of the first omitted mode can increase as its eigenvalue approaches the unit circle. Similarly, (29) further shows that stochastic fluctuations become larger when a non-constant Koopman eigenvalue approaches the unit circle. Intuitively, this happens because perturbations decay more slowly, so noise persists for longer. Thus, near B-tipping, ResKMD can increase through two complementary mechanisms: poorer finite-dimensional spectral approximation and stronger stochastic amplification. This mechanism is formalized for the bifurcation setting by Theorem 2 of [34], which shows, under its stated assumptions, that

Control-aware residuals for R-tipping

We now apply the Koopman residual framework to the augmented state-control system introduced above. Throughout this subsection, the sub/superscript $c$ denotes the control setting. Using $u_{t+1}=g(u_t)\in\mathbb{R}^U$, the corresponding one-step map on the joint state-control space $\Omega\times\mathbb{R}^U$ has the form

This rewrites the control-dependent dynamics as an autonomous stochastic system on the joint state-control space. We assume that the control update $g:\mathbb{R}^U\to\mathbb{R}^U$ is deterministic, so that stochasticity enters only through the state transition $\Phi_\varepsilon(\omega,u)$. This keeps the one-step uncertainty collected in a single stochastic term for notational simplicity.

Let $\mu$ be a finite reference measure on $\Omega\times\mathbb{R}^U$, such as the empirical measure induced by an observed state-control trajectory. We define state-only, joint state-control, and control-only observable blocks by

The concatenated control-aware dictionary is

The three blocks need not have the same dimension. Their specific construction is implementation-dependent. In our numerical experiments, we construct the control variable directly i.e., using the identity observable, and the joint observable using simple products between the state observables and the control variable, following [35].

The control-aware first-moment Koopman operator acts on observables of the augmented state-control space as

Thus, $\mathcal{K}_{c}^{(1)}$ is simply the stochastic Koopman operator associated with the augmented map in (31). The control-aware formulation therefore applies the same Koopman residual idea as before, but now to the pair $(\omega,u)$ rather than to $\omega$ alone.

Let $\{(\lambda_{n}^{c},\phi_{n}^{c})\}_{n\geq1}$ denote the retained non-constant eigenvalue-eigenfunction pairs of $\mathcal{K}_{c}^{(1)}$, ordered by decreasing magnitude. We assume that such point-spectrum components exist. For a truncation level $m\in\mathbb{N}^{+}$, we define:

Following the standard Koopman mode decomposition, the truncated control-aware representation is

where $\mathbf{v}_{n}^{c}\in\mathbb{C}^{M}$ is the Koopman mode associated with $\phi_{n}^{c}$. We define the control-aware residual by

Because this is the usual ResKMD residual applied to the augmented state-control system, the same decomposition into approximation error and stochastic fluctuation carries over. This gives rise to Proposition 1, and detailed proof is given in Supplementary Information S2.2.

Remark 1: No universal residual divergence in the R-tipping setting. Unlike in the B-tipping case, one should not expect a universal residual-divergence result for R-tipping, because rate-induced transitions need not involve loss of local stability or Koopman eigenvalues approaching the unit circle. If the augmented state-control dynamics are represented exactly in the chosen observable space, the control-aware residual can remain zero even when the trajectory loses track of the moving equilibrium branch; a counterexample is given in Supplementary Information S2.3. Thus, the residual should be interpreted as detecting unresolved tracking dynamics in a finite-dimensional data-driven representation, rather than tracking failure in a perfectly resolved model.

We now connect the control-aware residual to the lower bound mechanism in Theorem 1. Recall the tracking error, $\bar{\omega}(\omega,u):=\omega-\omega^{\star}(u)$ and branch drift, $d(\omega,u):=\omega^{\star}(u)-\omega^{\star}(g(u))$ as defined earlier, but now adapted using the control dynamics $g$.

Assumption 1: Tracking error observability. There exists a bounded linear map $\mathcal{C}:\mathbb{C}^{M}\to\mathbb{R}^{N}$, such that $\bar{\omega}(\omega,u)=\mathcal{C}\mathbf{\Psi}(\omega,u)$.

Assumption 1 says that the chosen observables contain enough information to reconstruct the tracking error. In practice this may hold only approximately, but it makes explicit the link between residuals in the augmented observable space and tracking errors in state space.

Using the tracking error dynamics from (16), we write the first-moment evolution of $\bar{\omega}$ as

(37) separates the expected tracking error evolution into local restoration, branch drift, and remaining higher-order terms from linearization. The lower bound follows by projecting this decomposition onto the part not represented by the retained Koopman subspace. Let $\mathcal{P}_{m}$ denote projection onto $\mathcal{S}_{c}^{m}$, and let $\mathcal{Q}_{m}:=I-\mathcal{P}_{m}$ denote the unresolved projection. Under the tracking error observability Assumption 1, and assuming that the local restoration term is captured by the retained Koopman subspace, the control-aware residual satisfies

If $\mathcal{Q}_{m}\epsilon=0$, this further reduces to

The proof is provided in Supplementary Information S2.4. Theorem 1 shows that unresolved branch drift provides a lower-bound contribution to the control-aware residual. This does not mean that all residual growth is caused by branch drift: unresolved state dynamics, stochastic fluctuations, imperfect observability, and finite-dimensional approximation error may also contribute. Rather, the theorem identifies one mechanism by which rate-induced loss of tracking can appear as increased Koopman residual.

Extended dynamic mode decomposition

We estimate finite-dimensional autonomous and control-aware Koopman operators from trajectory data using extended dynamic mode decomposition (EDMD) [57]. Given a trajectory $\{\omega_t\}_{t=1}^{T+1}$, in the autonomous setting, we form one-step snapshot pairs $(\omega_t,\omega_{t+1})$ and lift each observation into a higher-dimensional feature space using a dictionary of observable $\psi\in\mathbb{C}^{M}$ constructing two data matrices:

The Koopman approximation is then given as $K\approx\Psi_{Y}\Psi_{X}^{\dagger}$, where $\dagger$ is the Moore-Penrose pseudoinverse [58]. In the control-aware setting, we proceed analogously on the augmented trajectory $\{(\omega_t,u_t)\}_{t=1}^{T+1}$ using dictionary $\Psi(\omega,u)$ (see Equation 32):

from which we estimate $K_c\approx\Psi_{Y_c}\Psi_{X_c}^{\dagger}$ jointly captures the state-control dynamics.

We consider both prescribed (random Fourier features (RFF) [36] and time-delay features (TDF) [37]) and deep learning-based (latent embeddings from multi-layer perceptron (MLP) autoencoders [59, 56]) for our choice of observables. Unless otherwise stated, we retain Koopman modes explaining 90% of the variability. By default, we also use RFF with 500-dimensional and TDF with 20-step delay, unless otherwise specified. For MLP observables, we use the default hyperparameters of learning rate $10^{-3}$, batch size 32, and training epochs of 1000. Training is performed once over the first window and the frozen encoder is applied over the subsequent sliding windows. In summary, this yields a fully data-driven procedure that does not require explicit knowledge of the governing equations while retaining the spectral interpretation of the Koopman framework.

Quantitative metrics

In addition to the qualitative rolling-window analysis of the indicators, we quantitatively benchmark the proposed Koopman EWS using Kendall’s $\tau$ [60] to measure monotonic trends prior to tipping, and receiver operating characteristic (ROC) analysis [61, 62] together with the corresponding area under the ROC curve (AUROC) [63, 64] to evaluate discrimination between tipping and non-tipping trajectories.

We consider a binary classification setting where each sample is labeled as either tipping ($y=1$) or non-tipping ($y=0$). Let $\hat{s}:\mathcal{D}\to\mathbb{R}$ denote a scoring function that assigns a confidence score to each observed trajectory, where higher values indicate a higher likelihood of a critical transition. A prediction $\hat{y}=1$ is made whenever $\hat{s}>\tau$ for a threshold $\tau$.

For each threshold $\tau$, we define the true positive rate (TPR) and false positive rate (FPR) as

The receiver operating characteristic (ROC) curve [61, 62] is obtained by plotting $\operatorname{TPR}(\tau)$ against $\operatorname{FPR}(\tau)$ as the threshold $\tau$ varies over all possible values. The ROC curve characterizes the trade-off between correctly detecting tipping events and incorrectly classifying non-tipping trajectories.

The Area Under the ROC Curve (AUROC) [63, 64] is defined as

which corresponds to the probability that a randomly chosen tipping trajectory is assigned a higher score than a randomly chosen non-tipping trajectory. An AUROC of 0.5 indicates random performance, while 1.0 corresponds to perfect discrimination.

Prototypical tipping systems

This section provides the data-generation protocols and model specifications for the prototypical systems and AMOC experiments presented in the main text.

For bifurcation-induced tipping, we consider the following systems:

We first model the May’s single-species harvesting model for a biomass variable $x$, corresponding to the one-dimensional case $\omega= x \in\mathbb{R}_{+}$. Its continuous-time dynamics are

where $r$ is the intrinsic growth rate, $k$ is the carrying capacity, $s$ sets the saturation scale of harvesting, $h$ is the harvesting rate, and $\xi_x(t)$ denotes Gaussian white noise. The harvesting term counteracts logistic growth, and as $h$ increases the stable positive-biomass equilibrium eventually disappears in a fold bifurcation. In our experiments, we set $r = 1$, $k = 1$, $s = 0.1$, and $\sigma= 0.01$, for which the deterministic system loses stability at the critical harvesting rate $h^\star= 0.26$. Tipping runs are generated by ramping $h$ linearly from 0.15 to 0.27, while non-tipping runs use the same initial ramp but capped at $h = 0.25$, keeping the system below the fold.

We next consider the Rosenzweig-MacArthur consumer-resource model, for which the state is two-dimensional, $\omega= (x, y)^{\mathsf{T}} \in\mathbb{R}_{+}^{2}$, with $x$ the resource population and $y$ the consumer population. The continuous-time dynamics are

where $r$ is the intrinsic growth rate of the resource, $k$ is its carrying capacity, $a(t)$ is the consumer attack rate, $e$ is the conversion efficiency, $h$ is the handling time, $m$ is the consumer mortality rate, and $\xi_x, \xi_y$ are independent Gaussian white noise processes. We use $r = 4$, $k = 1.7$, $e = 0.5$, $h = 0.15$, and $m = 2$, for which the coexistence equilibrium destabilizes at the Hopf threshold $a^\star= 15.69$. Tipping runs are produced by ramping $a$ from 12 to 17, thereby crossing $a^\star$, whereas non-tipping runs use the same ramp but capped at $a = 14.38$, so the system remains on the stable side of the Hopf bifurcation.

For both systems, we generate 50 tipping and 50 non-tipping stochastic trajectories with noise amplitude $\sigma= 0.01$. Integration uses the Euler–Maruyama scheme with internal step size $\Delta t = 0.01$. EWS are computed on pre-transition time series of length 500 using sliding windows of size 50% and stride 1. The data-generation protocol follows [42], except that for the non-tipping cases we use the same parameter schedule as the tipping cases before capping the bifurcating parameter slightly below the critical threshold, thereby making the discrimination task more challenging.

For rate-induced tipping, we consider the following systems:

We first consider a one-dimensional saddle-node normal form with state $\omega= x \in\mathbb{R}$ and control variable $u(t) = \lambda(t)$. The continuous-time dynamics are

where $\lambda_{\max}$ sets the forcing amplitude and $\epsilon$ controls the forcing rate. For each fixed $\lambda$, the frozen system has the usual saddle-node equilibrium structure, but the trajectory can still fail to track the attracting branch if $\lambda(t)$ increases too rapidly. In our experiments, we use $\lambda_{\max}=3$ and compare a tipping case with $\epsilon=1.25$ against a non-tipping case with the same forcing profile at half the rate, $\epsilon=0.625$. We set $\sigma_x=0.008$. We use the scalar observable $x(t)$ for EWS estimation, and identify tipping when the trajectory crosses the threshold $x \ge0$.

We next consider a stochastic Bautin system with complex state $z_t=x_t+\mathrm{i}y_t$ and forcing variable $u(t)=\Lambda(t)$. The continuous-time dynamics are

where $a$ is the linear growth parameter, $\omega$ is the angular frequency, $b$ is the cubic coefficient, $\sigma_z$ sets the noise magnitude, and $r$ controls the forcing rate. In this example the forcing translates the underlying Bautin normal form through state space, and sufficiently rapid variation in $\Lambda(t)$ can trigger rate-induced escape even when the dynamics remain locally stable. We use the parameter values $a=0.1$, $\omega=3$, $b=1$, $\sigma_z=0.2$, and $\Lambda_{\max}=8$. Tipping runs use $r=0.10$, whereas non-tipping runs use the same forcing shape with half the rate, $r=0.05$. For EWS estimation we use the scalar radius observable $\rho(t)=\sqrt{x_t^2+y_t^2}$, and identify tipping when $\rho(t)\ge10$.

For the reduced AMOC experiments, we use the box-model analyzed in [23, 48]. This system is obtained from the original five-box model by treating the Southern Ocean (S) and bottom water (B) salinities as slow or fixed and diagnosing the Indo-Pacific (IP) salinity from salt conservation, yielding a two-dimensional dynamical system with state $\omega=(S_N,S_T)^\top\in\mathbb{R}^2$, where $S_N$ and $S_T$ denote the North Atlantic and tropical Atlantic salinity variables, respectively. The AMOC strength is diagnosed as

and the observable supplied to the EWS estimators is the scalar time series $Q(t)$. The model is driven by freshwater hosing through the surface freshwater fluxes. In the reduced formulation, the salinity dynamics are piecewise-defined according to the sign of $Q$.

For $Q\ge0$,

For $Q<0$,

Here, $F_N$ and $F_T$ denote the freshwater-flux terms induced by the hosing protocol $H(r,t)$, with parameters given in [23]. Different from earlier formulation, we transform the deterministic box-model as stochastic system where noise amplitudes are set to $\sigma_N=\sigma_T=1.0$.

For the B-tipping configuration, we use a linear schedule for the bifurcation parameter $H$: tipping runs ramp $H$ from $0.0$ to $0.4736$, whereas non-tipping runs follow the same ramp but are capped at $H=0.3736$. The deterministic fold occurs at $H^*=0.4236$ [48]. For the R-tipping configuration, we use

so that tipping and non-tipping trajectories differ only through the forcing rate $r$. We use $r=0.017$ for tipping and $r=0.005$ for non-tipping, with $t_{\mathrm{crit}}=500$.

Empirical observations

microcosm We analyze a cyanobacteria population in chemostats subjected to dilution events under gradually increasing light levels [43]. Population density is inferred from the light attenuation coefficient, computed from continuous measurements of outgoing light intensity. The full dataset comprises 7,784 observations over 28.86 days at a sampling interval of 5 minutes. The time series is divided into six segments separated by dilution events, which are external perturbations and not part of the intrinsic population dynamics. To ensure consistency, we restrict analysis to recovery phases between dilution events. Specifically, for each segment, we use the final 250 time points (approximately one day) prior to the next dilution event, yielding a total of 1,500 data points across all segments.

voice We analyze experimental time series of phonation onset, corresponding to the emergence of vocal fold oscillations, which can be interpreted as a Hopf bifurcation [44, 65]. The data are obtained from a physical replica of human vocal folds (EPI model) with a multilayer body–cover structure that reproduces realistic tissue mechanics. Oscillations are induced by gradually increasing airflow, while subglottal pressure is recorded using a pressure transducer. The available time series (from [41]) consists of 7,501 data points over 0.375 s with a sampling interval of $5 \times10^{-5}$ s. Tipping corresponds to the transition from a stable, non-oscillatory state to sustained oscillations as the flow rate crosses a critical threshold. We use the measured pressure signal as the observable for tipping detection.

cellular_atp We analyze experimental time series of cytosolic ATP (i.e., the readily available energy pool) dynamics in living plant cells, measured via a genetically encoded FRET sensor sensitive to $\mathrm{MgATP}^{2-}$ [45]. The data capture the response of leaf tissue under gradually increasing hypoxia, induced by sealing samples in a dark environment where respiration depletes available oxygen. ATP levels remain relatively stable during oxygen decline and then undergo a sudden collapse, indicating a critical transition in cellular energy state. The dataset consists of 271 observations over 3 minutes with a sampling interval of 0.011 minutes. We use the fluorescence ratio signal as the observable, with tipping identified as the abrupt drop in ATP concentration.

greenhouse_earth We analyze a paleoclimate time series of calcium carbonate ($\mathrm{CaCO}_{3}$) associated with the end of greenhouse Earth, marking the transition from an ice-free state to the formation of polar ice caps [46]. The data exhibit increasing autocorrelation prior to the climate shift, consistent with critical slowing down. The dataset consists of 462 observations spanning 5.9 million years with a sampling interval of 0.013 million years. Tipping is identified as the transition to a glaciated climate state.

blackout_frequency We analyze bus voltage frequency data preceding the Western Interconnect blackout of August 1996, measured within the Bonneville Power Administration network [47]. The time series captures system dynamics leading up to grid separation and exhibits critical fluctuations prior to the blackout, consistent with EWS of instability. The dataset consists of 11,301 observations over 565 seconds with a sampling interval of 0.05 s. Tipping is identified as the transition to system-wide failure, i.e., blackout.

Declarations

Acknowledgements

JN, HF, PG acknowledge funding, computing, and storage resources from the NSF Science and Technology Center (STC) Learning the Earth with Artificial Intelligence and Physics (LEAP) (Award #2019625). The material is based upon work supported by NASA under award No 80NSSC25K0062. VL acknowledges the partial support provided by the Horizon Europe Projects Past2Future (Grant No. 101184070) and ClimTIP (Grant No. 100018693), by the ARIA SCOP-PR01-P003 - Advancing Tipping Point Early Warning AdvanTip project, and by the European Space Agency Project PREDICT (Contract 4000146344/24/I-LR). Further support was provided by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research and Office of Biological and Environmental Research, Scientific Discovery through Advanced Computing (SciDAC) program under Award Number 9233218CNA000001. PK acknowledges the support from U.S. Department of Energy (DOE), Office of Science, Office of Biological and Environmental Research, Regional and Global Model Analysis program area as part of the HiLAT-RASM project.

Competing interests

All authors declare no financial or non-financial competing interests.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Data Availability

The toy datasets have been generated with custom code (see the Code Availability statement) and simulation data to reproduce the AMOC results is freely available from the reference in the text.

Materials availability

Not applicable.

Code Availability

The code for this study will be available at https://github.com/juannat7/kews.

Author contributions

JN and CR conceived of the study and performed the analysis under the guidance of PG. DD, PK, and VL assisted in the methodological formulation. HF and AR provided domain science expertise to frame and interpret the results. The first draft of the manuscript was written by JN and CR, and all authors commented on the manuscript.

References

Spectral view of classical indicators

The classical CSD interpretation of AR1 and variance derived below is inherently equilibrium-centered: one linearizes the dynamics around a stable fixed point and studies the dominant eigenvalue of the resulting Jacobian. This provides the mechanism by which slowing recovery implies increasing autocorrelation and variance near bifurcation. However, this reasoning assumes that the local stability structure is well approximated by a nearly fixed equilibrium over the analysis window. For methods to estimate a local, time-varying Jacobian from data around the observed state itself, see discussion in [41]. This class of methods extend Jacobian-based warning analysis beyond the strict requirement of first locating the equilibrium. The underlying bifurcation intuition remains the same, but the estimation is performed locally along the trajectory rather than exclusively at an equilibrium.

Spectral view of AR1

where the second line applies stationarity to simplify the denominator, and the third uses uncorrelatedness of innovation to eliminate the cross term. As $\rho(\mathcal{J}(\beta)) \to1$, perturbations decay increasingly slowly, indicating greater persistence. For a real dominant eigenvalue, this implies $|\operatorname{Corr}(z_t,z_{t+1})| \to1$. In particular, $\operatorname{Corr}(z_t,z_{t+1}) \to1$ when the leading eigenvalue approaches $+1$, as in a saddle-node (fold) bifurcation, whereas $\operatorname{Corr}(z_t,z_{t+1}) \to-1$ when the leading eigenvalue approaches $-1$, as in a flip (period-doubling) bifurcation. A complete derivation of the link between closure of the spectral gap and occurrence of $|\operatorname{Corr}(z_t,z_{t+1})| \to1$ is presented in [15, 31, 32].

Spectral view of variance

where the second line uses the uncorrelatedness of innovation to drop the cross term, and the third applies stationarity. As $\rho(\mathcal{J}(\beta)) \to1$, the variance term increases and diverges in the linear approximation. A more general derivation of this property is presented in [66].

Derivation and Proof

Derivation of the R-tipping’s tracking error

Given $\bar{\omega}_t := \omega_t - \omega^\star(u_t)$, we have

where the fourth line adds and subtracts $\Phi(\omega^\star(u_t),u_t)$ and the fifth line uses the equilibrium identity $\Phi(\omega^\star(u_t),u_t) = \omega^\star(u_t)$. Linearizing $\Phi(\cdot,u_t)$ around $\omega^\star(u_t)$ gives

and hence

Proof of Proposition 1

Proof. First we define

For fixed $(\omega,u)$, add and subtract $\mathbb{E}[Y \mid\omega,u]$

Expanding the squared Euclidean norm gives

Conditioning on $(\omega,u)$ and taking expectations yields

The cross term vanishes because

Therefore,

Since $\mathbb{E}[Y \mid\omega,u] = [K^{(1)}_{c}\Psi](\omega,u)$ by the definition of the control-aware first-moment Koopman operator in (33), the first term in (64) is exactly the squared truncation error at $(\omega,u)$.

The second term equals the trace of the conditional covariance,

Integrating (64) over the joint state-control space $\Omega\times\mathbb{R}^{U}$ with respect to measure $\mu$ gives

which proves Proposition 1. □

Throughout, expectations and covariances in the control-aware setting are understood with respect to the one-step stochastic transition from the current pair $(\omega,u)$, and are integrated over $(\omega,u)$ against the reference measure $\mu$.

Counterexample for Remark 1

Consider the linear forced system

For fixed $u$, the frozen equilibrium $\omega^{\star}(u)$ satisfies

Next, we define the tracking error

Then

Thus the frozen dynamics remain stable as $0<a<1$, but sufficiently large $r$ prevents tracking.

We now consider the following observable dictionary

Its one-step lifted evolution is

with

Hence the augmented dynamics close exactly in the chosen observable space. Because the system is deterministic, both the truncation error and stochastic fluctuation vanish, and therefore

even though the trajectory fails to track the moving equilibrium branch.

Proof of Theorem 1

Proof. Since each component of the local restoration, $\mathcal{J}(u)\bar{\omega}(\omega,u)$ lies in $\mathcal{S}^{c}_{m}$, its unresolved projection onto $Q_{m}$ gives zero vector. Applying $Q_{m}$ to (37) therefore gives

Hence, by the reverse triangle inequality,

From $\bar{\omega}=C\Psi$ and linearity of $C$ and $\mathcal{K}^{(1)}_{c}$, we have

Since $\mathcal{K}^{(1)}_{c}\Psi\in(\mathcal{S}^{c}_{m})^{M}$, each component of $C(\mathcal{K}^{(1)}_{c}\Psi)$ also lies in $\mathcal{S}^{c}_{m}$. By the best-approximation property of $\mathcal{P}_{m}$,

Combining the two inequalities gives

By Proposition 1, the residual is bounded below by the truncation error because the stochastic fluctuation term is nonnegative. This proves (38). If each component of $\epsilon$ also lies in $\mathcal{S}^{c}_{m}$, then $Q_{m}\epsilon=0$, then (39) follows immediately as a special case.

Prototypical Experiments

Below represents the tipping and non-tipping trajectories across prototypical B-tipping and R-tipping systems.

Bifurcation-induced tipping

Rate-induced tipping

Empirical Experiments

We provide more empirical results showcasing our Koopman-based EWS ability to correctly identify critical transitions.

Additional Experiments

Coupled AMOC results at $26^{\circ}\mathrm{N}$

The main text focuses on $48^{\circ}\mathrm{N}$. For completeness, Figure S5 reports the corresponding coupled AMOC results at $26^{\circ}\mathrm{N}$ where similar results and conclusions can still be drawn. The control-aware ResKMD again provides better discrimination and earlier detection of early warning signals than AR1 and variance.

AMOC control variable ablation

To assess sensitivity of the control-aware Koopman residual to the choice of forcing channel, we repeat the coupled AMOC analysis using several alternative control variables. In place of the Denmark Strait sea-ice flux used in the main text, we consider the maximum mixed-layer depth (MLD) in the Labrador Sea, Irminger Sea, Scotland-Iceland Basin, and Greenland-Iceland-Norway Seas. For each forcing variable, we apply the same preprocessing and compute control-aware ResKMD with TDF observables for AMOC strength at both 26$^\circ$N and 48$^\circ$N.

Figure S6 compares the forcing trajectories and the corresponding control-aware residuals across the tipping and non-tipping ensembles. The results show that the performance of control-aware ResKMD is sensitive to forcing choice: some forcings produce clearer separation between tipping and non-tipping members than others. This supports the interpretation that control-aware Koopman indicators are most informative when the supplied control captures the dominant mechanism associated with AMOC weakening.