# EFFECTIVE REDUCED MODELS FROM DELAY DIFFERENTIAL EQUATIONS: BIFURCATIONS, TIPPING SOLUTION PATHS, AND ENSO VARIABILITY

MICKAËL D. CHEKROUN AND HONGHU LIU

**ABSTRACT.** Conceptual delay models have played a key role in the analysis and understanding of El Niño-Southern Oscillation (ENSO) variability. Based on such delay models, we propose in this work a novel scenario for the fabric of ENSO variability resulting from the subtle interplay between stochastic disturbances and nonlinear invariant sets emerging from bifurcations of the unperturbed dynamics.

To identify these invariant sets we adopt an approach combining Galerkin-Koornwinder (GK) approximations of delay differential equations and center-unstable manifold reduction techniques. In that respect, GK approximation formulas are reviewed and synthesized, as well as analytic approximation formulas of center-unstable manifolds. The reduced systems derived thereof enable us to conduct a thorough analysis of the bifurcations arising in a standard delay model of ENSO. We identify thereby a saddle-node bifurcation of periodic orbits co-existing with a subcritical Hopf bifurcation, and a homoclinic bifurcation for this model. We show furthermore that the computation of unstable periodic orbits (UPOs) unfolding through these bifurcations is considerably simplified from the reduced systems.

These dynamical insights enable us in turn to design a stochastic model whose solutions—as the delay parameter drifts slowly through its critical values—produce a wealth of temporal patterns resembling ENSO events and exhibiting also decadal variability. Our analysis dissects the origin of this variability and shows how it is tied to certain transition paths between invariant sets of the unperturbed dynamics (for ENSO’s interannual variability) or simply due to the presence of UPOs close to the homoclinic orbit (for decadal variability). In short, this study points out the role of solution paths evolving through tipping “points” beyond equilibria, as possible mechanisms organizing the variability of certain climate phenomena.

**Keywords.** Center Manifold Reduction | Galerkin-Koornwinder Approximations | Stochastic Modeling | Transition Paths

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. Galerkin-Koornwinder (GK) approximations of DDEs</td><td>3</td></tr><tr><td>3. Center-Unstable Manifold Approximations from GK systems</td><td>6</td></tr><tr><td>4. Bifurcations Analysis of a Delay El Niño-Southern Oscillation (ENSO) Model</td><td>11</td></tr><tr><td>5. Tipping Solution Paths and ENSO variability</td><td>21</td></tr><tr><td>6. Concluding Remarks</td><td>26</td></tr><tr><td>Acknowledgements</td><td>27</td></tr><tr><td>Appendix A. Determination of P<sub>N</sub> and Q<sub>N</sub> in Eq. (4.10)</td><td>27</td></tr><tr><td>References</td><td>28</td></tr></table>

## 1. INTRODUCTION

Since its inception in the 80s and 90s as conceptual models to study the El Niño-Southern Oscillation (ENSO) phenomenon [SS88, BH89, TSCJ94, NBH<sup>+</sup>98], delay models have attracted a growing attention in climate modeling; see [GT00, GZT08, KF11, RCC<sup>+</sup>14, GCS15, KTF17, KKP17, CGN18, BCL<sup>+</sup>17, FQS<sup>+</sup>19].

*Date:* January 2, 2024.
However, only recently the bifurcation analysis of these conceptual delay models arising in climate has been undertaken with modern tools dedicated to delay models; see e.g. [KS14, KKP15, KKP16, KKD19, CKL20].

In this work, we pursue this endeavor by providing a detailed analysis of the bifurcations arising in the ENSO delay model of Suarez and Schopf [SS88]. Partial bifurcation results have been reported about this model. These results were often treated in a subsidiary fashion within a work of more general scope [FQS<sup>+19</sup>, AR23] relegating details about certain calculations and dynamical insights to the reader. Here, we provide a thorough analysis of these bifurcations and associated calculations to help the reader build up intuitions in view of applications, in particular regarding the stochastic modeling of ENSO as discussed also in this study. Our approach employs the Galerkin method based on the Koornwinder basis functions [Koo84] as introduced in [CGLW16] to derive rigorous low-dimensional ordinary differential equation (ODE) approximations (GK systems for short), to which a reduction to the center-unstable manifold is applied, following [CKL20]. Here, both the derivation of GK systems and the center-unstable manifold calculations are made self-contained and more transparent than in [CGLW16, CKL20] to better serve the “practitioner” interested in applications.

For the Suarez and Schopf delay model, the resulting two-dimensional reduced ODE system demonstrates remarkable skills in predicting the local and global bifurcations including a subcritical Hopf bifurcation, a homoclinic bifurcation, and a saddle-node bifurcation of periodic orbits (SNO bifurcation) when the delay parameter $\tau$ is varied. Following [CKL20, Theorem III.1], the subcritical Hopf bifurcation is proved by analyzing the sign of the Lyapunov coefficient of the reduced system. Theorem III.1 of [CKL20] simplifies the calculation of this coefficient to that of basic operations consisting of solving linear algebraic systems with triangular matrices, computing the GK spectrum at the critical delay parameter, and forming the involved inner products; see (4.16) below. The detection of homoclinic and SNO bifurcations benefits from the stable and accurate approximation skills of the Unstable Periodic Orbits (UPOs) to the Suarez and Schopf model, by simple backward integration of the reduced system. The computation of UPOs is known to be challenging for high-dimensional problems and to require a careful formulation and implementation of Newton methods and the like [Gri08, Gri13]. Here, our reduced system allows thus for bypassing these difficulties.

There is a relatively vast body of literature dealing with Galerkin methods exploiting other basis functions to approximate delay model by ODEs, including step functions [BB78, KS78], splines [BK79, BRI84], linear and sine functions [Vya12, WC05], and Legendre polynomials [Kap86, IT86]. Techniques based on pseudospectral collocation methods have also emerged as an alternative route to approximate delay differential equations (DDEs) by ODEs [BDG<sup>+16</sup>] and have shown relevance in computing characteristic roots [BMV05], Lyapunov exponents [BVV14], and bifurcation analysis [ABL<sup>+22</sup>] using ODE bifurcation software packages such as AUTO [Aut] and MatCont [DGK<sup>+08</sup>, Mat]. These methods, although showing great promises, do not have yet developed into toolboxes such as DDE-BIFTOOL [ELR02, SEL<sup>+14</sup>] and KNUT [Knu] which provide currently the most evolved software packages to analyze bifurcations arising in DDEs with discrete delays. Nevertheless, ODE approximation methods offer a framework that is more flexible to analyze bifurcations, covering a broader range of applications as encompassing delay models with (possibly nonlinear) functional of distributed delays [CGLW16], or renewal equations [BDG<sup>+16</sup>].

Yet, the effective derivation of reduced systems from the ODE approximations produced by such methods, seems to have been poorly exploited. As shown in the case of GK systems, to dispose of such reduced systems has many potential attributes that await to be tested in applications involving DDEs. Steps in this direction include e.g. the effective computation of low-dimensional solutions to (nearly) optimal controls in feedback form, enabling to avoid solving the infinite-dimensional Hamilton-Jacobi-Bellman (HJB) equation associated with the optimal control problem of the original DDE [CKL18].

Reduced models derived from GK systems have also shown recently their relevance to improve the realism of DDE's solutions via stochastic modeling, as illustrated in cloud physics [CKLL22]. Conceptual delay models were proposed in [KF11, KTF17] to envision open cellular cells formed by marine stratocumulus clouds as an oscillatory predator-prey mechanism of clouds (prey) scavenged by rain (predator). Such conceptual models produce oscillations that, albeit grounded on physical intuitions, are too regular (periodic) to bear the
comparison with e.g. satellite observations or high-end model simulations. In [CKLL22], a significant step was taken to bridge the gap between the latter and these conceptual models. To do so, the reduced systems derived in [CKLL22] from GK approximations of the DDE model of [KF11] allowed for identifying nonlinear structures providing the level set of the oscillation's phase function, known as isochrons [Guc75, ACN16]. Current algorithms to compute isochrons [MRMM14, DDMG16] suffer the “curse of dimensionality,” making the direct computation of isochrons out of reach, for DDEs. In [CKLL22], the knowledge of the isochrons in the reduced state space was shown to be sufficient to design stochastic parameterizations interacting with the genuine DDE's isochrons, leading in turn to stochastic chaotic solutions with new virtues including enhanced time-variability mimicking that of clouds' oscillations.

In this study, benefiting from the high-accuracy approximation skills achieved by our 2D reduced GK system, we reach a detailed understanding of the DDE phase portrait when the delay parameter $\tau$ is varied; see Fig. 4 below. This understanding allows us to dissect the response of the Suarez and Schopf model when subject to stochastic disturbances while the delay parameter drifts slowly through the aforementioned bifurcations. As shown in Sec. 5.1, we obtain indeed tipping solution path (TSP) whose certain time episodes resemble ENSO-like temporal patterns that are explained in terms of transition paths between invariant sets of the DDE's unperturbed dynamics, such as locally stable/unstable steady states and UPOs. By conducting in Sec. 5.2 a spectral analysis of the TSP's frequency content, we show furthermore that the TSP's irregular behavior on longer timescales is dominated by decadal variability such as documented in recent ENSO studies [DCP+21]. The dynamical insights gained in Sec. 4 allow us to trace back the origin of this decadal variability in our stochastic model which results from the presence of UPOs located close to the homoclinic orbit.

The remainder of this paper is organized as follows. We first recall the key analytic and algebraic elements for the formal construction of the Galerkin-Koornwinder (GK) approximations of general scalar DDEs in Sec. 2. We survey in Sec. 3 analytic approximation formulas of center-unstable manifolds—including leading-order (Theorem 3.1) and higher-order formulas (Sec. 3.2)—that we present for the reduction of GK systems experiencing a loss of stability at a critical value of the delay parameter $\tau$. These formulas are then applied to the Suarez and Schopf model in Sec. 4 to derive an effective 2D reduced GK system with $\tau$-dependent coefficients (Eq. (4.14)). Sections 4.2 and 4.3 present the reduced system skills in predicting for the delay model a subcritical Hopf bifurcation, and an SNO bifurcation co-existing with an homoclinic bifurcation, respectively. Section 4.4 provides numerical evidences of the highly-accurate approximation skills achieved by the reduced systems, in particular regarding the computation of UPOs and limit cycles. Section 4.5 gives then model error estimates complementing these numerical results. Finally, in Sec. 5, we take advantage of the dynamical insights gained in Sec. 4 to design a stochastic model whose solutions exhibit ENSO-like patterns and decadal variability. Some final remarks and potential future directions are then outlined in Sec. 6.

## 2. GALERKIN-KOORNWINDER (GK) APPROXIMATIONS OF DDEs

We consider nonlinear scalar DDEs of the form

$$ (2.1) \qquad \frac{dx(t)}{dt} = ax(t) + bx(t - \tau) + c \int_{t-\tau}^{t} x(s) \, ds + F(x(t), x(t-\tau), \int_{t-\tau}^{t} x(s) \, ds), $$

where $a, b$ and $c$ are real numbers, $\tau > 0$ is the delay parameter, and $F$ is a nonlinear function. We restrict ourselves to the scalar case to simplify the presentation, but the approach extends to systems of nonlinear DDEs involving possibly several delays as detailed in [CGLW16] and illustrated in [CKL20] for the cloud-rain delay model of [KF11].

### 2.1. Koornwinder polynomials.
First, let us recall that Koornwinder polynomials $K_n$ are obtained from Legendre polynomials $L_n$, for any nonnegative integer $n$, according to the relation

$$ (2.2) \qquad K_n(s) = -(1+s) \frac{d}{ds} L_n(s) + (n^2 + n + 1) L_n(s), \quad s \in [-1, 1], $$
see [CGLW16, Eq. (3.3)].

Koornwinder polynomials are known to form an orthogonal set for the following weighted inner product on [-1, 1] with a point-mass, $\mu(dx) = \frac{1}{2}dx + \delta_1$, where $\delta_1$ denotes the Dirac point-mass at the right endpoint $x = 1$; see [Koo84]. In other words, the following orthogonality property holds:

$$ (2.3) \qquad \begin{aligned} & \int_{-1}^{1} K_n(s) K_p(s) \, d\mu(s) = \frac{1}{2} \int_{-1}^{1} K_n(s) K_p(s) \, ds + K_n(1) K_p(1) \\ & = 0, \quad \text{if } p \neq n. \end{aligned} $$

It is also worthwhile noting that Koornwinder polynomials augmented by the right endpoint values as follows

$$ (2.4) \qquad \mathcal{K}_n = (K_n, K_n(1)), $$

form an orthogonal basis of the product space

$$ \mathcal{E} = L^2([-1, 1); \mathbb{R}) \times \mathbb{R}, $$

endowed with the inner product:

$$ (2.5) \qquad \left\langle (f, p), (g, q) \right\rangle_{\mathcal{E}} = \frac{1}{2} \int_{-1}^{1} f(s)g(s) \, ds + pq, \quad \text{for all } (f, p), (g, q) \in \mathcal{E}. $$

The norm induced by this inner product is denoted by $\| \cdot \|_{\mathcal{E}}$. The basis function $\mathcal{K}_n$ has then its $\| \cdot \|_{\mathcal{E}}$-norm given by [CGLW16, Prop. 3.1]:

$$ (2.6) \qquad \|\mathcal{K}_n\|_{\mathcal{E}} = \sqrt{\frac{(n^2 + 1)((n + 1)^2 + 1)}{2n + 1}}. $$

This is a useful property for calculating the GK approximations, as it will come apparent below.

Finally, since the original Koornwinder basis is given on the interval [-1, 1] and the state space of a DDE such as Eq. (2.1) is made of functions defined on (-τ, 0), we will work mainly with the following rescaled version of Koornwinder polynomials. The rescaled Koornwinder polynomials $K_n^\tau$ are defined as follows

$$ (2.7) \qquad \begin{aligned} K_n^\tau: [-\tau, 0] &\to \mathbb{R}, \\ \theta &\mapsto K_n \left(1 + \frac{2\theta}{\tau}\right). \end{aligned} $$

They form orthogonal polynomials on the interval [-τ, 0] for the $L^2$-inner product on (-τ, 0) with a Dirac point-mass at the right endpoint 0.

The following family of rescaled Koornwinder polynomials augmented with a right endpoint value,

$$ (2.8) \qquad \mathcal{K}_n^\tau = (K_n^\tau, K_n^\tau(0)), $$

forms then an orthogonal basis for the Hilbert space $\mathcal{H} = L^2([-\tau, 0); \mathbb{R}) \times \mathbb{R}$ endowed with the inner product $\langle \cdot, \cdot \rangle_{\mathcal{H}}$ given by (2.5) in which the integral is taken on (-τ, 0) and weighted by $1/\tau$. Note that since $K_n(1) = 1$ [CGLW16, Prop. 3.1], we have also

$$ (2.9) \qquad K_n^\tau(0) = 1. $$

Finally, observe that $\|\mathcal{K}_n^\tau\|_{\mathcal{H}} = \|\mathcal{K}_n\|_{\mathcal{E}}$.

## 2.2. Space-time representation and GK approximations.
It is well-known that a DDE such as Eq. (2.1) is an infinite-dimensional dynamical system [Hal88] in which a history segment has to be specified over the interval [-τ, 0) to apprehend properly the existence and computation problems of its solutions [HVL93, CZ95, BZ13].

Thus, given a function of time, $x$, solving Eq. (2.1) one distinguishes between the history segment $\{x(t + \theta) : \theta \in [-\tau, 0)\}$ and the current state, $x(t)$. Denoting by $u(t, \theta)$ the history segment, one can then rewrite the DDE (2.1) as the transport equation

$$ (2.10) \qquad \partial_t u = \partial_\theta u, \quad -\tau \le \theta < 0, $$
subject to the *nonlocal and nonlinear boundary condition*

$$ (2.11) \qquad \partial_{\theta} u|_{\theta=0} = au(t, 0) + bu(t, -\tau) + c \int_{t-\tau}^{t} u(s, 0) \, ds \\ \qquad \qquad + F(u(t, 0), u(t, -\tau), \int_{t-\tau}^{t} u(s, 0) \, ds), \quad \text{for } t \ge 0. $$

This reformulation is often called the space-time representation in the literature; see Fig. 1 for an illustration.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Panel A shows a solution $u(t, \theta)$ to a transport problem of the form (2.10)–(2.11) obtained as a reformulation of the DDE $\dot{x} = ax(t - \tau) - bx(t - \tau)^3$ from [CGLW16]. Panel B shows the time series $u(t, 0)$ at the right endpoint $\theta = 0$. It coincides with the time series $x(t)$ obtained by solving the original DDE. The solution is shown for $a = 0.5$, $b = 20$ and $\tau = 0.4$.</figcaption></figure>

Using the rescaled Koornwinder polynomials (2.7), one can then seek for approximations of $u(t, \theta)$ solving Eqns. (2.10)-(2.11) under the form

$$ (2.12) \qquad u_N(t, \theta) = \sum_{j=0}^{N-1} y_j(t) K_j^{\tau}(\theta). $$

Given the property (2.9), the current state $x(t) = u(t, 0)$ is approximated by

$$ (2.13) \qquad x_N(t) = \sum_{j=0}^{N-1} y_j(t). $$

Note that in (2.12) the Koornwinder polynomials are indexed according to their degree $j$.

The question arises then of determining the equations that the coefficients $y_j(t)$ in (2.12) must satisfy in order to have that $u_N(t, \theta)$ provides a rigorous approximation that converges to $u(t, \theta)$ as $N \to \infty$. This nontrivial problem has been solved in [CGLW16]. The next section recalls the structure of these equations forming what we call a GK system.
2.3. **The formulas of GK approximations.** The analysis of [CGLW16] shows that $y_j(t)$ solving the $N$-dimensional system of ODEs (2.14) provides a rigorous method for approximating the solutions to a broad class of DDEs; see [CGLW16, Sec. 5]. Such systems are given by:

$$ (2.14) \qquad \begin{aligned} \frac{dy_j}{dt} = \frac{1}{\|\mathcal{K}_j\|_{\mathcal{E}}^2} \sum_{n=0}^{N-1} & \left( a + bK_n(-1) + c\tau(2\delta_{n,0} - 1) \right. \\ & \left. \qquad + \frac{2}{\tau} \sum_{k=0}^{n-1} a_{n,k} (\delta_{j,k} \|\mathcal{K}_j\|_{\mathcal{E}}^2 - 1) \right) y_n(t) \\ & + \frac{1}{\|\mathcal{K}_j\|_{\mathcal{E}}^2} F \left( \sum_{n=0}^{N-1} y_n(t), \sum_{n=0}^{N-1} y_n(t)K_n(-1), \tau y_0(t) - \tau \sum_{n=1}^{N-1} y_n(t) \right), \\ 0 \le j \le N-1. \end{aligned} $$

Here the Kronecker symbol $\delta_{j,k}$ has been used, and the coefficients $a_{n,k}$ are obtained by solving a triangular linear system in which the right-hand side has explicit coefficients depending on $n$ [CGLW16, Prop. 5.1]; see Appendix A. The values of the $K_n(-1)$'s are known and recalled in (A.1) below, and we recall that the formula for the $\|\mathcal{K}_n\|_{\mathcal{E}}$ is given by (2.6).

We rewrite now the GK system (2.14) into the following compact form:

$$ (2.15) \qquad \frac{dy_j}{dt} = A_N(\tau)y_j + G_2(y_j, \tau), $$

where $y_j = (y_0, \dots, y_{N-1})^T$, and in which the terms $A_N(\tau)y_j$ and $G_2(y_j, \tau)$ group the linear and nonlinear terms in Eq. (2.14), respectively. The entries of the $N \times N$ matrix $A_N(\tau)$ are given by [CGLW16, Eq. (5.20)]

$$ (2.16) \qquad \left( A_N(\tau) \right)_{i,j} = \frac{1}{\|\mathcal{K}_i\|_{\mathcal{E}}^2} \left( a + bK_j(-1) + c\tau(2\delta_{j,0} - 1) \right. \\ \left. \qquad + \frac{2}{\tau} \sum_{k=0}^{j-1} a_{j,k} (\delta_{i,k} \|\mathcal{K}_i\|_{\mathcal{E}}^2 - 1) \right), \quad i,j = 0, \dots, N-1. $$

Here, the coefficients $a_{j,k}$ are independent of the DDE model and are determined by [CGLW16, Prop. 5.1]; see Proposition 1 in Appendix A. Re-writing $A_N(\tau)$ as $A_N(\tau) = 2/\tau P_N + Q_N$, one observes that only $Q_N$ depends on the model's parameters. More generally, the matrix $A_N(\tau)$ accounts for the approximation of the transport equation $\partial_t u = \partial_\theta u$ and the contribution of the linear terms arising in the nonlocal boundary condition (2.11); see [CGLW16].

The nonlinear mapping $G_2$ has its components $F_N^j$ given by [CGLW16, Eq. (5.21)]

$$ (2.17) \qquad F_N^j(\mathbf{y}, \tau) = \frac{1}{\|\mathcal{K}_j\|_{\mathcal{E}}^2} F \left( \sum_{n=0}^{N-1} y_n(t), \sum_{n=0}^{N-1} y_n(t)K_n(-1), \tau y_0(t) - \tau \sum_{n=1}^{N-1} y_n(t) \right), $$

for each $0 \le j \le N-1$.

Due to rigorous convergence results of GK approximations [CGLW16], the GK formulas above provide a powerful apparatus to analyze DDEs by means of ODE approximations.

### 3. CENTER-UNSTABLE MANIFOLD APPROXIMATIONS FROM GK SYSTEMS

We survey in this section classical techniques of approximations of center-unstable manifolds of a steady state for systems of autonomous ordinary differential equations (ODEs) in $\mathbb{R}^N$, for which the right-hand side (RHS) is the RHS of a GK system such as given by (2.15).
3.1. **The leading-order approximation theorem.** Invariant manifold theory allows for the rigorous derivation of low-dimensional surrogate systems from which not only the system's qualitative behavior near e.g. a steady state is preserved, but also quantitative features of the nonlinear dynamics are reasonably well approximated such as the solution's amplitude or possible dominant periods.

As a common practice in invariant manifold theory and in view of applications considered in Sec. 4, we work with the perturbed equation of Eq. (2.1) (such as Eq. (4.4) below) about a steady state, and its GK approximation of the form

$$ (3.1) \qquad \frac{d\mathbf{y}}{dt} = A(\tau)\mathbf{y} + G(\mathbf{y}), $$

dropping the dependence on $N$. We assume that nonlinear mapping, $G: \mathbb{R}^N \to \mathbb{R}^N$, satisfies the following tangency condition

$$ (3.2) \qquad G(0) = 0, \quad \text{and} \quad DG(0) = 0. $$

The reduction formulas presented below, are derived for the case where $G$ does not depend on $\tau$. This choice is made to simplify the notations. The general case where $G$ depends on $\tau$ (see (2.15)) leads to the same type of formulas as long as the tangency condition (3.2) is satisfied for all $\tau$.

Assuming that $G$ is sufficiently smooth, then $G(\mathbf{y})$ admits the following Taylor expansion for $\mathbf{y}$ near the origin:

$$ (3.3) \qquad G(\mathbf{y}) = G_k(\underbrace{\mathbf{y}, \dots, \mathbf{y}}_{k \text{ times}}) + O(\|\mathbf{y}\|^{k+1}), $$

where

$$ (3.4) \qquad G_k: \underbrace{\mathbb{R}^N \times \dots \times \mathbb{R}^N}_{k \text{ times}} \to \mathbb{R}^N $$

denotes a homogenous polynomial of order $k \ge 2$. Often, $G_k(\mathbf{y})$ is used hereafter as a compact notation for $G_k(\mathbf{y}, \dots, \mathbf{y})$. We label the elements of the spectrum $\sigma(A(\tau))$ of $A(\tau)$ according to the lexicographical order. According to this rearrangement, the eigenvalues are labeled by a single positive integer $n$, so that

$$ (3.5) \qquad \sigma(A(\tau)) = \{\lambda_n(\tau), 1 \le n \le N\}, $$

with, for any $1 \le n < n'$, either

$$ (3.6) \qquad \operatorname{Re} \lambda_n(\tau) > \operatorname{Re} \lambda_{n'}(\tau), $$

or

$$ (3.7) \qquad \operatorname{Re} \lambda_n(\tau) = \operatorname{Re} \lambda_{n'}(\tau), \quad \text{and} \quad \operatorname{Im} \lambda_n(\tau) \ge \operatorname{Im} \lambda_{n'}(\tau). $$

In this convention, an eigenvalue of algebraic multiplicity $m_c$, is repeated $m_c$ times.

We are concerned with describing how linear instabilities translate to the nonlinear dynamics. To do so, the onset of instability is described in terms of the Principle of Exchange of Stability (PES) [MW05], concerned with the loss of stability of the basic steady state. More precisely, the PES describes situations for which the spectrum of $A(\tau)$ experiences the following change at a critical parameter $\tau_c$:

$$ (3.8) \qquad \begin{array}{ll} \operatorname{Re} \lambda_j(\tau) & \begin{cases} < 0 & \text{if } \tau < \tau_c, \\ = 0 & \text{if } \tau = \tau_c, \quad 1 \le j \le m_c, \\ > 0 & \text{if } \tau > \tau_c, \end{cases} \\ \\ \operatorname{Re} \lambda_j(\tau_c) < 0, & j \ge m_c + 1, \end{array} $$

for some $m_c > 0$, and for $\tau$ in some neighborhood $U$ of $\tau_c$. Of course, the PES holds also when one crosses $\tau_c$ from above with $\operatorname{Re} \lambda_j(\tau) < 0$ for $\tau > \tau_c$, while $\operatorname{Re} \lambda_j(\tau) > 0$ when $\tau < \tau_c$.
Whatever the situation (destabilization by crossing from above or below a critical value), one associates to the PES condition the following decomposition of the spectrum $\sigma(A(\tau))$:

$$ (3.9) \qquad \begin{aligned} \sigma_c(A(\tau)) &= \{\lambda_j(\tau) \mid j = 1, \dots, m_c\}, \\ \sigma_s(A(\tau)) &= \{\lambda_j(\tau) \mid j = m_c + 1, \dots, N\}. \end{aligned} $$

The PES condition prevents other eigenvalues in $\sigma_s(A(\tau))$ from crossing the imaginary axis as $\tau$ varies in $\mathcal{U}$. Hence, no eigenvalues other than those of $\sigma_c(A(\tau))$ change sign in the neighborhood $\mathcal{U}$. Furthermore, the PES condition implies, by a continuity argument, the following uniform spectral gap by possibly reducing $\mathcal{U}$ accordingly [CLW15a, Lemma 6.1],

$$ (3.10) \qquad 0 > 2k\eta_c > \eta_s, $$

where $k$ is the leading order of $G$, and

$$ \eta_c = \inf_{\tau \in \mathcal{U}} \inf_{j=1, \dots, m_c} \mathrm{Re}(\lambda_j(\tau)), \quad \eta_s = \sup_{\tau \in \mathcal{U}} \sup_{j \ge m_c+1} \mathrm{Re}(\lambda_j(\tau)). $$

To the modes that lose their stability according to the PES, we associate the reduced state space, $H_c$, given by

$$ (3.11) \qquad H_c = \operatorname{span}\{e_1, \dots, e_{m_c}\}, $$

while a mode $e_n$ with $n \ge m_c + 1$ denotes a stable mode. Throughout this article we chose not to make explicit the $\tau$-dependence of the eigenmodes but this dependence should be kept in mind. We denote by $H_s$ the subspace spanned by these stable modes. The projector onto the subspace $H_s$ (resp. $H_c$) is denoted by $\Pi_s$ (resp. $\Pi_c$). Similarly, $A_c(\tau)$ (resp. $A_s(\tau)$) denotes $\Pi_c A(\tau)$ (resp. $\Pi_s A(\tau)$), and $y_c$ (resp. $y_s$) denotes the vector in $H_c$ (resp. $H_s$) of the low-mode amplitudes (resp. stable-mode amplitudes). The inner product in $\mathbb{C}^N$ is denoted by $\langle \cdot, \cdot \rangle$ and is defined by

$$ (3.12) \qquad \langle \boldsymbol{a}, \boldsymbol{b} \rangle = \sum_{i=1}^{N} a_i b_i^*, \quad \boldsymbol{a}, \boldsymbol{b} \in \mathbb{C}^N. $$

In what follows we also denote by $e_j^*$ the eigenmodes of the conjugate transpose $A(\tau)^*$ of $A(\tau)$.

Condition (3.10) ensures in particular that the following spectral gap holds for $\tau$ in $\mathcal{U}$

$$ \gamma_{m_c}(\tau) = \mathrm{Re}(\lambda_{m_c}(\tau)) - \mathrm{Re}(\lambda_{m_c+1}(\tau)) > 0. $$

It is well known that the existence of a (local) exact parameterization or say in other words, of a local $m_c$-dimensional invariant manifold, is subject to the following spectral gap condition:

$$ (3.13) \qquad \gamma_{m_c}(\tau) \ge CL\mathrm{Lip}(G|_{\mathcal{V}}), $$

where $\mathrm{Lip}(G|_{\mathcal{V}})$ denotes the Lipschitz constant of the nonlinearity $G$ restricted to some neighborhood $\mathcal{V}$ of the origin in $\mathbb{C}^N$, and $C > 0$ is typically independent on $\mathcal{V}$. Due to the tangency condition (3.2), condition (3.13) always holds once $\mathcal{V}$ is chosen sufficiently small. The theory of local invariant manifolds makes thus sense if solutions with sufficiently small amplitudes lie within the appropriate neighborhood $\mathcal{V}$.

In the context of e.g. nonlinear oscillations that bifurcate from a steady state, the local invariant manifold provides an exact parameterization<fnref n="1" /> of the stable limit cycle near criticality in the case of a supercritical Hopf bifurcation. In the case of subcritical Hopf bifurcation, it provides the parameterization of the unstable limit cycle that emerges in a continuous fashion from the steady state. In Sec. 4, we show that the approximation formulas provided by Theorem 3.1 and in Sec. 3.2 below may allow for approximating not only such unstable limit cycles that unfold continuously from the origin but also the stable limit cycles that are distant from the origin, corresponding to a “jump” transition [MW14] associated with a subcritical Hopf bifurcation.

<footnote n="1">¹As provided for instance by a center manifold or the unstable manifold of the origin.</footnote>
**Theorem 3.1.** Assume that $G$ and $A(\tau)$ satisfy the assumptions recalled above, and that the PES condition (3.8) is satisfied. Then for each $\tau$ in a neighborhood $\mathcal{U}$ of $\tau_c$, Eq. (3.1) admits a local invariant manifold, $\mathfrak{M}_\tau = \text{graph}(h_\tau)$, with $h_\tau$ that maps $H_c$ into the stable subspace $H_\mathfrak{s}$ such that $h_\tau(0) = 0$ and $Dh_\tau(0) = 0$.

Assume that the following non-resonance condition is satisfied:

$$ (NR) \qquad \forall (j_1, \dots, j_k) \in (1, \dots, m_c)^k, n \ge m_c + 1, \text{ it holds that} \\ \qquad \qquad \left( \langle G_k(e_{j_1}, \dots, e_{j_k}), e_n^* \rangle \ne 0 \right) \implies \left( \text{Re} \left( \sum_{\ell=1}^k \lambda_{j_\ell} - \lambda_n \right) \ne 0 \right), $$

where $G_k$ denotes the leading-order term in the Taylor expansion of $G$.

Then, the Lyapunov-Perron integral, $\Phi_\tau(\mathbf{X}) = \int_{-\infty}^0 e^{-sA_\mathfrak{s}(\tau)} \Pi_\mathfrak{s} G_k(e^{sA_c(\tau)}\mathbf{X}) \mathrm{d}s$, is well defined and is a solution to the following homological equation

$$ (3.14) \qquad \mathcal{L}_{A(\tau)}[\psi](\mathbf{X}) = \Pi_{\mathfrak{s}} G_k(\mathbf{X}), $$

where $\mathcal{L}_A$ denotes the differential operator acting on smooth mappings $\psi$ from $H_c$ into $H_\mathfrak{s}$, defined as:

$$ (3.15) \qquad \mathcal{L}_{A(\tau)}[\psi](\mathbf{X}) = D\psi(\mathbf{X})A_c(\tau)\mathbf{X} - A_\mathfrak{s}(\tau)\psi(\mathbf{X}), \quad \mathbf{X} \in H_c. $$

Moreover, $\Phi_\tau(\mathbf{X})$ provides the leading-order approximation of the invariant manifold function $h_\tau$ characterizing $\mathfrak{M}_\tau$, in the sense that

$$ (3.16) \qquad \|h_\tau(\mathbf{X}) - \Phi_\tau(\mathbf{X})\| = o(\|\mathbf{X}\|^k), $$

as long as $\mathbf{X}$ lies in a neighborhood $\mathcal{N}(\tau)$ of the origin in $H_c$ spanned by the $m_c$ eigenmodes $e_j$ losing their stability (see (3.11)), as $\tau$ crosses $\tau_c$.

Finally, $\Phi_\tau(\mathbf{X})$ possesses the following explicit formula:

$$ (3.17) \qquad \Phi_\tau(\mathbf{X}) = \sum_{n \ge m_c + 1} \Phi_{n,\tau}(\mathbf{X}) e_n, \quad \mathbf{X} \in H_c, $$

where

$$ (3.18) \qquad \Phi_{n,\tau}(\mathbf{X}) = \sum_{1 \le j_1, \dots, j_k \le m_c} G_{j_1 \dots j_k}^n M_{j_1 \dots j_k}^{n,\tau} X_{j_1} \dots X_{j_k}, \quad X_j = \langle \mathbf{X}, e_{j_1} \rangle, $$

with

$$ (3.19) \qquad M_{j_1 \dots j_k}^{n, \tau} = \left( \sum_{\ell=1}^{k} \lambda_{j_\ell}(\tau) - \lambda_n(\tau) \right)^{-1}, $$

and the $G_{j_1 \dots j_k}^n$ denoting the coefficients accounting for the magnitude carried out by $e_n^*$, of the nonlinear interactions through the leading-order term $G_k$, between the low modes $e_{j_1}, \dots, e_{j_k}$ (in $H_c$), namely:

$$ (3.20) \qquad G_{j_1 \dots j_k}^n = \langle G_k(e_{j_1}, \dots, e_{j_k}), e_n^* \rangle, \quad 1 \le j_1, \dots, j_k \le m_c. $$

Conditions similar to (NR) arise in the smooth linearization of dynamical systems near an equilibrium [Sel85]. Here, condition (NR) implies that the eigenvalues of the stable part satisfy a Sternberg condition of order $k$ [Sel85] with respect to the eigenvalues associated with the modes spanning the reduced state space $H_c$.

This theorem is essentially a consequence of [CLM20, Theorems 1 and 2], in which condition (NR) is a stronger version of that used for [CLM20, Theorem 2]; see also [CLM20, Remark 1 (iv)]. This condition is necessary and sufficient here for $\int_{-\infty}^0 e^{-sA_\mathfrak{s}(\tau)} \Pi_\mathfrak{s} G_k(e^{sA_c(\tau)}\mathbf{X}) \mathrm{d}s$ to be well defined. For the derivation of (3.14) see that of Eq. (4.6) in [CLM20]. The error estimate (3.16) is a corollary of the more general results [CLW15a, Theorem 6.1] and [CLW15a, Corollary 6.1] proved for stochastic Partial Differential Equations (PDEs). See also [Hen81, Thm. 6.2.3] and [Hen81, Lemma 6.2.4].

The assumptions of Theorem 3.1 ensure that for any $\tau$ in $\mathcal{U}$, the following reduction principle holds:
(i) Any solution $\mathbf{y}(t)$ to Eq. (3.1) such that $\mathbf{y}(t_0)$ belongs to $\mathfrak{M}_\tau$ for some $t_0$, stays on $\mathfrak{M}_\tau$ over an interval of time $[t_0, t_0 + \sigma)$, $\sigma > 0$, i.e.

$$ (3.21) \qquad \mathbf{y}(t) = \mathbf{y}_c(t) + h_\tau(\mathbf{y}_c(t)), \quad t \in [t_0, t_0 + \sigma), $$
where $\mathbf{y}_c(t)$ denotes the projection of $\mathbf{y}(t)$ onto the subspace $H_c$.

(ii) If there exists a trajectory $t \mapsto \mathbf{y}(t)$ such that $\mathbf{y}_c(t)$ belongs, for all $-\infty < t < \infty$, to the neighborhood $\mathfrak{B}(\tau)$ (in $H_c$) over which $\mathfrak{M}_\tau$ is well defined, then the trajectory must lie on $\mathfrak{M}_\tau$.

Such a formulation of the reduction principle can be inferred for instance from [SY02, Theorem 71.4] and [Cra91]. Property (ii) implies that an invariant set $\Sigma$ to Eq. (3.1) of any type, e.g., equilibria, periodic orbits, invariant tori, must lie in $\mathfrak{M}_\tau$ if its projection onto $H_c$ is contained in $\mathfrak{B}(\tau)$, i.e. if $\Pi_c\Sigma \subset \mathfrak{B}(\tau)$. Property (3.21) holds then globally in time for the solutions that composed such invariant sets, and thus the knowledge of the $m_c$-dimensional variable, $\mathbf{y}_c(t)$, is sufficient to entirely determine any solution $\mathbf{y}(t)$ that belongs to such an invariant set; see [SY02, Theorem 71.4]. Furthermore, $\mathbf{y}_c(t)$ is obtained as the solution of the following reduced $m_c$-dimensional problem

$$ (3.22) \qquad \hat{\mathbf{y}} = A_c(\tau)\hat{\mathbf{y}} + \Pi_c G(\hat{\mathbf{y}} + h_\tau(\hat{\mathbf{y}})), \quad \hat{\mathbf{y}}(0) = \mathbf{y}_c(0) \in \mathfrak{B}(\tau), $$

which in turn characterizes the solution $\mathbf{y}(t)$ in $\Sigma$, since the slaving relationship $\mathbf{y}_s(t) = h_\tau(\mathbf{y}_c(t))$ holds for any solution $\mathbf{y}(t)$ that belongs to an invariant set $\Sigma$ for which $\Pi_c\Sigma \subset \mathfrak{B}(\tau)$.

Based on the approximation Theorem 3.1 and the reduction principle, it is thus reasonable to anticipate that the following effective reduced equation,

$$ (3.23) \qquad \dot{x} = A_c(\tau)x + \Pi_c G_k(x + \Phi_\tau(x)), \quad x(0) = \mathbf{y}_c(0) \in \mathfrak{B}(\tau), $$

provides an approximation of the invariant set $\Sigma$ for which $\Pi_c\Sigma \subset \mathfrak{B}(\tau)$. Since Eq. (3.23) is derived in practice from a GK system, it will be referred to as an *effective reduced GK system* or *mD reduced GK system*, where $m = \dim(H_c)$. Eq. (3.23) is thus particularly suited for effective approximations of invariant sets to Eq. (3.1), made of solutions of sufficiently small amplitudes, and that bifurcate when one crosses a critical parameter $\tau_c$ that destabilizes an equilibrium point. Due to the rigorous convergence results as $N \to \infty$ to DDE solutions from solutions to GK systems [CGLW16], it is expected that the invariant set made of bifurcating solutions to the original DDE (2.1), can also be well approximated from the solutions of the effective reduced system (3.23), when $N$ is sufficiently large. In some applications, one may need though to go to higher-order approximations than the leading-order one in order to obtain more accurate approximations of $\hat{\mathbf{y}}(t)$ (and thus of $\mathbf{y}(t)$). This point is described next. Numerical evidences of such efficient approximations are then presented in Sec. 4 below; see e.g. Fig. 5.

**Remark 3.1.** In certain applications, we can encounter that, although designed for $\tau$ close to a critical parameter value $\tau_c$, the effective reduced system (3.23) may show skills in approximating bifurcating solutions to the DDE with not necessarily small amplitudes (order 1), as $\tau$ is crossing other critical values. Sections 4.3 and 4.4 below provide numerical illustrations of such situations in the case of a delay model of ENSO.

**3.2. Analytic formulas for higher-order approximations.** We provide here a simple derivation of higher-order approximations of an invariant manifold. The *invariance equation* is insightful in that respect. It is the equation satisfied by the local invariant manifold, $h_\tau$, to Eq. (3.1) in a neighborhood of the origin, namely

$$ (3.24) \qquad Dh(\mathbf{X})[A_c\mathbf{X} + \Pi_c G(\mathbf{X} + h(\mathbf{X}))] - A_s h(\mathbf{X}) = \Pi_s G(\mathbf{X} + h(\mathbf{X})), $$

in which we have dropped the dependence on $\tau$ to simplify the notations; see [Hen81, Corollary 6.2.2].

This functional equation is a nonlinear system of first order PDEs that cannot be solved in closed form except in special cases. Many methods exist to seek for approximate solutions to Eq. (3.24); see e.g. [BKF98, EvP04, HCF+16]. Here, we adopt a standard use of power series expansions that we blend with energy content arguments to neglect certain terms and thus simplify the expressions. Indeed, instead of keeping all the monomials at a given degree arising from such an expansion, we filter out terms that carries significantly less
energy compared with those that are kept. This elimination procedure relies on the observation that, close to criticality, the projected ODE dynamics onto the resolved subspace $H_c$ contains typically a large fraction of the solution's total energy. We illustrate this idea to the case where $G(\mathbf{y}) = G_2(\mathbf{y}, \mathbf{y}) + G_3(\mathbf{y}, \mathbf{y}, \mathbf{y})$ and a cubic approximation is sought. The higher-order cases proceed in a same fashion and are omitted here for the sake of concision. We refer to [HCF$^+$16] for mathematical details about higher-order approximations of invariant manifolds.

To determine the third-order approximation, we replace $h$ in the invariance equation (3.24) by $\Phi_{\tau} + \psi$, where $\psi$ denotes the homogeneous cubic terms in the power expansion of $h$, to be determined. By identifying the terms of order two, we recover Eq. (3.14) satisfied by $\Phi_{\tau}$. By identifying the terms of order three, we obtain the following equation for $\psi$:

$$ (3.25) \qquad \mathcal{L}_A[\psi](\mathbf{X}) = -D\Phi_{\tau}(\mathbf{X})\Pi_c G_2(\mathbf{X}) + \Pi_s G_2(\mathbf{X}, \Phi_{\tau}(\mathbf{X})) + \Pi_s G_2(\Phi_{\tau}(\mathbf{X}), \mathbf{X}) + \Pi_s G_3(\mathbf{X}), $$
where $\mathcal{L}_A[\psi](\mathbf{X}) = D\psi(\mathbf{X})A_c\mathbf{X} - A_s\psi(\mathbf{X})$; see (3.15).

Note that the RHS is a homogeneous cubic polynomial in the $X$-variable. As recalled above, close to criticality, a large fraction of the energy is contained in the low modes and therefore the energy carried by $y_s$ is typically much smaller than $\|y_c\|^2$. It is then reasonable to expect that the energy carried by $\Phi_{\tau}(\mathbf{X})$ is much smaller than $\|\mathbf{X}\|^2$ for $\mathbf{X} = y_c$. Thus, one may assume that in the RHS of (3.25), the term $\Pi_s G_3(\mathbf{X})$ dominates the other three terms provided that $\|G_2(\mathbf{y})\|/\|\mathbf{y}\|^2$ is on the same order of magnitude as $\|G_3(\mathbf{y})\|/\|\mathbf{y}\|^3$. As a consequence, it is reasonable to hope for reasonably accurate approximations of $\psi$ with $h_3$ by simply solving the equation:

$$ (3.26) \qquad Dh_3(\mathbf{X})A_c\mathbf{X} - A_s h_3(\mathbf{X}) = \Pi_s G_3(\mathbf{X}). $$

Note that this equation is exactly Eq. (3.14) with $k=3$. In virtue of Theorem 3.1, the existence of $h_3$ is guaranteed by non-resonance condition (NR) (with $k=3$), and $h_3$ is thus given by (3.17)–(3.18) with $k=3$. We arrive then at the following “high-mode” parameterization

$$ (3.27) \qquad \begin{aligned} \Psi_{\tau}(\mathbf{X}) = \Phi_{\tau}(\mathbf{X}) + h_3(\mathbf{X}) &= \sum_{n=m_c+1}^{N} \Psi_{n,\tau}(X)e_n \\ &= \sum_{n=m_c+1}^{N} \left( \begin{aligned} &\sum_{1 \le j_1, j_2 \le m_c} \frac{\langle G_2(e_{j_1}, e_{j_2}), e_n^* \rangle}{\lambda_{j_1}(\tau) + \lambda_{j_2}(\tau) - \lambda_n(\tau)} X_{j_1} X_{j_2} \\ &+ \sum_{1 \le j_1, j_2, j_3 \le m_c} \frac{\langle G_3(e_{j_1}, e_{j_2}, e_{j_3}), e_n^* \rangle}{\lambda_{j_1}(\tau) + \lambda_{j_2}(\tau) + \lambda_{j_3}(\tau) - \lambda_n(\tau)} X_{j_1} X_{j_2} X_{j_3} \end{aligned} \right) e_n. \end{aligned} $$

In Sec. 4 below, we show that such parameterizations enable us to reach accurate approximations of bifurcating solutions in the case of a DDE with cubic nonlinearity.

## 4. BIFURCATIONS ANALYSIS OF A DELAY EL NIÑO-SOUTHERN OSCILLATION (ENSO) MODEL

### 4.1. GK approximations of the Suarez and Schopf ENSO model.
The Suarez and Schopf model is given by the following DDE [SS88]

$$ (4.1) \qquad \frac{dT}{dt} = T(t) - \alpha T(t - \tau) - T^3(t), $$

where the unknown $T$ represents the non-dimensionalized sea surface temperature (SST) anomalies at the eastern equatorial Pacific Ocean, $\tau$ and $\alpha$ are positive constants, and the physically relevant range of $\alpha$ used in [SS88] is (0, 1). For this given range of $\alpha$, Eq. (4.1) admits three fixed points:

$$ (4.2) \qquad T_0 = 0, \quad T_+ = \sqrt{1 - \alpha}, \quad T_- = -\sqrt{1 - \alpha}. $$
The first term in Eq. (4.1) reflects an instantaneous positive feedback, whereby an SST perturbation heats the atmosphere, whose wind response drives ocean currents to reinforce the original perturbation. The second term in Eq. (4.1) accounts for the transit time of equatorially trapped oceanic waves to cross the Pacific ocean [BTR07]. The third term accounts then for saturation effects which limit the instability growth due to the positive feedback by effects tied to advective processes in the ocean and moist processes in the atmosphere. A linear stability analysis of this simple model was conducted in [SS88] to show that the steady-state solution loses stability for certain parameter values. The resulting periodic solutions have periods of at least twice the length of delay, supporting that this simple feedback mechanism can be, in theory, consistent with ENSO's oscillatory behaviour on an interannual timescale. In Sec. 5 below, we analyze in greater details the nonlinear structures at play that organize the time-variability of solutions to the Suarez and Schopf model when subject to noise disturbances. These nonlinear structures are invariant sets that emerge through local and global bifurcations.

In what follows, we conduct thus a bifurcation analysis of the Suarez and Schopf model. To do so, we apply the GK approximation framework of Sec. 2 and the center-unstable manifold reduction formulas of Sec. 3. In that respect, to fit within the framework of Sec. 3, we first derive from Eq. (4.1), the DDE satisfied by the perturbed variable about the steady state $T_+$,

$$ (4.3) \qquad \hat{T} = T - T_+, $$

namely,

$$ (4.4) \qquad \frac{d\hat{T}}{dt} = (1 - 3T_+^2) \hat{T}(t) - \alpha \hat{T}(t - \tau) - 3T_+ \hat{T}^2(t) - \hat{T}^3(t). $$

The above equation fits into Eq. (2.1) with

$$ (4.5) \qquad a = 1 - 3T_+^2, \quad b = -\alpha, \quad c = 0, \text{ and } F(\hat{T}(t)) = -3T_+ \hat{T}^2(t) - \hat{T}^3(t). $$

By applying Eqns. (2.15)–(2.17) to Eq. (4.4), we obtain the following N-dimensional GK system

$$ (4.6) \qquad \frac{d\mathbf{y}}{dt} = A(\tau)\mathbf{y} + G(\mathbf{y}), \qquad \mathbf{y} = (y_1, \dots, y_N)^T, $$

where the entries of the matrix $A(\tau)$ are given here by

$$ (4.7) \qquad \begin{aligned} (A(\tau))_{j,n} = \frac{1}{\|\mathcal{K}_j\|_{\varepsilon}^2} & \sum_{n=0}^{N-1} (1 - 3T_+^2 - \alpha K_n(-1) \\ & + \frac{2}{\tau} \sum_{k=0}^{n-1} a_{n,k} (\delta_{j,k} \|\mathcal{K}_j\|_{\varepsilon}^2 - 1)), \quad \text{for } j, n = 0, \dots, N-1. \end{aligned} $$

As in Sec. 2.3, the $a_{n,k}$ are determined thanks to Proposition 1 recalled in Appendix A, and the $\|\mathcal{K}_j\|_{\varepsilon}$ and $K_n(-1)$ are given by (2.6) and (A.1), respectively.

The nonlinearity $G$ is given as the following sum of N-dimensional mappings

$$ G(\mathbf{y}) = G_2(\mathbf{y}) + G_3(\mathbf{y}), $$

with

$$ (4.8) \qquad G_2(\mathbf{y}) = -3T_+ \left( \sum_{n=0}^{N-1} y_n \right)^2 \nu_N, \quad \text{and} \quad G_3(\mathbf{y}) = - \left( \sum_{n=0}^{N-1} y_n \right)^3 \nu_N, $$

where $\nu_N$ denotes the N-dimensional column vector

$$ (4.9) \qquad \nu_N = \left( \frac{1}{\|\mathcal{K}_0\|_{\varepsilon}^2}, \dots, \frac{1}{\|\mathcal{K}_{N-1}\|_{\varepsilon}^2} \right)^T $$
We arrive finally at the following N-dimensional GK approximation for the DDE (4.4)

$$ (4.10) \qquad \frac{d\mathbf{y}}{dt} = \left(\frac{2}{\tau}\mathcal{P}_N + \mathcal{Q}_N\right) \mathbf{y} - \left(\sum_{n=1}^{N} y_n\right)^2 \left(3T_+ + \sum_{n=1}^{N} y_n\right) \boldsymbol{\nu}_N. $$

We refer to Appendix A for the computation of $\mathcal{P}_N$ and $\mathcal{Q}_N$. By construction, the nonlinear term $G$ satisfies the tangency condition (3.2) and thus one can unpack the center-unstable manifold framework of Sec. 3. This is done in Sec. 4.2 below. To prepare for the resulting bifurcation analysis, we consider the case $\alpha = 0.75$. For this value of $\alpha$, the steady state is given by $T_+ = \sqrt{1 - \alpha} = 0.5$.

We conclude this section with Table 1 that shows the $L^\infty$-error achieved by various N-dimensional GK systems in approximating stable limit cycles to the DDE (4.4) for $\tau$-values that will be considered later on.

<table><caption>**Table 1.** $L^\infty$-error achieved by the N-dimensional GK approximations to the DDE (4.4). Here, the $L^\infty$-errors are computed over one period of the stable limit cycle for each tested $\tau$-value, with the DDE (4.4) initialized using a segment on the stable limit cycle, and the GK systems initialized using the projection of the DDE initial data onto the first N Koornwinder basis functions. The timestep size is set to $\delta t = \tau/2^{18}$.</caption><thead><tr><th>GK dimension</th><th>&tau; = 1.562</th><th>&tau; = 1.6</th><th>&tau; = 1.7</th><th>&tau; = 1.9</th></tr></thead><tbody><tr><td>N = 4</td><td>5.66 &times; 10<sup>-2</sup></td><td>3.87 &times; 10<sup>-2</sup></td><td>4.29 &times; 10<sup>-2</sup></td><td>5.71 &times; 10<sup>-2</sup></td></tr><tr><td>N = 6</td><td>1.44 &times; 10<sup>-4</sup></td><td>1.80 &times; 10<sup>-4</sup></td><td>7.75 &times; 10<sup>-4</sup></td><td>3.60 &times; 10<sup>-3</sup></td></tr><tr><td>N = 8</td><td>1.54 &times; 10<sup>-4</sup></td><td>6.95 &times; 10<sup>-5</sup></td><td>2.13 &times; 10<sup>-5</sup></td><td>4.37 &times; 10<sup>-4</sup></td></tr><tr><td>N = 10</td><td>1.31 &times; 10<sup>-4</sup></td><td>5.60 &times; 10<sup>-5</sup></td><td>3.65 &times; 10<sup>-5</sup></td><td>8.13 &times; 10<sup>-5</sup></td></tr></tbody></table>

## 4.2. Subcritical Hopf Bifurcation.
In this section, we apply the center-unstable manifold reduction framework of Sec. 3 to the the GK system (4.10), and in particular the high-mode parameterization $\Psi_\tau$ given by (3.27), to derive our effective reduced GK system and conduct a bifurcation analysis of the (perturbed) Suarez and Schopf model (4.4).

Thus, our $m_c$-dimensional effective reduced GK system based on $\Psi_\tau$ writes in compact form

$$ (4.11) \qquad \dot{x} = A_c(\tau)x + \Pi_c G(x + \Psi_\tau(x)), $$

for which $m_c$ is determined by analyzing the modes that destabilize once a critical delay parameter $\tau$ is crossed.

In the eigenbasis coordinate system of $A(\tau)$, the GK system (4.10) writes

$$ (4.12) \qquad \dot{y}_j = \lambda_j(\tau)y_j + \langle G\left(\sum_{\ell=1}^{N} y_\ell e_\ell\right), e_j^* \rangle, \quad j = 1, \dots, N, $$

where $y_j = \langle \mathbf{y}, e_j^* \rangle$, with inner product defined in (3.12).

We observe that a dominant pair of complex eigenvalues crosses the imaginary axis as $\tau$ crosses from below the critical value $\tau_c \approx 1.74$ (see Table 2), whereas the other pairs of eigenvalues stay within the left half complex plane, with a clear spectral gap, even for $\tau$ above the critical value; see Fig. 2. We choose thus $m_c = 2$ and the subspaces $H_c$ and $H_s$ as follows:

$$ (4.13) \qquad H_c = \operatorname{span}\{e_1, e_2\}, \qquad H_s = \operatorname{span}\{e_3, \dots, e_N\}, $$

in which $(e_1, e_2)$ denotes the conjugate pair of modes that destabilize as $\tau$ crosses $\tau_c$. Setting $m_c = 2$ in (3.27), the 2D reduced GK system (4.11) (in $\mathbb{C}^2$) writes in the eigenbasis coordinates:

$$ (4.14) \qquad \dot{x}_j = \lambda_j(\tau)x_j + \langle G(x + \Psi_\tau(x)), e_j^* \rangle, \quad j = 1, 2, $$

with $x = x_1e_1 + x_2e_2$, and where the parameterization $\Psi_\tau$ of the neglected modes in $H_s$ is given by (3.27).
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Eigenvalues dependence as $\tau$ crosses its critical value $\tau_c$. Are shown here, the first 10 pairs $(\lambda_j(\tau), \bar{\lambda}_j(\tau))$ as $\tau$ is increased from $\tau = 1.3$ to $\tau = 2.5$, for $\alpha = 0.75$. These pairs move all from left to right. The pair of eigenvalues that crosses the imaginary axis for the critical delay parameter $\tau = \tau_c$ is shown in red, while the stable pairs are shown in black. These eigenvalues are computed from the GK linear part $A(\tau)$ given in (4.7) in dimension $N = 50$. It is noteworthy that these 10 pairs of GK eigenvalues (for $N = 50$) satisfy the actual characteristic equation associated with the linear part of the DDE (4.4), $\lambda = (1 - 3T_+^2) - \alpha e^{-\lambda\tau}$, up to a maximal error of $10^{-4}$, as $\tau$ varies from $\tau = 1.3$ to $\tau = 2.5$.</figcaption></figure>

The behavior of the eigenvalues shown in Fig. 2 is an illustration of the PES condition (3.8) for the $N$-dimensional GK approximation (4.12) of the DDE (4.4). In particular, the critical crossing of the dominant pair (red curves in Fig. 2) from the left to the right half plane is a manifestation of the well-known transversality condition, for the underlying GK system to admit a Hopf bifurcation [CKL20, Eqns. (32)-(33)]. Theorem III.1 of [CKL20] provides then the precise conditions for a Hopf bifurcation to take place for the GK system and ensures that its type, either supercritical or subcritical, is determined by calculating the Lyapunov coefficient from the 2D reduced GK system (4.14) put into its Stuart-Landau (SL) form [CKL20, Eq. (36)]. As examined in [CKL20], the same type of Hopf bifurcation such as predicted by the SL equation, occurs then for the DDE as well, as long as the coefficients of the SL equation are calculated from a sufficiently large dimensional GK system.

In that respect, we determine the Lyapunov coefficient as given by the analytic formula [CKL20, Eq. (40)] from the coefficients of the $N$-dimensional GK approximation (4.12) of the DDE (4.4). To this end, we rewrite $G_2$ and $G_3$ given in (4.8) as the following bilinear and trilinear forms:

$$ (4.15) \qquad \begin{aligned} G_N^{(2)}(\mathbf{p}, \mathbf{q}) &= -3T_+ \left(\sum_{n=0}^{N-1} p_n\right) \left(\sum_{n=0}^{N-1} q_n\right) \boldsymbol{\nu}_N \\ G_N^{(3)}(\mathbf{p}, \mathbf{q}, \mathbf{r}) &= -\left(\sum_{n=0}^{N-1} p_n\right) \left(\sum_{n=0}^{N-1} q_n\right) \left(\sum_{n=0}^{N-1} r_n\right) \boldsymbol{\nu}_N. \end{aligned} $$

The Lyapunov coefficient $\ell_1^N(\tau_c)$ is then given by

$$ (4.16) \qquad \ell_1^N(\tau_c) = \operatorname{Re}\left(\frac{a_{20}^N a_{11}^N \sqrt{-1}}{\operatorname{Im}(\lambda_1(\tau_c))} + a_{21}^N\right), $$
where

$$ (4.17) \qquad \begin{aligned} a_{20}^N &= \langle G_N^{(2)}(e_1, e_1), e_1^* \rangle, \\ a_{11}^N &= \langle G_N^{(2)}(e_1, e_2), e_1^* \rangle + \langle G_N^{(2)}(e_2, e_1), e_1^* \rangle, \end{aligned} $$

and

$$ (4.18) \qquad \begin{aligned} a_{21}^N &= \langle G_N^{(3)}(e_1, e_1, e_2), e_1^* \rangle + \langle G_N^{(3)}(e_1, e_2, e_1), e_1^* \rangle + \langle G_N^{(3)}(e_2, e_1, e_1), e_1^* \rangle \\ & \quad + \sum_{n=3}^{Nd} \frac{\langle G_N^{(2)}(e_1, e_2), e_n^* \rangle + \langle G_N^{(2)}(e_2, e_1), e_n^* \rangle}{2\mathrm{Re}(\lambda_1(\tau_c)) - \lambda_n(\tau_c)} \left[ \langle G_N^{(2)}(e_1, e_n), e_1^* \rangle + \langle G_N^{(2)}(e_n, e_1), e_1^* \rangle \right] \\ & \quad + \sum_{n=3}^{Nd} \frac{\langle G_N^{(2)}(e_1, e_1), e_n^* \rangle}{2\lambda_1(\tau_c) - \lambda_n(\tau_c)} \left[ \langle G_N^{(2)}(e_2, e_n), e_1^* \rangle + \langle G_N^{(2)}(e_n, e_2), e_1^* \rangle \right]. \end{aligned} $$

Theorem III.1 of [CKL20] ensures then that subcritical Hopf bifurcation occurs for the N-dimensional GK system (4.12) if $\ell_1^N(\tau_c) > 0$ and a supercritical Hopf bifurcation occurs if $\ell_1^N(\tau_c) < 0$.

Table 2 reports on the $\tau$-value $\tau_c$ at which the dominant pair $(\lambda_1(\tau), \overline{\lambda_1(\tau)})$ crosses the imaginary axis, along with its corresponding Lyapunov coefficient $\ell_1^N(\tau_c)$, in terms of the GK dimension $N$. As can be observed, both $\tau_c$ and $\ell_1^N(\tau_c)$ converges quickly as the GK dimension $N$ increases. Since $\ell_1^N(\tau_c)$ converges to a positive value, we infer that the original Suarez and Schopf model (4.4) admits a subcritical Hopf bifurcation for the parameter regime considered here ($\alpha = 0.75$).

The Lyapunov coefficient $\ell_1^N(\tau_c)$ can obviously be computed for other values of $\alpha$ in $(0.5, 1)$, once $\tau_c$ is determined as the critical $\tau$-value at which the dominant pair of eigenvalues of $A(\tau)$ crosses the imaginary axis as in Fig. 2. The results are shown in Fig. 3 for $N = 20$. The positivity of $\ell_1^N(\tau_c)$ shows that the steady state $T_+$ undergoes always a subcritical Hopf bifurcation for $\alpha$ in $(0.5, 1)$. Note that this is the maximum range of $\alpha$ values for which $T_+$ can lose stability, since when $\alpha > 1$, $T_+ = \sqrt{1 - \alpha}$ ceases to exist, and when $\alpha < 0.5$, $T_+$ is always linearly stable for any $\tau$.

One might wonder though whether $\tau_c$ as predicted from the GK systems converges to the actual value of $\tau_c$, obtained by analyzing the characteristic equation associated with the DDE (4.4). We answer this questioning by the affirmative. Recall indeed that one can infer the following analytic expression $\tau_c^{\text{true}}$ given by [BTR07, Eq. (7)] from the characteristic equation,

$$ (4.19) \qquad \tau_c^{\text{true}} = \frac{\arccos\left(\frac{3\alpha-2}{\alpha}\right)}{\sqrt{\alpha^2 - (3\alpha - 2)^2}}. $$

For $\alpha = 0.75$ considered here, we have $\tau_c^{\text{true}} \approx 1.740839502734206$. As shown in Table 2, the 10-dimensional GK system already provides a highly accurate approximation of $\tau_c^{\text{true}}$ with precision $10^{-7}$. A similar precision is obtained for other values of $\alpha$ in $(0.5, 1)$.

<table><caption>Table 2. The estimated $\tau_c$ and Lyapunov coefficient $\ell_1^N(\tau_c)$ as GK dimension $N$ increases.</caption><thead><tr><th></th><th>$N = 4$</th><th>$N = 6$</th><th>$N = 8$</th><th>$N = 10$</th><th>$N = 12$</th></tr></thead><tbody><tr><td>$\tau_c$</td><td>1.7343471</td><td>1.7408640</td><td>1.7408394</td><td>1.7408395</td><td>1.7408395</td></tr><tr><td>$\ell_1^N(\tau_c)$</td><td>2.1884711</td><td>2.2246906</td><td>2.2247594</td><td>2.2247568</td><td>2.2247568</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Lyapunov coefficient $\ell_1^N(\tau_c)$ given in Eq. (4.16): $\alpha$-dependence. Here, the critical value $\tau_c$ and the entries in Eq. (4.16) are computed from a high-dimensional GK approximation ($N = 20$) of the (perturbed) Suarez and Schopf model (4.4).</figcaption></figure>

**4.3. Saddle-Node bifurcation of periodic Orbits (SNO) and homoclinic orbit.** We show in this Section, that the 2D reduced GK system (4.14) is not only useful to predict the subcritical Hopf bifurcation occurring for the Suarez and Schopf model (4.4) at the parameter value $\tau = \tau_c$, but also bifurcations taking place at other critical values of $\tau$. For instance, a simple backward integration of Eq. (4.14) allows for computing a whole family of unstable periodic orbits (UPOs) unfolding from subcritical Hopf bifurcations from $T_+$ and $T_-$ (red points in Fig. 4A) as $\tau$ is further decreased to a value $\tau^\ddagger$ for which the UPOs (dashed black curves in Fig. 4A) hit a saddle equilibrium from both sides, resulting into a homoclinic orbit when $\tau$ reaches a critical value $\tau^\ddagger$; see cyan curve in Fig. 4A. The presence of this homoclinic orbit is manifested by the jump displayed by the curve of UPO's amplitudes<fnref n="2" />, $\tau \mapsto \max(\hat{T}(t)) - \min(\hat{T}(t))$, in the diagram of Fig. 4. This jump is explained by our way of calculating them as these amplitudes are indeed computed over $(\tau^\ddagger, \tau_c)$ only from the branching family of UPOs emanating from $T_+$, i.e. located in the lobe encircling $T_+$ of this homoclinic orbit.

After this jump, as $\tau$ is further decreased from $\tau^\ddagger$, the UPOs encompass the homoclinic orbit (blue curves in Fig. 4B), and their amplitude keeps increasing, until eventually, losing their instability into a Saddle-Node bifurcation of periodic Orbits (SNO bifurcation) for $\tau = \tau^*$. At this critical value the periodic orbit shows mixed stability, with its basin of attraction corresponding to the exterior of the closed curve. Finally, as $\tau$ is increased from $\tau^*$, stable limit cycles unfold with increasing amplitude (red curves in Fig. 4C) and pronounced nonlinear features expressed by their ovaloid shape.

As explained in the next section, this bifurcation diagram obtained from the 2D reduced GK equation (4.14) describes precisely the local and global bifurcations occurring for the (perturbed) Suarez and Schopf model (4.4) due to the remarkable approximation skills shown by the reduced equation.

**4.4. Approximation results of the stable and unstable DDE limit cycles.** We now illustrate numerically that the GK approximation skills coupled with the accuracy of our effective reduced systems as $\tau$ varies, lead to accurate approximations of the UPOs to the DDE that unfold through the subcritical Hopf bifurcation, on one hand, and the stable DDE limit cycles that are produced from the SNO bifurcation, on the other. To prepare for these results, recall that any solution $\hat{T}(t)$ to the DDE (4.4) emanating from $T_{\text{init}}(s)$ ($-\tau \le s \le 0$), is approximated by the solution, $y(t) = (y_1(t), \dots, y_N(t))^T$, to the GK system Eq. (4.10) that emanates from the projection of $T_{\text{init}}$ onto the Koornwinder polynomials $K_n^\tau$ (see (2.8)), according to the formula

$$ (4.20) \qquad \hat{T}(t) \approx \sum_{j=1}^{N} y_j(t), $$

<footnote n="2">²Obtained by approximating $\hat{T}(t)$ solving Eq. (4.4) by the UPOs from the reduced equation according to (4.21) below.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 4. Bifurcation Diagram from the 2D Reduced GK System (4.14). This diagram describes precisely the local and global bifurcations occurring for the (perturbed) Suarez and Schopf model (4.4), given the approximation skills of the reduced system (4.14) in approximating the DDE’s periodic orbits; see Sec. 4.4 below. As for Fig. 5 below, it turned out to be sufficient to use $N = 6$ in the construction of the parameterization $\Psi_{\tau}$ involved in Eq. (4.14), to reach such skills. In each inset, the stable steady states $T_+$ and $T_-$ given by Eq. (4.2) are marked by red dots with $T_+$ corresponding to (0, 0) since this diagram is computed for the perturbed DDE (4.14). This diagram reads as follows. The steady state $T_+$ (resp. $T_-$) is locally stable for all $\tau < \tau_c$, and loses its stability through a subcritical Hopf bifurcation at $\tau_c \approx 1.7408$; see Table 2 and Eq. (4.19). As $\tau$ approaches $\tau^\ddagger \approx 1.5906$ from above, the bifurcating UPOs (black dashed curves) merge into a homoclinic orbit (cyan curve) connecting the stable and unstable directions of the saddle steady state $T_0$, marked by an empty blue circle in inset A. This merging of UPOs terminates the branch of UPOs shown as dashed black line in the diagram. The diminishing of $\tau$ below $\tau^\ddagger$ leads to new UPOs that encompass the homoclinic orbit, with amplitude that grows as $\tau$ approaches $\tau^\ddagger$; see blue dashed curves in inset B. This latter branch of UPOs terminates at $\tau^* \approx 1.5592$ through a SNO bifurcation that gives rise to a branch of stable periodic orbits shown by the red curves in inset C.</figcaption></figure>

see Eq. (2.13) and [CGLW16, Sec. 6].

Given a GK system with its dimension, $N$, sufficiently large so that the approximation error in $L^\infty$ is small in (4.20) (see Table 1), we can expect to reach a good approximation of $\hat{T}(t)$, at least for $\tau$ sufficiently close to $\tau_c$ (due to Theorem 3.1), by replacing the low-mode amplitudes $y_1(t)$ and $y_2(t)$ by their approximations $x_1(t)$ and $x_2(t)$ from the 2D reduced system (4.14), and the $y_n(t)$ by $\Psi_{n,\tau}(x_1(t)e_1 + x_2(t)e_2)$, for $n \ge 3$.

We consider thus the following approximation formula of $\hat{T}(t)$ that is built solely from the solutions of the 2D reduced GK system Eq. (4.14) and the high-mode parameterization $\Psi_{\tau}$ given by (3.27):

$$ (4.21) \qquad \begin{gathered} T^*(t) = \sum_{j=1}^{N} \left( x_1(t)e_1^j + x_2(t)e_2^j + \sum_{n=3}^{N} \Psi_{n,\tau}(x(t))e_n^j \right), \\ \text{with } x(t) = x_1(t)e_1 + x_2(t)e_2 \text{ solving Eq. (4.14),} \end{gathered} $$
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Approximation results: DDE vs Effective Reduced GK system. The solutions are shown in lagged coordinates for $\alpha = 0.75$. For different values of $\tau$ as indicated, the solutions to the DDE (4.4) are compared to those obtained from the formula of $T^*$ given by (4.21), built from the solutions of the 2D reduced GK system Eq. (4.14) and the parameterization $\Psi_{\tau}$ given by (3.27). The stable (resp. unstable) DDE limit cycles are shown in black (resp. blue). The stable (resp. unstable) limit cycles, obtained from the 2D reduced system after lifting through (4.21), are shown by the orange (resp. cyan) dashed curves. The stable (resp. unstable) equilibria are shown as filled (resp. empty) circles. Here, it is sufficient to use a GK system of dimension $N = 6$ to build the parameterization $\Psi_{\tau}$ and the modes and eigenvalues used in (4.21), in order to achieve such approximation skills.</figcaption></figure>

and where $e_{\ell}^j$ denotes the $j^{th}$ component of the eigenmodes $e_{\ell}$ for $1 \le \ell \le N$. Recall that the latter depend on $\tau$ and $N$ as eigenvectors of $A(\tau)$ given by (4.7).

Figure 5 shows, for $N = 6$, the approximation skills by $T^*(t)$ of periodic orbits to the (perturbed) Suarez and Schopf model (4.4) for four different cases. Three of these cases exhibit UPOs coexisting with a stable limit cycle and cover situations for which (i) $\tau$ is close to $\tau_c$ from below ($\tau = 1.7$, Fig. 5A), (ii) $\tau$ is close (from below) to $\tau_c^*$ where the homoclinic orbit emerges ($\tau = 1.6$, Fig. 5B), and (iii) $\tau$ is close to $\tau_c^*$ at which the SNO bifurcation occurs ($\tau = 1.562$, Fig. 5C). The last case shown corresponds $\tau$ away from $\tau_c$ from above ($\tau = 1.9$, Fig. 5D) for which only a stable limit cycle exists as periodic orbit. In each case, the approximation formula (4.21) based on the 2D reduced GK system (4.14) and the parameterization $\Psi_{\tau}$ enables us to achieve highly accurate approximations of the bifurcating DDE periodic orbits, including the corresponding UPOs and stable limit cycles along with their pronounced nonlinear features, that unfold from the subcritical Hopf bifurcation and the SNO bifurcation, respectively; see Fig. 4.

The reason behind these impressive skills achieved by the reduced system (4.14) when $\tau$ varies, lies in a simple property satisfied by the DDE periodic solutions: they sit very close (almost slaved) to the manifold given as graph of $\Psi_{\tau}$ built from the GK system (4.10), already for $N = 6$. This property remains valid even for solution's amplitude (for $\hat{T}(t)$) of order one such as shown for $\tau = 1.9$ in Fig. 6, which corresponds to the furthest value $\tau$ to $\tau_c$ examined in Fig. 5.
These comments are made visual in Fig. 6. This figure shows the graph of the norm of $\Psi_{\tau}$,

$$ (4.22) \qquad \varphi : (z_1, z_2) \mapsto \|\Psi_{\tau}(z)\|_{s}, \quad z = z_1 e_1 + z_2 e_2, $$

as a function of Re($z_1$) and Im($z_1$), where

$$ (4.23) \qquad \|w\|_s = \sqrt{\sum_{j=m_c+1}^{N} |\langle w, e_j^* \rangle|^2}, \quad \text{for any } w \in H_s. $$

The representation of $\varphi$ as a function of Re$z_1$ and Im$z_1$ is made possible since $H_c$ is spanned by eigenvectors that form a complex conjugate pair. The resulting surface, (Re($z_1$), Im($z_1$)) $\mapsto$ $\varphi$(Re($z_1$), Im($z_1$)), is then intercepted by two vertical planes for a better visual inspection of the distance between the stable periodic trajectory from Eq. (4.10) and this surface. The left and right panels of Fig. 6 show the curves (in red) obtained as cross sections with, respectively, the vertical plane corresponding to Re($z_1$) $\equiv$ 0, and that corresponding to Im($z_1$) $= a$ Re($z_1$) + b, with $a$ $=$ $-0.63$ and $b$ $=$ 0.21. The intersections between these planes and the stable periodic solution $\Gamma(t)$ to the DDE (4.4) are shown by the black dots shown in each panel of Fig. 6.

<figure><img src="image_6.png" /><figcaption>**Figure 6.** Visualization of the high-mode parameterization $\Psi_{\tau}$ used in the effective reduced Eq. (4.14). The curves (in red) are obtained as intersections between the graph of $\varphi$ given by (4.22) (for $N$ $=$ 6) with respectively, the vertical plane corresponding to Re$z_1$ $\equiv$ 0 (cross section A), and that corresponding to Im$z_1$ $= a$ Re$z_1$ + b with $a$ $=$ $-0.63$ and $b$ $=$ 0.21 (cross section B). The black dots correspond to the intersections points of the DDE stable limit cycle from (4.4) (for $\tau$ $=$ 1.9), with these vertical planes.</figcaption></figure>

It is noteworthy to emphasize that the good parameterization of the stable periodic solutions to Eq. (4.10) are not limited to the choice of cross sections picked up here. The parameterization skill can be indeed assessed quantitatively by inspecting e.g. the parameterization defect

$$ R(t) = \Psi_{\tau}(\Gamma_1(t), \Gamma_2(t)) - \Gamma_s(t), $$

where $\Gamma_j$ denotes the projection of stable periodic solution $\Gamma(t)$ to the DDE (4.4) onto $e_j$ ($j$ $=$ 1, 2), while $\Gamma_s(t)$ denotes the projection of $P(t)$ onto $H_s$ $=$ span{$e_3, \dots, e_N$}. For the case shown in Fig. 6 ($\tau$ $=$ 1.9), the time-averaged over one period of the ratio, $\|R(t)\|^2/\|\Gamma_s(t)\|^2$, is approximately equal to $2.5 \times 10^{-2}$, indicating that $\Gamma(t)$ is almost slaved to the manifold $\Psi_{\tau}$. This observation explains that the discrepancies observed in Fig. 6 between the DDE solution and the manifold $\Psi_{\tau}$ do not affect the approximation skills of the reduced system (4.14) shown in Fig. 5. Indeed, the fraction of the energy contained in the parameterized modes is so small compared to that contained in the resolved ones, that such discrepancies do not impact the reduced system's approximation skills.
**4.5. Model error estimates.** To complement the numerical results of the previous section, we present now some basic error estimates for the effective reduced GK system (4.12). As shown in Proposition 4.1 below, the parameterization defect associated with the manifold $\Psi_{\tau}$, is the main controlling factor. Although we present the result in the context of the Suarez-Schopf model, the same type of calculations can be performed for GK approximations of DDEs with polynomial nonlinearities. The result is also not limited to reduced systems constructed from the center-unstable manifold parameterizations such as given by (3.17) or (3.27). For this reason, we present the calculation for reduced system based on an arbitrary parameterization function $\varphi_{\tau}: H_c \to H_s$.

Let us first recall that the $N$-dimensional GK system (4.6) of the Suarez-Schopf DDE (4.4) is given by

$$ (4.24) \qquad \frac{d\mathbf{y}}{dt} = A(\tau)\mathbf{y} + G_2(\mathbf{y}) + G_3(\mathbf{y}), $$

with the matrix $A(\tau)$ and the nonlinear terms given by (4.7) and (4.8), respectively. As in Sec. 3, we assume that $A(\tau)$ is diagonalizable over $\mathbb{C}$, which is observed to be true for all the numerical results reported in this study for this DDE model. Recall that $H_c$ and $H_s$ are the eigen subspaces spanned by the low modes and the high modes, respectively.

We will also make use of the following notations. Given a solution $\mathbf{y}(t)$ of (4.24) evolving in $\mathbb{R}^N$, denote by $\mathbf{y}_c = \Pi_c \mathbf{y} = \sum_{j=1}^{m_c} z_j e_j$ the low-mode projection of $\mathbf{y}$ and by $\mathbf{y}_s = \Pi_s \mathbf{y} = \sum_{j=m_c+1}^{N} z_j e_j$ the high-mode projection of $\mathbf{y}$. Note that the vector $\mathbf{z} = (z_1, \dots, z_N)^T$ is nothing else than the vector $\mathbf{y}$ under the eigenbasis of $A(\tau)$. Let $P = [e_1, \dots, e_N]$ be the $N \times N$ matrix whose columns consist of the eigenvectors of $A(\tau)$. We have $\mathbf{y} = P\mathbf{z}$.

We introduce also $P_c$ as the matrix made of the first $m_c$ columns of $P$, namely $P_c = [e_1, \dots, e_{m_c}]$. Similarly, $P_s = [e_{m_c+1}, \dots, e_N]$, $z_c = (z_1, \dots, z_{m_c})^T$ and $z_s = (z_{m_c+1}, \dots, z_N)^T$. We have then $\mathbf{y}_c = P_c \mathbf{z}_c$ and $\mathbf{y}_s = P_s \mathbf{z}_s$. Finally, let $Q$ be the $N \times N$ matrix whose columns consist of the eigenvectors of $A^*$. We have then

$$ (4.25) \qquad Q^*P = I_{N \times N}, $$

where $Q^*$ denotes the complex conjugate transpose. Similarly as for $P_c$, the matrix $Q_c$ denotes the matrix made of the first $m_c$ columns of $Q$.

**Proposition 4.1.** Consider $\mathbf{y}(t)$ to be a solution to Eq. (4.24) over a given interval $[0, t_m]$. Assume that $A(\tau)$ is diagonalizable over $\mathbb{C}$. Let $\varphi_{\tau}: H_c \to H_s$ be a parameterization of $\mathbf{y}_s$ in terms of $\mathbf{y}_c$. Assume that $\mathbf{y}(t)$ satisfies

$$ (4.26) \qquad \|\mathbf{y}(t)\| \le R, \quad \|\mathbf{y}_c(t) + \varphi_{\tau}(\mathbf{y}_c(t))\| \le R, \quad t \in [0, t_m]. $$

Then there exists a constant $C(R) > 0$ such that

$$ (4.27) \qquad \left\| \dot{\mathbf{z}}_c - \Lambda_c \mathbf{z}_c - Q_c^* G_2 \left( P_c \mathbf{z}_c + \varphi_{\tau} (P_c \mathbf{z}_c) \right) \right\|^2 \le C(R) \overline{\left\| P_s \mathbf{z}_s - \varphi_{\tau} (P_s \mathbf{z}_s) \right\|^2}, $$

where $\overline{(\cdot)}$ denotes the time average over $(0, t_m)$ and $\Lambda_c$ denotes the $m_c \times m_c$ diagonal matrix made of the eigenvalues $\lambda_1, \dots, \lambda_{m_c}$ of $A(\tau)$.

*Proof.* The desired result (4.27) follows from a direct calculation by comparing the vector field of the reduced system associated with $\varphi_{\tau}$ (see (4.30) below) and the projection of the GK vector field in (4.24) onto $H_c$. Indeed, note that the GK system (4.24) can be rewritten as:

$$ (4.28) \qquad \dot{\mathbf{z}} = \Lambda(\tau)\mathbf{z} + Q^*(G_2(P\mathbf{z}) + G_3(P\mathbf{z})), $$

where $\Lambda(\tau) = Q^*A(\tau)P$.

By projecting the above equation onto the low modes, we have

$$ (4.29) \qquad \dot{\mathbf{z}}_c = \Lambda_c(\tau)\mathbf{z}_c + Q_c^*(G_2(P_c\mathbf{w}) + G_3(P_c\mathbf{w})). $$

On the other hand, the reduced system associated with $\varphi_{\tau}$ is

$$ (4.30) \qquad \dot{\mathbf{w}} = \Lambda_c(\tau)\mathbf{w} + Q_c^*(G_2(P_c\mathbf{w}) + \varphi_{\tau}(P_c\mathbf{w})) + G_3(P_c\mathbf{w} + \varphi_{\tau}(P_c\mathbf{w})). $$
We define the model error as

$$
\begin{align*}
R_T(\mathbf{z}_c, \varphi_\tau) &= \int_0^{t_m} \left\| \frac{\mathrm{d}\mathbf{z}_c}{\mathrm{d}t} - \Lambda_c \mathbf{z}_c - Q_c^* \left( G_2(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) + G_3(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) \right) \right\|^2 \mathrm{d}t \\
(4.31) \qquad &= \int_0^{t_m} \left\| Q_c^* \left( G_2(\mathbf{y}) + G_3(\mathbf{y}) - G_2(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) - G_3(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) \right) \right\|^2 \mathrm{d}t \\
&\le C \int_0^{t_m} \left\| G_2(\mathbf{y}) + G_3(\mathbf{y}) - G_2(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) - G_3(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) \right\|^2 \mathrm{d}t.
\end{align*}
$$

Let us introduce the notation

$$
(4.32) \qquad \mathbf{w}_5 = \varphi_{\tau}(P_c \mathbf{z}_c),
$$

and recall also that $\mathbf{y}_c = P_c \mathbf{z}_c$ and $\mathbf{y}_5 = P_5 \mathbf{z}_5$. Then, by the definitions of $G_2$ and $G_3$ given in (4.8), we have

$$
\begin{align}
G_2(\mathbf{y}) - G_2(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) &= -3T_+ \left( \left( \sum_{n=1}^N y_{c,n} + \sum_{n=1}^N y_{5,n} \right)^2 - \left( \sum_{n=1}^N y_{c,n} + \sum_{n=1}^N w_{5,n} \right)^2 \right) \boldsymbol{\nu}_N, \tag{4.33} \\
&= -3T_+ \left( \sum_{n=1}^N (2y_{c,n} + y_{5,n} + w_{5,n}) \right) \left( \sum_{n=1}^N (y_{5,n} - w_{5,n}) \right) \boldsymbol{\nu}_N, \notag
\end{align}
$$

and

$$
\begin{align}
G_3(\mathbf{y}) - G_3(P_c \mathbf{z}_c + \varphi_\tau(P_c \mathbf{z}_c)) &= -\left(\left(\sum_{n=1}^N y_{c,n} + \sum_{n=1}^N y_{5,n}\right)^3 - \left(\sum_{n=1}^N y_{c,n} + \sum_{n=1}^N w_{5,n}\right)^3\right) \boldsymbol{\nu}_N \tag{4.34} \\
&= -\mathcal{N}(\mathbf{y}, \mathbf{w}_5) \left(\sum_{n=1}^N (y_{5,n} - w_{5,n})\right) \boldsymbol{\nu}_N, \notag
\end{align}
$$

where

$$
\mathcal{N}(\mathbf{y}, \mathbf{w}_5) = 3 \left( \sum_{n=1}^{N} y_{c,n} \right)^2 + \left( \sum_{n=1}^{N} y_{5,n} \right)^2 + \left( \sum_{n=1}^{N} w_{5,n} \right)^2 \\
+ 3 \left( \sum_{n=1}^{N} y_{c,n} \right) \left( \sum_{n=1}^{N} (y_{5,n} + w_{5,n}) \right) + \left( \sum_{n=1}^{N} y_{5,n} \right) \left( \sum_{n=1}^{N} w_{5,n} \right).
$$

Finally, by using (4.33) and (4.34) in (4.31), we conclude, in virtue of assumption (4.26), that

$$
R_T(\mathbf{z}_c, \varphi_\tau) \le C(R) \int_0^{t_m} \left| \sum_{n=1}^N (y_{5,n} - w_{5,n}) \right|^2 \mathrm{d}t \le C(R) \int_0^{t_m} \| \mathbf{y}_5(t) - \mathbf{w}_5(t) \|^2 \mathrm{d}t. \quad \square
$$

5. TIPPING SOLUTION PATHS AND ENSO VARIABILITY

5.1. Transition paths and nonlinear building blocks of temporal variability. The plausibility of a stochastic forcing as a mechanism for ENSO irregularity has been argued in many studies. Physical origins of such a forcing include large-scale synoptic atmospheric transients such as the Madden Julian Oscillation [BH05] or westerly wind bursts [Fed02]. Dynamically, the idea is to explicitly separate the slow and fast modes in the atmosphere and add the latter to simple deterministic models as a stochastic forcing term. Other sources of stochasticity include processes associated with atmospheric/moist convective disturbances whose timescales can vary from hours to weeks. Typically, the additional fast-mode random forcing disrupts the slow scales and convert the original periodic or damped oscillation supported by the deterministic model into an irregular one. Such mechanisms have been previously advocated through a combination of observational data and a
hierarchy of ENSO models including data-driven models [PS95, KKRG05, CKG11, CCH$^{+}$16]; PDE models [BNG97, CMT18, EL97, TMCS16, RN00, ZGMPK03], and conceptual models [CFY22, CZ23].

Here, we propose a novel scenario for the fabric of ENSO variability. Our approach is based on (i) the dynamical insights gained from the bifurcation diagram and phase portrait of the deterministic dynamics as summarised by Fig. 4, and (ii) the design of stochastic model’s disturbances interacting with the underlying nonlinear invariant sets (UPOs, stable limit cycles, and homoclinic orbit) and the topology they form across a $\tau$-interval $[\tau_0, \tau_1]$ containing $[\tau^*, \tau^\#]$, where $\tau^*$ (resp. $\tau^\#$) corresponds to the SNO (resp. homoclinic) bifurcation point. Since these invariant sets are $\tau$-dependent, as dictated by the bifurcation diagram of Fig. 4, we naturally consider these rapidly varying disturbances to be superimposed to slowly varying lagged effects. Physically, such slow variations can be envisioned as resulting from slight variations occurring in key properties of the equatorially trapped oceanic waves such as their transit time of the Pacific ocean.

In that respect, we propose the following stochastic model that we write in the perturbed variable $\theta(t) = T - T_+$:

$$ (5.1a) \qquad d\theta(t) = (a\theta(t) - \alpha\theta(t - \tau(t)) - b\theta^2(t) - \theta^3(t)) dt + \frac{\sigma}{1 + |\theta^2(t)|} dW_t, $$

$$ (5.1b) \qquad \text{with } \tau(t) = \tau_0 + \epsilon t, \quad \tau(t) \in [\tau_0, \tau_1], $$

where $a = 1 - 3T_+^2$, $b = 3T_+$, and $\epsilon, \sigma \ge 0$ while $W_t$ denotes a Brownian motion. The noise term is a Lorentzian function favouring larger noise when $\theta \approx 0$ while damping the noise effects when $\theta$ gets too large. It is interpreted in the sense of Itô to fix ideas [Gar04, Chapter 4]. Due to this noise term, a stochastic path solving (5.1) is more likely to meander away from $T_+$ than around it, as time flows. Since $T_+$ corresponds for the Suarez and Schopf model to the basic steady state associated with an El Niño event, using such a nonlinear noise favours, intuitively, a less frequent occurrence of metastable El Niño events corresponding to visiting a neighborhood of $T_+$ than if Eq. (5.1) would be driven by a pure white noise.<fnref n="3" /> As discussed below, Eq. (5.1) supports actually the occurrence of other types of El Niño events that do not correspond to a metastable visit of $T_+$ but rather to experiencing excursions across the bifurcating solutions of the deterministic equation as identified in Sec. 4.3.

To show evidence of these other types of El Niño events, we performed an integration of Eq. (5.1) for $\alpha = 0.75$, $\sigma = 0.2$, with $\tau$ evolving according to Eq. (5.1b) over an interval $[0, t_m]$, with $\tau_0 = 1.45$, $\epsilon = 8.4 \times 10^{-4}$, and $t_m = 237.8$ in time unit of the model. With these parameters, we have $\tau(t_m) = \tau_1 = 1.65$, and thus the interval $[\tau_0, \tau_1]$, within which $\tau(t)$ varies, contains both the SNO critical value $\tau^*$ and the homoclinic critical value $\tau^\#$; see caption of Fig. 4 for numerical values of $\tau^*$ and $\tau^\#$. The rationale behind this numerical setup is guided by a simple intuition. Indeed, as the delay parameter $\tau$ slowly drifts across a portion of the bifurcation diagram shown in Fig. 4 and when the noise strength parameter $\sigma$ is not too small, it is expected that in the course of integration of Eq. (5.1), its solution would “hop” through the different branches of the bifurcation diagram and thus reveal fingerprints of the various invariant sets (UPOs, stable limit cycles, steady states) of the deterministic flow. We call such a solution a Tipping Solution Path (TSP).

In the literature, different tipping mechanisms have been identified that include bifurcation-induced tipping, noise-induced tipping, and rate-induced tipping [AWVC12, FPS18, Kue11]. It is known that under different setups of the external forcing and/or parameter drifting scenarios, a given model may experience different types of tipping; see e.g. [AWVC12, Section 4]. A mathematical characterization of the different types of tipping phenomenon in terms of the model parameters in Eq. (5.1) and drifting scenarios, requires a separate study. Instead, our focus is on the impact on the solution’s variability of the presence of tipping points involving UPOs and in that regard, the parameter setting and drifting scenario (5.1b) chosen above have been found to be physically insightful. In particular, as explained below (see Figs. 7 and 8), we observed that even during time

<footnote n="3">³The nonlinear noise favours also a less frequent occurrence of El Niño events with unrealistic long duration visits that would consist of e.g. meandering near $T_+$ for a too long amount of time.</footnote>
intervals preceding the crossing of the SNO critical value $\tau^*$, the solution's variability is affected by invariant sets such as limit cycles that emerge only for $\tau > \tau^*$.

We numerically observed that this ability of the stochastic model to “feel” the deterministic model’s dynamics at a given $\tau$-value prior $\tau(t)$ crosses this value, is actually independent of how the noise is interpreted (Itô vs. Stratonovich).<fnref n="4" /> Thus, noise can e.g. excites parameter regimes with oscillatory dynamics as $\tau(t)$ approaches $\tau$-values supporting oscillations for the unperturbed model. In terms of solution’s variability, this can be reflected on trajectory segments of the stochastic solution exhibiting then “fingerprint” of these underlying nonlinear oscillations.

Thus, noise-sustained oscillations associated with the SNO critical transition play an important role in the fabric of the temporal variability of a TSP solving Eq. (5.1). Their role is not exclusive though as other nonlinear invariant sets (steady states, homoclinic orbits) have also their share in the TSP’s variability. In that respect, we discuss below how the whole set of these nonlinear building blocks allow us to break down ENSO-like events from the model (5.1).

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Tipping solution path and phase portrait visualization of two of its transition snippets. One transition snippet emanates near the La Niña steady state $T_-$ (of negative value) and spends some time meandering around it (turquoise path in panel A) after escaping it via a nearly cyclic noisy trajectory (blue path in panel A). The corresponding temporal pattern is shown on $\theta(t)$ in panel B with the same color coding. The other transition snippet emanates from a location marked by the empty square in panel A and meanders around the UPO surrounding $T_-$ before escaping via a nearly cyclic noisy trajectory (red path in panel A). The corresponding episode is shown in the time domain with the same color coding in panel B. The two vertical green dashed lines in panel B mark the time-instants at which $\tau(t) = \tau^*$ ($t = 130$) and $\tau(t) = \tau^#$ ($t = 167.38$). In panel A, besides the stochastic trajectory, are shown several invariant sets (UPOs, stable limit cycles, and homoclinic orbit) for the deterministic model, for different $\tau$-values as pointed out in the text.</figcaption></figure>

Such a TSP is shown in Fig. 7B. To understand the dynamical origins behind the time-variability displayed by this TSP, we focus on a few episodes (snippets) as marked in color in Fig. 7B that we map, with the same color coding, to the phase portrait shown in Fig. 7A. To help interpretation, these snippets are superimposed on a few invariant curves of the deterministic model (4.4)<fnref n="5" /> encountered as $\tau$ evolves from $\tau_0$ to $\tau_1$ according

<footnote n="4">⁴In the Stratonovich sense ([Gar04, Section 4.3.6]), an additional term $\sigma^2\theta(t)(1 + \theta^2(t))^{-3}dt$ is added to Eq. (5.1a). In the present parameter setting, we have checked that this extra term has little influence on the dynamics as it scales here with $\sigma^2$.</footnote>

<footnote n="5">⁵as calculated from our 2D reduced equation; see Secns. 4.3 and 4.4.</footnote>
to (5.1b): a UPO at $\tau = 1.5607$ encompassing the homoclinic orbit, two UPOs symmetric with respect to the saddle point (for $\tau = \tau_1 = 1.65$), and two stable limit cycles, shown as two closed black (plain) curves whose one with smaller diameter corresponds to $\tau = 1.65$ while the other corresponds to $\tau = 2$.

The snippets marked by the turquoise and blue consecutive segments of the TSP as shown in Fig. 7B, correspond to a quiet episode (turquoise) followed by a large excursion (blue) reminiscent with what has been observed for certain El Niño events as recorded by the Niño 3.4 SST index. Such a large excursion following a certain stillness evokes the major El Niño event that occurred in Dec-Jan 2016 [LTW$^+$17]. It is worthwhile noting that the Niño 3.4 SST index during the four years preceding this event exhibited indeed temporal patterns resembling the episode marked in turquoise; cf. e.g. the global reanalysis product of Niño 3.4 index as documented by the E.U. Copernicus Marine Service Information [CME18]. The phase portrait reveals that these consecutive events correspond, in our model, to a meandering of the solution path near what corresponds to a La Niña steady state $T_-$ (turquoise path in Fig. 7A)$^6$ followed by a nearly cyclic orbit (blue path in Fig. 7A).

The snippets shown in yellow and red in Fig. 7B emphasize a different mechanism leading to an El Niño event. Indeed, by inspecting the corresponding path in Fig. 7A, we observe that while the El Niño excursion corresponds still to a nearly cyclic orbit (red), its ignition is now preceded by a path's wandering around an UPO surrounding itself the La Niña steady state. Due to their behavior starting either from the La Niña steady state or an UPO surrounding it, before engaging into a quasi cyclic excursion, we call such episodes as Transition Snippets.

Such transition snippets provide the generic building blocks at work for the production of El Niño-like patterns displayed by the solution to Eq. (5.1) such as shown in Fig. 7B. Our analysis reveals thus that a solution path wandering around an UPO surrounding the La Niña steady state or having a closer meandering around it can be interpreted as an early warning signal [SBB$^+$09] for an El Niño-like pattern to occur in this model.

Other episodes displayed by the TSP's temporal patterns correspond to a path connecting the El Niño with the La Niña equilibria. Such an event is shown in Fig. 8 with color coding marking the following three stages: (i) path escape towards the El Niño equilibrium (blue path), (ii) path meandering in a neighborhood of that equilibrium (yellow path), and (iii) path escape from that neighborhood towards the La Niña equilibrium (red path).

<figure><img src="image_8.png" /><figcaption>**Figure 8.** A transition snippet going through a neighborhood of $T_+$ before landing near $T_-$. The empty square in panel A marks the beginning of the blue path. See Text.</figcaption></figure>

<footnote>$^6$To avoid multiple notations, we still denote by $T_-$ the steady state to Eq. (4.4) that corresponds to $T_-$ given by (4.2).</footnote>
The mechanisms of fabric of the TSP’s temporal patterns reported here are therefore diverse, involving various scenarios such as noise-sustained oscillations associated with the SNO critical transition and its constitutive UPOs and stable limit cycles, or transitions between steady states. Due to the presence of (nonlinear) periodic orbits these mechanisms are characteristic of more intricate interactions between noise and nonlinear effects than classically encountered in noise-induced tipping dynamics between multiple equilibria; see e.g. [GHJM98] and [HL84, Chapter 6]. In particular, the good capture of such interactions by stochastic reduced models remains to be explored. We mention that stochastic invariant manifolds techniques [CLW15b, CLMW23] and their extension [CLM23] should be useful in that respect.

Finally, it worth mentioning that such mechanisms are not dependent on the noise path used to drive Eq. (5.1a) but rather conditioned to the choice of $τ_0$ in Eq. (5.1b) and $σ$ in Eq. (5.1a). Choosing $τ_0$ too far below $τ^*$ and/or $σ$ too small could lead to completely different solution’s temporal behavior with for instance more frequent transitions between equilibria and very rare excitation of oscillatory behavior. While allowing $τ$ in Eq. (5.1) to drift is motivated by physical considerations as mentioned above, it is not excluded though, for well-calibrated model’s parameters, that stochastic solutions for frozen $τ$-value could share similar temporal variability attributes than those reported here. Regardless, drifting or not, the nonlinear invariant sets (UPOs, stable limit cycles, steady states) of the deterministic flow play a central role in shaping the temporal patterns exhibited by the TSP. We have been though voluntarily evasive on the role of the homoclinic orbit in the fabric of the TSP variability. The next section clarifies this point.

## 5.2. Decadal variability and homoclinic orbit.
To decipher the role of the homoclinic orbit in the fabric of the TSP variability, we need to quantify this variability in more precise terms. To do so, we analyze the spectral content of a TSP produced over a long-time integration of Eq. (5.1) consisting of $10^6$ time steps with $δt = 2 \times 10^{-3}$, in which $τ(t)$ is now allowed to oscillate in a periodic fashion between $τ_0$ and $τ_1$ in (5.1b). The values of these and other model’s parameters are the same as in Sec. 5.1. To make a physical sense of the variability displayed by the TSP, we adopt the following conversion formula [BTR07] $s = tΔ/τ$ with $Δ = 349$, and $τ = 1.7$, in which $s$ denotes the physical time in year, while $t$ denotes the non-dimensional time unit of the DDE model. To analyse the frequency content, we use the data-adaptive harmonic decomposition (DAHD) method [CK17] that is a signal processing method which has been employed in many applications to analyze the temporal variability and extract coherent patterns of complex (multivariate) time series; see e.g. [KCG18, KCYG18, KC18, KCB18]. Said in simple terms, the method allows for performing a Fourier-like analysis from time-lagged correlations while extracting empirical modes of variability that are naturally ranked per Fourier frequency [CK17, Theorem V.1]. A simple projection of the signal onto a selected group of such modes over a given frequency range allows in turn for calculating in the time domain the contribution of that range within the original signal.

We applied DAHD to the aforementioned TSP simulated over a time window of approximately 1, 100 yrs. The DAH power spectrum reveals two distinct broadband peaks: a 21-yr peak associated with decadal variability, and a 5.5-yr peak associated with ENSO interannual variability, corresponding in Fig. 9A, to the small bump made of blue dots for decadal variability, and to the more energetic one located to its right, for the interannual one. The projection of the TSP onto the subspace spanned by the DAH modes (DAHMs) associated with the decadal variability (blue dots in Fig. 9A) is shown by the blue curve in Fig. 9B.

Decadal variability of ENSO has been documented in recent studies [DCP$^+$21]. The origin of this decadal variability in our model (5.1) can be traced back to UPOs that are located close to the homoclinic orbit, i.e. for $τ$ close to $τ^*$ from below; see blue curves in inset B of Fig. 4. Recall that the 2D reduced GK system (4.14) allows for a simple computation of the UPOs of the DDE (4.4) by a simple backward integration of Eq. (4.14) followed by a lifting of the orbit via the formula (4.21). Using this approach, UPOs closer to the homoclinic orbit such as shown in Fig. 4B can be easily probed. Among these UPOs, we found that the UPO with periodicity of 18.23-yr (blue curve in Fig. 10-lower panel, corresponding to $τ = 1.5827$ in the deterministic model) provides a good synchronization with the temporal (modulated) patterns exhibited by the decadal variability extracted from the
<figure><img src="image_9.png" /><figcaption>**Figure 9.** DAH power spectrum and decadal variability. The decadal peak located at ≈ 21 yr is marked by the vertical grey dashed line in panel A. The ENSO peak is located around ≈ 5.5 yr. The projection of the TSP onto the subspace spanned by the DAH modes (DAHMs) associated with the decadal variability (blue dots in panel A) is shown by the blue curve in panel B.</figcaption></figure>

DAHD (black curve in Fig. 10-lower panel). Thus, one can infer that the homoclinic bifurcation occurring at $\tau = \tau^{\ddagger}$ is responsible for the emergence of UPOs that provide the backbone of the decadal variability exhibited by our TSP. As shown in the top panel of Fig. 10, the stable limit cycles are responsible for the more regular oscillatory events occurring on the TSP on an interannual timescale.

<figure><img src="image_10.png" /><figcaption>**Figure 10.** Interannual and Decadal variability. The black curve shown in the lower panel corresponds to the decadal variability time series shown in blue in the right panel of Fig. 9. The blue curve shown here corresponds to a UPO of period 18.23-yr obtained for $\tau = 1.5827$, which lies close to the homoclinic orbit. The upper panel shows the TSP in black and a stable limit cycle in red of period 5.78-yr obtained for $\tau = 1.7689$.</figcaption></figure>

## 6. CONCLUDING REMARKS

Thus, our approach based on GK approximations and higher-order approximations to center-unstable manifolds to analyze bifurcations in the Suarez and Schopf model allowed us to exhibit the nonlinear building blocks (UPOs, stable limit cycles, homoclinic orbits, steady states) at play in the fabric of a rich temporal variability exhibited by solution paths drifting through the underlying bifurcations in presence of noise. The drift corresponds to slow changes in e.g. the traveling time of equatorially trapped waves (i.e. slow drift in $\tau$ through its critical values) while the noise term in (5.1) can be interpreted as a nonlinear coupling term between ocean and atmosphere.
This new mechanism identified in this study to produce interannual and decadal variability of ENSO-like patterns deserves more explorations, in particular in ENSO models including more physics while still tractable mathematically, such as the Cane-Zebiak models and the like; see [JN93, NJ93, CCHT19, TCDN20]. In that respect, exploring the effects of distributed delays in the Suarez and Schopf could be insightful as suggested in [FQS+19]. The approach presented in this work applies to such models where the integral terms of the form $\int_{t-\tau}^{t} x(s) ds$ are simply approximated by $\tau y_0(t) - \tau \sum_{n=1}^{N-1} y_n(t)$, in the corresponding GK approximations; see (2.17). Similarly the inclusion of multiple delays does not constitute a limitation to the approach presented here [CGLW16, CKL20], nor the approach is limited to the used of center-unstable manifolds formulas. In that respect, more general parameterizations of the disregarded modes such as introduced in [CLM20] could be relevant to derive effective reduced GK systems able to handle more complex bifurcations [CLM23].

Questions about the role of the annual cycle on this new route to complexity, in particular regarding its interactions with the underlying UPOs is a topic that should also lead to instructive dynamical insights in the tradition of the so-called slow-mode chaos approach to ENSO variability [CWLJ94, JNG94, TSCJ94, JNG96, GT00].

Finally, we would like to mention that strong evidences about the prominent role of UPOs in structuring atmospheric events, corresponding to zonal and blocking events in a low resolution quasi-geostrophic model, have been documented in [LG20]. Our study pointed out the subtle interactions that UPOs may have with stochastic disturbances to produce a wealth of temporal events. We hope that more studies about the role of UPOs in the understanding of climate variability will be pursued in the future.

## ACKNOWLEDGEMENTS

The authors thank Tetsushi Saburi for constructive comments on the GK method. The authors also express their gratitude to the anonymous referees for their insightful comments that helped improve the presentation and discussion of the results. This work has been partially supported by the Office of Naval Research (ONR) Multidisciplinary University Research Initiative (MURI) grant N00014-20-1-2023, by the National Science Foundation grant DMS-2108856, and by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 810370).

## APPENDIX A. DETERMINATION OF $P_N$ AND $Q_N$ IN Eq. (4.10)

The calculations of $P_N$ and $Q_N$ are based on the RHS of (4.7). $Q_N$ depends on the model parameters, while $P_N$ is independent of those as explained below.

First, the terms

$$ \frac{1}{\|\mathcal{K}_j\|_{\mathcal{E}}^2} \sum_{n=0}^{N-1} (1 - 3T_+^2 - \alpha K_n(-1)) $$

when collected for each $j$ and $n$ allow us to build up $Q_N$ readily once one recalls (2.6) and that due to (2.2),

$$ (A.1) \qquad K_n(-1) = (n^2 + n + 1)(-1)^n, $$

since $L_n(-1) = (-1)^n$ from the properties of the Legendre polynomials.

The determination of $P_N$ which results from the collection of the terms $\sum_{k=0}^{n-1} a_{n,k} (\delta_{j,k} \|\mathcal{K}_j\|_{\mathcal{E}}^2 - 1)$, requires the knowledge of the coefficients $a_{n,k}$. These coefficients arise in the expression of the derivatives of the Koornwinder polynomials in terms of these polynomials themselves; see [CGLW16, Appendix B] for more details.

The coefficients $a_{n,k}$ are then obtained by simply solving the triangular system (A.3) described in the following proposition.
**Proposition 1.** The Koornwinder polynomial $K_n$ of degree $n$ defined in (2.2) satisfies the differential relation

$$ (A.2) \qquad \frac{dK_n}{ds}(s) = \sum_{k=0}^{n-1} a_{n,k} K_k(s), \quad s \in (-1, 1), $$

where the vector $\mathbf{a}_n$ made of the $a_{n,k}$-coefficients,

$$ \mathbf{a}_n = (a_{n,0}, \dots, a_{n,n-1})^T, $$

solves the following upper triangular system

$$ (A.3) \qquad \mathbf{T} \mathbf{a}_n = \mathbf{b}_n, $$

with $\mathbf{T} = (\mathbf{T}_{i,j})_{n \times n}$ and $\mathbf{b}_n = (b_{n,0}, \dots, b_{n,n-1})^T$ given by

$$ (A.4) \qquad \mathbf{T}_{i,j} = \begin{cases} 0, & \text{if } j < i, \\ i^2 + 1, & \text{if } j = i, \\ -(2i + 1), & \text{if } j > i, \end{cases} \quad \text{where } 0 \le i, j \le n-1, \\ b_{n,i} = \begin{cases} -\frac{1}{2}(2i+1)(n+i+1)(n-i), & \text{if } n+i \text{ is even,} \\ (n^2+n)(2i+1) - \frac{i}{2}(n+i)(n-i+1), \\ -\frac{1}{2}(i+1)(n-i-1)(n+i+2), & \text{if } n+i \text{ is odd.} \end{cases} $$

Finally, the rescaled Koornwinder polynomials satisfy

$$ (A.5) \qquad \frac{dK_n^\tau}{d\theta}(\theta) = \frac{2}{\tau} \sum_{k=0}^{n-1} a_{n,k} K_k^\tau(\theta), \quad \theta \in (-\tau, 0). $$

**Remark 1.** We note that the formula for $b_{n,i}$ that appeared in [CGLW16] and in [CKL20] contains two typos (a sign error and a factor 1/2 missing) that are here rectified in the formula (A.4). We would like to express our gratitude to Tetsushi Saburi for pointing out these typos to us. We emphasize though that it is only a typo as the correct formulas were implemented in our codes, thus not questioning the diverse numerical results published in [CGLW16], [CKL18], [CKL20], and [CKLL22].

REFERENCES

[ABL$^{+22}$] A. Andò, D. Breda, D. Liessi, S. Maset, F. Scarabel, and R. Vermiglio, 15 years or so of pseudospectral collocation methods for stability and bifurcation of delay equations, Accounting for Constraints in Delay Systems, Springer, 2022, pp. 127–149.

[ACN16] P. Ashwin, S. Coombes, and R. Nicks, Mathematical frameworks for oscillatory network dynamics in neuroscience, Journal of Mathematical Neuroscience 6 (2016), no. 2, 1–92.

[AR23] M. Anikushin and A. Romanov, Hidden and unstable periodic orbits as a result of homoclinic bifurcations in the Suarez–Schopf delayed oscillator and the irregularity of ENSO, Physica D: Nonlinear Phenomena 445 (2023), 133653.

[Aut] AUTO, https://github.com/auto-07p/, Accessed: 2023-09-09.

[AWVC12] P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox, Tipping points in open systems: Bifurcation, noise-induced and rate-dependent examples in the climate system, Phil. Trans. Roy. Soc. A 370 (2012), no. 1962, 1166–1184.

[BB78] H. T. Banks and J. A. Burns, Hereditary control problems: Numerical methods based on averaging approximations, SIAM J. Control Optim. 16 (1978), 169–208.

[BCL$^{+17}$] N. Boers, M. D. Chekroun, H. Liu, D. Kondrashov, D.-D. Rousseau, A. Svensson, M. Bigler, and M. Ghil, Inverse stochastic–dynamic models for high-resolution Greenland ice core records, Earth Syst. Dynam. 8 (2017), 1171–1190.

[BDG$^{+16}$] D. Breda, O. Diekmann, M. Gyllenberg, F. Scarabel, and R. Vermiglio, Pseudospectral discretization of nonlinear delay equations: new prospects for numerical bifurcation analysis, SIAM Journal on Applied Dynamical Systems 15 (2016), 1–23.

[BH89] D. S. Battisti and A. C. Hirst, Interannual variability in a tropical atmosphere-ocean model: Influence of the basic state, ocean geometry and nonlinearity, J. Atmos. Sci. 46 (1989), 1687–1712.
[BH05] C. Batstone and H.H. Hendon, Characteristics of stochastic variability associated with ENSO and the role of the MJO, Journal of Climate 18 (2005), no. 11, 1773-1789.

[BK79] H. T. Banks and F. Kappel, Spline approximations for functional differential equations, Journal of Differential Equations 34 (1979), no. 3, 496-522.

[BK98] W.-J. Beyn and W. Kleß, Numerical Taylor expansions of invariant manifolds in large dynamical systems, Numer. Math. 80 (1998), no. 1, 1-38.

[BMV05] D. Breda, S. Maset, and R. Vermiglio, Pseudospectral differencing methods for characteristic roots of delay differential equations, SIAM Journal on Scientific Computing 27 (2005), 482-495.

[BNG97] B. Blanke, J.D. Neelin, and D. Gutzler, Estimating the effect of stochastic wind stress forcing on ENSO irregularity, Journal of Climate 10 (1997), no. 7, 1473-1486.

[BRI84] H. T. Banks, I. G. Rosen, and K. Ito, A spline based technique for computing Riccati operators and feedback controls in regulator problems for delay equations, SIAM J. Sci. Stat. Comput. 5 (1984), 830-855.

[BTR07] I. Boutle, R.H.S. Taylor, and R.A. Römer, El Niño and the delayed action oscillator, American Journal of Physics 75 (2007), 15-24.

[BVV14] D. Breda and E. Van Vleck, Approximating Lyapunov exponents and Sacker-Sell spectrum for retarded functional differential equations, Numerische Mathematik 126 (2014), no. 2, 225-257.

[BZ13] A. Bellen and M. Zennaro, Numerical Methods for Delay Differential Equations, Numerical Mathematics and Scientific Computation, Oxford University Press, 2013.

[CCH$^{+}$16] C. Chen, M.A. Cane, N. Henderson, D. E. Lee, D. Chapman, D. Kondrashov, and M. D. Chekroun, Diversity, nonlinearity, seasonality, and memory effect in ENSO simulation and prediction using empirical model reduction, Journal of Climate 29 (2016), no. 5, 1809-1830.

[CCHT19] Y. Cao, M. D. Chekroun, A. Huang, and R. Temam, Mathematical Analysis of the Jin-Neelin Model of El Niño-Southern-Oscillation, Chinese Annals of Mathematics, Series B 40 (2019), no. 1, 1-38.

[CFY22] N. Chen, X. Fang, and J.-Y. Yu, A multiscale model for El Niño complexity, npj Climate and Atmospheric Science 5 (2022), 16.

[CGLW16] M. D. Chekroun, M. Ghil, H. Liu, and S. Wang, Low-dimensional Galerkin approximations of nonlinear delay differential equations, Disc. Cont. Dyn. Sys. A 36 (2016), no. 8, 4133-4177.

[CGN18] M.D. Chekroun, M. Ghil, and J. D. Neelin, Pullback attractor crisis in a delay differential ENSO model, Advances in Nonlinear Geosciences (A. Tsonis, ed.), Springer, 2018, pp. 1-33.

[CK17] M. D. Chekroun and D. Kondrashov, Data-adaptive harmonic spectra and multilayer Stuart-Landau models, Chaos 27 (2017), no. 9, 093110.

[CKG11] M. D. Chekroun, D. Kondrashov, and M. Ghil, Predicting stochastic systems by noise sampling, and application to the El Niño-Southern Oscillation, Proc. Natl. Acad. Sci. 108 (2011), 11766-11771.

[CKL18] M. D. Chekroun, A. Kröner, and H. Liu, Galerkin approximations for the optimal control of nonlinear delay differential equations, Hamilton-Jacobi-Bellman Equations. Numerical Methods and Applications in Optimal Control. D. Kalise, K. Kunisch, and Z. Rao (Eds.), vol. 21, Berlin, Boston: De Gruyter, 2018, pp. 61-96.

[CKL20] M. D. Chekroun, I. Koren, and H. Liu, Efficient reduction for diagnosing Hopf bifurcation in delay differential systems: Applications to cloud-rain models, Chaos 40 (2020), no. 8, 053130, doi:10.1063/5.0004697.

[CKLL22] M. D. Chekroun, I. Koren, H. Liu, and H. Liu, Generic generation of noise-driven chaos in stochastic time delay systems: Bridging the gap with high-end simulations, Science Advances 8 (2022), no. 46, eabq7137, doi.org/10.1126/sciadv.abq7137.

[CLM20] M. D. Chekroun, H. Liu, and J. C. McWilliams, Variational approach to closure of nonlinear dynamical systems: Autonomous case, Journal of Statistical Physics 179 (2020), 1073-1160.

[CLM23] M.D. Chekroun, H. Liu, and J.C. McWilliams, Optimal parameterizing manifolds for anticipating tipping points and higher-order critical transitions, Chaos 33 (2023), 093126.

[CLMW23] M. D. Chekroun, H. Liu, J. C. McWilliams, and S. Wang, Transitions in stochastic non-equilibrium systems: Efficient reduction and analysis, Journal of Differential Equations 346 (2023), 145-204, doi.org/10.1016/j.jde.2022.11.025.

[CLW15a] M. D. Chekroun, H. Liu, and S. Wang, Approximation of Stochastic Invariant Manifolds: Stochastic Manifolds for Nonlinear SPDEs I, Springer Briefs in Mathematics, Springer, New York, 2015.

[CLW15b] ______, Stochastic Parameterizing Manifolds and Non-Markovian Reduced Equations: Stochastic Manifolds for Nonlinear SPDEs II, Springer Briefs in Mathematics, Springer, New York, 2015.

[CME18] CMEMS, Niño 3.4 Sea Surface Temperature time series from Reanalysis, E.U. Copernicus Marine Service Information. Marine Data Store (MDS) (2018), doi:10.48670/moi-00219.

[CMT18] N. Chen, A. J. Majda, and S. Thual, Observations and mechanisms of a simple stochastic dynamical model capturing El Niño diversity, Journal of Climate 31 (2018), 449-471.

[Cra91] J. D. Crawford, Introduction to bifurcation theory, Reviews of Modern Physics 63 (1991), no. 4, 991-1036.
[CWLJ94] P. Chang, B. Wang, T. Li, and L. Ji, *Interactions between the seasonal cycle and the Southern Oscillation-Frequency entrainment and chaos in a coupled ocean-atmosphere model*, Geophysical Research Letters **21** (1994), no. 25, 2817-2820.

[CZ95] R.F. Curtain and H. Zwart, *An Introduction to Infinite-Dimensional Linear Systems Theory*, vol. 21, Springer, 1995.

[CZ23] N. Chen and Y. Zhang, *Rigorous derivation of stochastic conceptual models for the El Niño-Southern Oscillation from a spatially-extended dynamical system*, Physica D: Nonlinear Phenomena **453** (2023), 133842.

[DCP$^{+}$21] B. Dieppois, A. Capotondi, B. Pohl, K.P. Chun, P.-A. Monerie, and J. Eden, *ENSO diversity shows robust decadal variations that must be captured for accurate future projections*, Communications Earth & Environment **2** (2021), no. 1, 212.

[DDMG16] M. Detrixhe, M. Doubeck, J. Moehlis, and F. Gibou, *A fast Eulerian approach for computation of global isochrons in high dimensions*, SIAM Journal on Applied Dynamical Systems **15** (2016), no. 3, 1501-1527.

[DGK$^{+}$08] A. Dhooge, W. Govaerts, Yu. A. Kuznetsov, H. G. E. Meijer, and B. Sautois, *New features of the software MatCont for bifurcation analysis of dynamical systems*, Mathematical and Computer Modelling of Dynamical Systems **14** (2008), 147-175.

[EL97] C. Eckert and M. Latif, *Predictability of a stochastically forced hybrid coupled model of El Niño*, Journal of Climate **10** (1997), no. 7, 1488-1504.

[ELR02] K. Engelborghs, T. Luzyanina, and D. Roosevelt, *Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL*, ACM Transactions on Mathematical Software **28** (2002), 1-21.

[EvP04] T. Eriola and J. von Pfaler, *Numerical Taylor expansions for invariant manifolds*, Numer. Math. **99** (2004), no. 1, 25-46.

[Fed02] A.V. Fedorov, *The response of the coupled tropical ocean-atmosphere to westerly wind bursts*, Quarterly Journal of the Royal Meteorological Society **128** (2002), no. 579, 1-23.

[FPS18] U. Feudel, A. N. Pisarchik, and K. Showalter, *Multistability and tipping: From mathematics and physics to climate and brain-Minireview and preface to the focus issue*, Chaos **28** (2018), 033501.

[FQS$^{+}$19] S. K. J. Falkena, C. Quinn, J. Sieber, J. Frank, and H. A. Dijkstra, *Derivation of delay equation climate models using the Mori-Zwanzig formalism*, Proceedings of the Royal Society A **475** (2019), no. 2227, 20190075.

[Gar04] C. W. Gardiner, *Handbook of Stochastic Methods*, 3rd ed., Springer-Verlag Berlin, 2004.

[GCS15] M. Ghil, M. D. Chekroun, and G. Stepan, *A collection on 'Climate Dynamics: Multiple Scales and Memory Effects', Editorial*, R. Soc. Proc. A **471** (2015), 20150097.

[GHJM98] L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, *Stochastic resonance*, Reviews of Modern Physics **70** (1998), 223-287.

[Gri08] A.S. Gritsun, *Unstable periodic trajectories of a barotropic model of the atmosphere*, Russ. J. Numer. Anal. Math. Modelling **43** (2008), 345-367.

[Gri13] A. Gritsun, *Statistical characteristics, circulation regimes and unstable periodic orbits of a barotropic atmospheric model*, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences **371** (2013), no. 1991, 20120336.

[GT00] E. Galanti and E. Tziperman, *ENSO's phase locking to the seasonal cycle in the fast-SST, fast-wave, and mixed-mode regimes*, J. Atmos. Sci. **57** (2000), 2936-2950.

[Guc75] J. Guckenheimer, *Isochrons and phaseless sets*, Journal of Mathematical Biology **1** (1975), no. 3, 259-273.

[GZT08] M. Ghil, I. Zaliapin, and S. Thompson, *A delay differential model of ENSO variability: parametric instability and the distribution of extremes*, Nonlinear Processes in Geophysics **15** (2008), no. 3, 417-433.

[Hal88] J. K. Hale, *Asymptotic Behavior of Dissipative Systems*, Mathematical Surveys and Monographs, vol. 25, American Mathematical Society, Providence, RI, 1988.

[HCF$^{+}$16] À. Haro, M. Canadell, J.-L. Figueras, A. Luque, and J.-M. Mondelo, *The Parameterization Method for Invariant Manifolds: From Rigorous Results to Effective Computations*, vol. 195, Springer-Verlag, Berlin, 2016.

[Hen81] D. Henry, *Geometric Theory of Semilinear Parabolic Equations*, Lecture Notes in Mathematics, vol. 840, Springer-Verlag, Berlin, 1981.

[HL84] W. Horsthemke and R. Lefever, *Noise-induced Transitions*, Springer Series in Synergetics, vol. 15, Springer-Verlag, Berlin, 1984, Theory and applications in physics, chemistry, and biology.

[HVL93] J. K. Hale and S. M. Verduyn Lunel, *Introduction to Functional-Differential Equations*, Applied Mathematical Sciences, vol. 99, Springer-Verlag, New York, 1993.

[IT86] K. Ito and R. Teglas, *Legendre-tau approximations for functional-differential equations*, SIAM J. Control Optim. **24** (1986), no. 4, 737-759.

[JN93] F.-F. Jin and J. D. Neelin, *Modes of interannual tropical ocean-atmosphere interaction-A unified view. Part I: Numerical results*, Journal of the atmospheric sciences **50** (1993), no. 21, 3477-3503.

[JNG94] F.-F. Jin, J. D. Neelin, and M. Ghil, *El Niño on the Devil's Staircase: Annual subharmonic steps to chaos*, Science **274** (1994), 70-72.
[JNG96] F.-F. Jin, J.D. Neelin, and M. Ghil, *El Niño/Southern Oscillation and the annual cycle: Subharmonic frequency-locking and aperiodicity*, Physica D: Nonlinear Phenomena **98** (1996), no. 2-4, 442-465.

[Kap86] F. Kappel, *Semigroups and delay equations*, Semigroups, Theory and Applications, Vol. II (Trieste, 1984), Pitman Res. Notes Math. Ser., vol. 152, Longman Sci. Tech., Harlow, 1986, pp. 136-176.

[KC18] D. Kondrashov and M. D. Chekroun, *Data-adaptive harmonic analysis and modeling of solar wind-magnetosphere coupling*, Journal of Atmospheric and Solar-Terrestrial Physics **177** (2018), 179-189.

[KCB18] D. Kondrashov, M. D. Chekroun, and P. Berloff, *Multiscale Stuart-Landau emulators: Application to wind-driven ocean gyres*, Fluids **3** (2018), no. 1, 21.

[KCG18] D. Kondrashov, M.D. Chekroun, and M. Ghil, *Data-adaptive harmonic decomposition and prediction of Arctic sea ice extent*, Dynamics and Statistics of the Climate System **3** (2018), no. 1, 1-23.

[KCYG18] D. Kondrashov, M. D. Chekroun, X. Yuan, and M. Ghil, *Data-adaptive harmonic decomposition and stochastic modeling of Arctic sea ice*, Advances in Nonlinear Geosciences (A. Tsonis, ed.), Springer, 2018, pp. 179-205.

[KF11] I. Koren and G. Feingold, *Aerosol-cloud-precipitation system as a predator-prey problem*, Proc. Natl. Acad. Sci. USA **108** (2011), no. 30, 12227-12232.

[KKD19] A. Keane, B. Krauskopf, and H. A. Dijkstra, *The effect of state dependence in a delay differential equation model for the El Niño Southern Oscillation*, Philosophical Transactions of the Royal Society A **377** (2019), no. 2153, 20180121.

[KKP15] A. Keane, B. Krauskopf, and C. Postlethwaite, *Delayed feedback versus seasonal forcing: Resonance phenomena in an El Niño Southern Oscillation model*, SIAM Journal on Applied Dynamical Systems **14** (2015), no. 3, 1229-1257.

[KKP16] A. Keane, B. Krauskopf, and C. M. Postlethwaite, *Investigating irregular behavior in a model for the El Niño Southern Oscillation with positive and negative delayed feedback*, SIAM J. Appl. Dyn. Syst. **15** (2016), 1656-1689.

[KKP17] __________, *Climate models with delay differential equations*, Chaos: An Interdisciplinary Journal of Nonlinear Science **27** (2017), no. 11, 114309.

[KKRG05] D. Kondrashov, S. Kravtsov, A. W. Robertson, and M. Ghil, *A hierarchy of data-based ENSO models*, Journal of Climate **18** (2005), 4425-4444.

[Knu] *KNUT*, http://rs1909.github.io/knut/, Accessed: 2023-09-09.

[Koo84] T. H. Koornwinder, *Orthogonal polynomials with weight function $(1-x)^\alpha(1+x)^\beta + M\delta(x+1) + N\delta(x-1)$*, Canad. Math. Bull. **27** (1984), 205-214.

[KS78] F. Kappel and W. Schappacher, *Autonomous nonlinear functional differential equations and averaging approximations*, Nonlinear Analysis: Theory, Methods & Applications **2** (1978), no. 4, 391-422.

[KS14] B. Krauskopf and J. Sieber, *Bifurcation analysis of delay-induced resonances of the El-Niño Southern Oscillation*, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences **470** (2014), no. 2169, 20140348.

[KTF17] I. Koren, E. Tziperman, and G. Feingold, *Exploring the nonlinear cloud and rain equation*, Chaos: An Interdisciplinary Journal of Nonlinear Science **27** (2017), no. 1, 013107.

[Kue11] Ch. Kuehn, *A mathematical framework for critical transitions: Bifurcations, fast-slow systems and stochastic dynamics*, Physica D **240** (2011), no. 12, 1020-1035.

[LG20] V. Lucarini and A. Gritsun, *A new mathematical framework for atmospheric blocking events*, Climate Dynamics **54** (2020), no. 1-2, 575-598.

[LTW$^{+}$17] M. L. L'Heureux, K. Takahashi, A. B. Watkins, A. G. Barnston, E. J. Becker, T. E. Di Liberto, F. Gamble, J. Gottschalck, M. S. Halpert, B. Huang, et al., *Observing and predicting the 2015/16 El Niño*, Bulletin of the American Meteorological Society **98** (2017), no. 7, 1363-1382.

[Mat] *MatCont*, https://sourceforge.net/projects/matcont/, Accessed: 2023-09-09.

[MRMM14] A. Mauroy, B. Rhoads, J. Moehlis, and I. Mezic, *Global isochrons and phase sensitivity of bursting neurons*, SIAM Journal on Applied Dynamical Systems **13** (2014), no. 1, 306-338.

[MW05] T. Ma and S. Wang, *Bifurcation Theory and Applications*, World Scientific Series on Nonlinear Science. Series A: Monographs and Treatises, vol. 53, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2005.

[MW14] __________, *Phase Transition Dynamics*, Springer, 2014.

[NBH$^{+}$98] J. D. Neelin, D. S. Battisti, A. C. Hirst, F.-F. Jin, Y. Wakata, T. Yamagata, and S. E. Zebiak, *ENSO theory*, J. Geophys. Res. Oceans **103** (1998), 14261-14290.

[NJ93] J. D. Neelin and F.-F. Jin, *Modes of interannual tropical ocean-atmosphere interaction-a unified view. Part II: Analytical results in the weak-coupling limit*, Journal of the atmospheric sciences **50** (1993), no. 21, 3504-3522.

[PS95] C. Penland and P. D. Sardeshmukh, *The optimal growth of tropical sea surface temperature anomalies*, Journal of climate **8** (1995), no. 8, 1999-2024.

[RCC$^{+}$14] L. Roques, M. D. Chekroun, M. Cristofol, S. Soubeyrand, and M. Ghil, *Parameter estimation for energy balance models with memory*, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences **470** (2014), no. 2169, 20140349.
[RN00] M.S. Roulston and J.D. Neelin, *The response of an ENSO model to climate noise, weather noise and intraseasonal forcing*, Geophysical Research Letters **27** (2000), no. 22, 3723-3726.

[SBB<sup>+</sup>09] M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, S. R. Carpenter, V. Dakos, H. Held, E. H. Van Nes, M. Rietkerk, and G. Sugihara, *Early-warning signals for critical transitions*, Nature **461** (2009), 53-59.

[Sel85] G. R. Sell, *Smooth linearization near a fixed point*, American Journal of Mathematics (1985), 1035-1091.

[SEL<sup>+</sup>14] J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, and D. Roosevelt, *DDE-BIFTOOL v. 3.1.1 Manual—Bifurcation analysis of delay differential equations*, arXiv preprint arXiv:1406.7144 (2014).

[SS88] M. J. Suarez and P. S. Schopf, *A delayed action oscillator for ENSO*, J. Atmos. Sci. **45** (1988), 3283-3287.

[SY02] G. Sell and Y. You, *Dynamics of Evolutionary Equations*, Applied Mathematical Sciences, vol. 143, Springer-Verlag, New York, 2002.

[TCDN20] A. Tantet, M.D. Chekroun, H. Dijkstra, and J. D. Neelin, *Ruelle-Pollicott Resonances of Stochastic Systems in Reduced State Space. Part III: Application to the Cane-Zebiak Model of the El Nino-Southern Oscillation*, J. Stat. Phys. **179** (2020), 1449-1474, doi: 10.1007/s10955-019-02444-8.

[TMCS16] S. Thual, A. J. Majda, N. Chen, and S. N. Stechmann, *Simple stochastic model for El Niño with westerly wind bursts*, Proc. Natl. Acad. Sci. USA **113** (2016), 10245-10250.

[TSCJ94] E. Tziperman, L. Stone, M. A. Cane, and H. Jarosh, *El Niño chaos: Overlapping of resonances between the seasonal cycle and the Pacific ocean-atmosphere oscillator*, Science **264** (1994), 72-74.

[Vya12] C. P. Vyasarayani, *Galerkin approximations for higher order delay differential equations*, J. Comput. Nonlinear Dynamics **7** (2012), 031004.

[WC05] P. Wahi and A. Chatterjee, *Galerkin projections for delay differential equations*, ASME. J. Dyn. Syst., Meas., Control **127** (2005), 80-87.

[ZGMPK03] J. Zavala-Garay, A.M. Moore, C.L. Perez, and R. Kleeman, *The response of a coupled model of ENSO to observed estimates of stochastic forcing*, Journal of Climate **16** (2003), no. 17, 2827-2842.

(MC) DEPARTMENT OF EARTH AND PLANETARY SCIENCES, WEIZMANN INSTITUTE, REHOVOT 76100, ISRAEL, DEPARTMENT OF ATMOSPHERIC AND OCEANIC SCIENCES, UNIVERSITY OF CALIFORNIA, LOS ANGELES, CA 90095-1565, USA

Email address: mchekroun@atmos.ucla.edu

(HL) DEPARTMENT OF MATHEMATICS, VIRGINIA TECH, BLACKSBURG, VA 24061, USA

Email address: hhliu@vt.edu