# Effective Reduced Models from Delay Differential Equations:
Bifurcations, Tipping Solution Paths, and ENSO variability

Mickaël D. Chekroun  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>and  Honghu Liu
Department of Mathematics, Virginia Tech, Blacksburg, VA 24061, USA
<EMAIL_ADDRESS>
(Date: January 2, 2024)

###### 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

  1. 1 Introduction
  2. 2 Galerkin-Koornwinder (GK) approximations of DDEs
  3. 3 Center-Unstable Manifold Approximations from GK systems
  4. 4 Bifurcations Analysis of a Delay El Niño-Southern Oscillation (ENSO) Model
  5. 5 Tipping Solution Paths and ENSO variability
  6. 6 Concluding Remarks
  7. A Determination of $\mathcal{P}_{N}$ and $\mathcal{Q}_{N}$ in Eq. (4.10)

##  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+98],
delay models have attracted a growing attention in climate modeling; see
[GT00, GZT08, KF11, RCC+14, GCS15, KTF17, KKP17, CGN18, BCL+17, FQS+19].
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+19, 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+16] and have
shown relevance in computing characteristic roots [BMV05], Lyapunov exponents
[BVV14], and bifurcation analysis [ABL+22] using ODE bifurcation software
packages such as AUTO [Aut] and MatCont [DGK+08, Mat]. These methods, although
showing great promises, do not have yet developed into toolboxes such as DDE-
BIFTOOL [ELR02, SEL+14] 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+16].

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)
$\frac{\,\mathrm{d}x(t)}{\,\mathrm{d}t}=ax(t)+bx(t-\tau)+c\int_{t-\tau}^{t}x(s)\,\mathrm{d}s+F\Big{(}x(t),x(t-\tau),\int_{t-\tau}^{t}x(s)\,\mathrm{d}s\Big{)},$

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)
$K_{n}(s)=-(1+s)\frac{\,\mathrm{d}}{\,\mathrm{d}s}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(\,\mathrm{d}x)=\frac{1}{2}\,\mathrm{d}x+\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) $\displaystyle\int_{-1}^{1}K_{n}(s)K_{p}(s)\,\mathrm{d}\mu(s)$
$\displaystyle=\frac{1}{2}\int_{-1}^{1}K_{n}(s)K_{p}(s)\,\mathrm{d}s+K_{n}(1)K_{p}(1)$
$\displaystyle=0,\,\mbox{ if }p\neq n.$

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

(2.4) $\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)
$\Big{\langle}(f,p),(g,q)\Big{\rangle}_{\mathcal{E}}=\frac{1}{2}\int_{-1}^{1}f(s)g(s)\,\mathrm{d}s+pq,\qquad\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)
$\|\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 $(-\tau,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) $\displaystyle K^{\tau}_{n}\colon[-\tau,0]$
$\displaystyle\rightarrow\mathbb{R},$ $\displaystyle\theta$
$\displaystyle\mapsto K_{n}\Bigl{(}1+\frac{2\theta}{\tau}\Bigr{)}.$

They form orthogonal polynomials on the interval $[-\tau,0]$ for the
$L^{2}$-inner product on $(-\tau,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) $\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 $(-\tau,0)$ and weighted by $1/\tau$. Note that
since $K_{n}(1)=1$ [CGLW16, Prop. 3.1], we have also

(2.9) $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 $[-\tau,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) $\partial_{t}u=\partial_{\theta}u,\quad-\tau\leq\theta<0,$

subject to the nonlocal and nonlinear boundary condition

(2.11) $\displaystyle\partial_{\theta}u|_{\theta=0}$
$\displaystyle=au(t,0)+bu(t,-\tau)+c\int_{t-\tau}^{t}u(s,0)\,\mathrm{d}s$
$\displaystyle\qquad\;\;+F\Big{(}u(t,0),u(t,-\tau),\int_{t-\tau}^{t}u(s,0)\,\mathrm{d}s\Big{)},\qquad\mbox{for
$t\geq 0$.}$

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

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$.

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) $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) $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\rightarrow\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) $\displaystyle\frac{\,\mathrm{d}y_{j}}{\,\mathrm{d}t}$
$\displaystyle=\frac{1}{\|\mathcal{K}_{j}\|_{\mathcal{E}}^{2}}\sum_{n=0}^{N-1}\bigg{(}a+bK_{n}(-1)+c\tau(2\delta_{n,0}-1)$
$\displaystyle\hskip
80.00012pt+\frac{2}{\tau}\sum_{k=0}^{n-1}a_{n,k}\left(\delta_{j,k}\|\mathcal{K}_{j}\|^{2}_{\mathcal{E}}-1\right)\bigg{)}y_{n}(t)$
$\displaystyle\hskip
10.00002pt+\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),$ $\displaystyle\;\;0\leq j\leq
N-1.$

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)
$\boxed{\frac{\,\mathrm{d}\bm{y}}{\,\mathrm{d}t}=A_{N}(\tau)\bm{y}+G_{2}(\bm{y},\tau),}$

where $\bm{y}=(y_{0},\cdots,y_{N-1})^{T}$, and in which the terms
$A_{N}(\tau)\bm{y}$ and $G_{2}(\bm{y},\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) $\boxed{\begin{aligned}
\Big{(}A_{N}(\tau)\Big{)}_{i,j}&=\frac{1}{\|\mathcal{K}_{i}\|_{\mathcal{E}}^{2}}\Big{(}a+bK_{j}(-1)+c\tau(2\delta_{j,0}-1)\\\
&\hskip
80.00012pt+\frac{2}{\tau}\sum_{k=0}^{j-1}a_{j,k}\left(\delta_{i,k}\|\mathcal{K}_{i}\|^{2}_{\mathcal{E}}-1\right)\Big{)},\;\;i,j=0,\cdots,N-1.\end{aligned}}$

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\mathcal{P}_{N}+\mathcal{Q}_{N}$, one
observes that only $\mathcal{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)
$\boxed{F_{N}^{j}(\bm{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\leq j\leq 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) $\frac{\,\mathrm{d}\bm{y}}{\,\mathrm{d}t}=A(\tau)\bm{y}+G(\bm{y}),$

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

(3.2) $G(0)=0,\qquad\text{ and }\qquad 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(\bm{y})$ admits the
following Taylor expansion for $\bm{y}$ near the origin:

(3.3) $G(\bm{y})=G_{k}(\underbrace{\bm{y},\cdots,\bm{y}}_{k\text{
times}})+O(\|\bm{y}\|^{k+1}),$

where

(3.4) $\displaystyle
G_{k}\colon\underbrace{\mathbb{R}^{N}\times\cdots\times\mathbb{R}^{N}}_{k\text{
times}}\rightarrow\mathbb{R}^{N}$

denotes a homogenous polynomial of order $k\geq 2$. Often, $G_{k}(\bm{y})$ is
used hereafter as a compact notation for $G_{k}(\bm{y},\,\cdots\,,\bm{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) $\displaystyle\sigma(A(\tau))=\\{\lambda_{n}(\tau),\;1\leq n\leq N\\},$

with, for any $1\leq n<n^{\prime}$, either

(3.6)
$\displaystyle\,\mathrm{Re}\,\lambda_{n}(\tau)>\,\mathrm{Re}\,\lambda_{n^{\prime}}(\tau),$

or

(3.7)
$\displaystyle\,\mathrm{Re}\,\lambda_{n}(\tau)=\,\mathrm{Re}\,\lambda_{n^{\prime}}(\tau),\;\text{
and
}\;\,\mathrm{Im}\,\lambda_{n}(\tau)\geq\,\mathrm{Im}\,\lambda_{n^{\prime}}(\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) $\displaystyle\,\mathrm{Re}\,\lambda_{j}(\tau)\begin{cases}<0&\mbox{if
}\tau<\tau_{c},\\\ =0&\mbox{if }\tau=\tau_{c},\\\ >0&\mbox{if
}\tau>\tau_{c},\end{cases}$ $\displaystyle 1\leq j\leq m_{c},$
$\displaystyle\,\mathrm{Re}\,\lambda_{j}(\tau_{c})<0,$ $\displaystyle j\geq
m_{c}+1,$

for some $m_{c}>0$, and for $\tau$ in some neighborhood $\mathcal{U}$ of
$\tau_{c}$. Of course, the PES holds also when one crosses $\tau_{c}$ from
above with $\,\mathrm{Re}\,\lambda_{j}(\tau)<0$ for $\tau>\tau_{c}$, while
$\,\mathrm{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)
$\displaystyle\sigma_{\mathfrak{c}}(A(\tau))=\\{\lambda_{j}(\tau)\>|\>j=1,\cdots,m_{c}\\},$
$\displaystyle\sigma_{\mathfrak{s}}(A(\tau))=\\{\lambda_{j}(\tau)\>|\>j=m_{c}+1,\cdots,N\\}.$

The PES condition prevents other eigenvalues in
$\sigma_{\mathfrak{s}}(A(\tau))$ from crossing the imaginary axis as $\tau$
varies in $\mathcal{U}$. Hence, no eigenvalues other than those of
$\sigma_{\mathfrak{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) $0>2k\eta_{\mathfrak{c}}>\eta_{\mathfrak{s}},$

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

$\displaystyle\eta_{\mathfrak{c}}=\inf_{\tau\in\mathcal{U}}\inf_{j=1,\cdots,m_{c}}\mathrm{Re}(\lambda_{j}(\tau)),$
$\displaystyle\eta_{\mathfrak{s}}=\sup_{\tau\in\mathcal{U}}\sup_{j\geq
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_{\mathfrak{c}}$, given by

(3.11) $H_{\mathfrak{c}}=\mathrm{span}\\{\bm{e}_{1},\cdots,\bm{e}_{m_{c}}\\},$

while a mode $\bm{e}_{n}$ with $n\geq 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_{\mathfrak{s}}$ the subspace spanned by these stable modes. The projector
onto the subspace $H_{\mathfrak{s}}$ (resp. $H_{\mathfrak{c}}$) is denoted by
$\Pi_{\mathfrak{s}}$ (resp. $\Pi_{\mathfrak{c}}$). Similarly,
$A_{\mathfrak{c}}(\tau)$ (resp. $A_{\mathfrak{s}}(\tau)$) denotes
$\Pi_{\mathfrak{c}}A(\tau)$ (resp. $\Pi_{\mathfrak{s}}A(\tau)$), and
$\bm{y}_{c}$ (resp. $\bm{y}_{\mathfrak{s}}$) denotes the vector in
$H_{\mathfrak{c}}$ (resp. $H_{\mathfrak{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)
$\langle\bm{a},\bm{b}\rangle=\sum_{i=1}^{N}a_{i}b_{i}^{\ast},\;\bm{a},\bm{b}\in\mathbb{C}^{N}.$

In what follows we also denote by $\bm{e}_{j}^{\ast}$ the eigenmodes of the
conjugate transpose $A(\tau)^{\ast}$ 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) $\gamma_{m_{c}}(\tau)\geq C\mbox{Lip}(G|_{\mathcal{V}}),$

where $\mbox{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 parameterization111As
provided for instance by a center manifold or the unstable manifold of the
origin. 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.

######  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}=\mbox{graph}(h_{\tau})$, with $h_{\tau}$ that maps
$H_{\mathfrak{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) $\displaystyle\text{
}\forall\>\,(j_{1},\cdots,j_{k})\in(1,\cdots,m_{c})^{k},\;n\geq m_{c}+1,\text{
it holds that}$ $\displaystyle\Bigl{(}\langle
G_{k}(\bm{e}_{j_{1}},\cdots,\bm{e}_{j_{k}}),\bm{e}_{n}^{\ast}\rangle\neq
0\Bigr{)}\Longrightarrow\biggl{(}\,\mathrm{Re}\,(\sum_{\ell=1}^{k}\lambda_{j_{\ell}}-\lambda_{n})\neq
0\biggr{)},$

where $G_{k}$ denotes the leading-order term in the Taylor expansion of $G$.

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

(3.14) $\mathcal{L}_{A(\tau)}[\psi](\bm{X})=\Pi_{\mathfrak{s}}G_{k}(\bm{X}),$

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

(3.15)
$\mathcal{L}_{A(\tau)}[\psi](\bm{X})=D\psi(\bm{X})A_{\mathfrak{c}}(\tau)\bm{X}-A_{\mathfrak{s}}(\tau)\psi(\bm{X}),\;\bm{X}\in
H_{\mathfrak{c}}.$

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

(3.16) $\left\lVert
h_{\tau}(\bm{X})-\Phi_{\tau}(\bm{X})\right\rVert=o(\left\lVert\bm{X}\right\rVert^{k}),$

as long as $\bm{X}$ lies in a neighborhood $\mathcal{N}(\tau)$ of the origin
in $H_{\mathfrak{c}}$ spanned by the $m_{c}$ eigenmodes $\bm{e}_{j}$ losing
their stability (see (3.11)), as $\tau$ crosses $\tau_{c}$.

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

(3.17) $\Phi_{\tau}(X)=\sum_{n\geq
m_{c}+1}\Phi_{n,\tau}(\bm{X})\bm{e}_{n},\;\qquad\bm{X}\in H_{\mathfrak{c}},$

where

(3.18) $\Phi_{n,\tau}(\bm{X})=\sum_{1\leq j_{1},\cdots,j_{k}\leq
m_{c}}G_{j_{1}\cdots j_{k}}^{n}M_{j_{1}\cdots j_{k}}^{n,\tau}X_{j_{1}}\cdots
X_{j_{k}},\;\qquad X_{j}=\langle\bm{X},\bm{e}_{j}\rangle,$

with

(3.19) $M_{j_{1}\cdots
j_{k}}^{n,\tau}=\Big{(}\sum_{\ell=1}^{k}\lambda_{j_{\ell}}(\tau)-\lambda_{n}(\tau)\Big{)}^{-1},$

and the $G_{j_{1}\cdots j_{k}}^{n}$ denoting the coefficients accounting for
the magnitude carried out by $\bm{e}_{n}^{\ast}$, of the nonlinear
interactions through the leading-order term $G_{k}$, between the low modes
$\bm{e}_{j_{1}}$, $\cdots$, $\bm{e}_{j_{k}}$ (in $H_{c}$), namely:

(3.20) $G_{j_{1}\cdots j_{k}}^{n}=\langle
G_{k}(\bm{e}_{j_{1}},\cdots,\bm{e}_{j_{k}}),\bm{e}_{n}^{\ast}\rangle,\;\;1\leq
j_{1},\cdots,j_{k}\leq 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_{\mathfrak{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_{\mathfrak{c}}(\tau)}\bm{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 $\bm{y}(t)$ to Eq. (3.1) such that $\bm{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)
$\bm{y}(t)=\bm{y}_{\mathfrak{c}}(t)+h_{\tau}(\bm{y}_{\mathfrak{c}}(t)),\;t\in[t_{0},t_{0}+\sigma),$

where $\bm{y}_{\mathfrak{c}}(t)$ denotes the projection of $\bm{y}(t)$ onto
the subspace $H_{\mathfrak{c}}$.

  * (ii)

If there exists a trajectory $t\mapsto\bm{y}(t)$ such that
$\bm{y}_{\mathfrak{c}}(t)$ belongs, for all $-\infty<t<\infty$, to the
neighborhood $\mathfrak{B}(\tau)$ (in $H_{\mathfrak{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_{\mathfrak{c}}$ is contained in $\mathfrak{B}(\tau)$, i.e. if
$\Pi_{\mathfrak{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,
$\bm{y}_{\mathfrak{c}}(t)$, is sufficient to entirely determine any solution
$\bm{y}(t)$ that belongs to such an invariant set; see [SY02, Theorem 71.4].
Furthermore, $\bm{y}_{\mathfrak{c}}(t)$ is obtained as the solution of the
following reduced $m_{c}$-dimensional problem

(3.22)
$\dot{\widehat{\bm{y}}}=A_{\mathfrak{c}}(\tau)\widehat{\bm{y}}+\Pi_{\mathfrak{c}}G(\widehat{\bm{y}}+h_{\tau}(\widehat{\bm{y}})),\qquad\widehat{\bm{y}}(0)=\bm{y}_{\mathfrak{c}}(0)\in\mathfrak{B}(\tau),$

which in turn characterizes the solution $\bm{y}(t)$ in $\Sigma$, since the
slaving relationship
$\bm{y}_{\mathfrak{s}}(t)=h_{\tau}(\bm{y}_{\mathfrak{c}}(t))$ holds for any
solution $\bm{y}(t)$ that belongs to an invariant set $\Sigma$ for which
$\Pi_{\mathfrak{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)
$\dot{\bm{x}}=A_{\mathfrak{c}}(\tau)\bm{x}+\Pi_{\mathfrak{c}}G_{k}(\bm{x}+\Phi_{\tau}(\bm{x})),\qquad\bm{x}(0)=\bm{y}_{\mathfrak{c}}(0)\in\mathfrak{B}(\tau),$

provides an approximation of the invariant set $\Sigma$ for which
$\Pi_{\mathfrak{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=\mbox{dim}(H_{\mathfrak{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\rightarrow\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 $\widehat{\bm{y}}(t)$ (and thus of $\bm{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)
$Dh(\bm{X})[A_{\mathfrak{c}}\bm{X}+\Pi_{\mathfrak{c}}G(\bm{X}+h(\bm{X}))]-A_{\mathfrak{s}}h(\bm{X})=\Pi_{\mathfrak{s}}G(\bm{X}+h(\bm{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. [BK98, 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_{\mathfrak{c}}$ contains typically a large fraction of
the solution’s total energy. We illustrate this idea to the case where
$G(\bm{y})=G_{2}(\bm{y},\bm{y})+G_{3}(\bm{y},\bm{y},\bm{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)
$\mathcal{L}_{A}[\psi](\bm{X})=-D\Phi_{\tau}(\bm{X})\Pi_{\mathfrak{c}}G_{2}(\bm{X})+\Pi_{\mathfrak{s}}G_{2}(\bm{X},\Phi_{\tau}(\bm{X}))+\Pi_{\mathfrak{s}}G_{2}(\Phi_{\tau}(\bm{X}),\bm{X})+\Pi_{\mathfrak{s}}G_{3}(\bm{X}),$

where
$\mathcal{L}_{A}[\psi](\bm{X})=D\psi(\bm{X})A_{\mathfrak{c}}\bm{X}-A_{\mathfrak{s}}\psi(\bm{X})$;
see (3.15).

Note that the RHS is a homogeneous cubic polynomial in the $\bm{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
$\bm{y}_{\mathfrak{s}}$ is typically much smaller than
$\|\bm{y}_{\mathfrak{c}}\|^{2}$. It is then reasonable to expect that the
energy carried by $\Phi_{\tau}(\bm{X})$ is much smaller than $\|\bm{X}\|^{2}$
for $\bm{X}=\bm{y}_{\mathfrak{c}}$. Thus, one may assume that in the RHS of
(3.25), the term $\Pi_{\mathfrak{s}}G_{3}(\bm{X})$ dominates the other three
terms provided that $\|G_{2}(\bm{y})\|/\|\bm{y}\|^{2}$ is on the same order of
magnitude as $\|G_{3}(\bm{y})\|/\|\bm{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)
$Dh_{3}(\bm{X})A_{\mathfrak{c}}\bm{X}-A_{\mathfrak{s}}h_{3}(\bm{X})=\Pi_{\mathfrak{s}}G_{3}(\bm{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) $\displaystyle\Psi_{\tau}(\bm{X})$
$\displaystyle=\Phi_{\tau}(\bm{X})+h_{3}(\bm{X})=\sum_{n=m_{c}+1}^{N}\Psi_{n,\tau}(X)\bm{e}_{n}$
$\displaystyle=\sum_{n=m_{c}+1}^{N}\bigg{(}\sum_{1\leq j_{1},j_{2}\leq
m_{c}}\frac{\langle
G_{2}(\bm{e}_{j_{1}},\bm{e}_{j_{2}}),\bm{e}_{n}^{\ast}\rangle}{\lambda_{j_{1}}(\tau)+\lambda_{j_{2}}(\tau)-\lambda_{n}(\tau)}X_{j_{1}}X_{j_{2}}$
$\displaystyle\hskip 56.9055pt+\sum_{1\leq j_{1},j_{2},j_{3}\leq
m_{c}}\frac{\langle
G_{3}(\bm{e}_{j_{1}},\bm{e}_{j_{2}},\bm{e}_{j_{3}}),\bm{e}_{n}^{\ast}\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}}\bigg{)}\bm{e}_{n}.$

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) $\frac{\,\mathrm{d}T}{\,\mathrm{d}t}=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) $T_{0}=0,\qquad T_{+}=\sqrt{1-\alpha},\qquad 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) $\,\widehat{T}=T-T_{+},$

namely,

(4.4)
$\frac{\,\mathrm{d}\,\widehat{T}}{\,\mathrm{d}t}=(1-3T_{+}^{2})\,\widehat{T}(t)-\alpha\,\widehat{T}(t-\tau)-3T_{+}\,\widehat{T}^{2}(t)-\,\widehat{T}^{3}(t).$

The above equation fits into Eq. (2.1) with

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

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

(4.6)
$\frac{\,\mathrm{d}\bm{y}}{\,\mathrm{d}t}=A(\tau)\bm{y}+G(\bm{y}),\qquad\bm{y}=(y_{1},\cdots,y_{N})^{T},$

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

(4.7) $\displaystyle(A(\tau))_{j,n}$
$\displaystyle=\frac{1}{\|\mathcal{K}_{j}\|_{\mathcal{E}}^{2}}\sum_{n=0}^{N-1}\Bigl{(}1-3T_{+}^{2}-\alpha
K_{n}(-1)$
$\displaystyle\qquad+\frac{2}{\tau}\sum_{k=0}^{n-1}a_{n,k}\left(\delta_{j,k}\|\mathcal{K}_{j}\|^{2}_{\mathcal{E}}-1\right)\Bigr{)},\qquad\mbox{for
$j,n=0,\cdots,N-1$.}$

As in Sec. 2.3, the $a_{n,k}$ are determined thanks to Proposition 1 recalled
in Appendix A, and the $\|\mathcal{K}_{j}\|_{\mathcal{E}}$ 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(\bm{y})=G_{2}(\bm{y})+G_{3}(\bm{y}),$

with

(4.8)
$G_{2}(\bm{y})=-3T_{+}\left(\sum_{n=0}^{N-1}y_{n}\right)^{2}\bm{\nu}_{N},\mbox{
and }\;G_{3}(\bm{y})=-\left(\sum_{n=0}^{N-1}y_{n}\right)^{3}\bm{\nu}_{N},$

where $\bm{\nu}_{N}$ denotes the $N$-dimensional column vector

(4.9)
$\bm{\nu}_{N}=\Big{(}\frac{1}{\|\mathcal{K}_{0}\|_{\mathcal{E}}^{2}},\cdots,\frac{1}{\|\mathcal{K}_{N-1}\|_{\mathcal{E}}^{2}}\Big{)}^{T}$

We arrive finally at the following $N$-dimensional GK approximation for the
DDE (4.4)

(4.10)
$\displaystyle\frac{\,\mathrm{d}\bm{y}}{\,\mathrm{d}t}=\left(\frac{2}{\tau}\mathcal{P}_{N}+\mathcal{Q}_{N}\right)\bm{y}-\left(\sum_{n=1}^{N}y_{n}\right)^{2}\left(3T_{+}+\sum_{n=1}^{N}y_{n}\right)\bm{\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 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}$. GK dimension | $\tau=1.562$ | $\tau=1.6$ | $\tau=1.7$ | $\tau=1.9$
---|---|---|---|---
$N=4$ | $5.66\times 10^{-2}$ | $3.87\times 10^{-2}$ | $4.29\times 10^{-2}$ | $5.71\times 10^{-2}$
$N=6$ | $1.44\times 10^{-4}$ | $1.80\times 10^{-4}$ | $7.75\times 10^{-4}$ | $3.60\times 10^{-3}$
$N=8$ | $1.54\times 10^{-4}$ | $6.95\times 10^{-5}$ | $2.13\times 10^{-5}$ | $4.37\times 10^{-4}$
$N=10$ | $1.31\times 10^{-4}$ | $5.60\times 10^{-5}$ | $3.65\times 10^{-5}$ | $8.13\times 10^{-5}$

###  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)
$\dot{\bm{x}}=A_{\mathfrak{c}}(\tau)\bm{x}+\Pi_{\mathfrak{c}}G\bigl{(}\bm{x}+\Psi_{\tau}(\bm{x})\bigr{)},$

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)
$\dot{y}_{j}=\lambda_{j}(\tau)y_{j}+\Bigl{\langle}G\Bigl{(}\sum_{\ell=1}^{N}y_{\ell}\bm{e}_{\ell}\Bigr{)},\bm{e}_{j}^{\ast}\Bigr{\rangle},\;\;j=1,\cdots,N,$

where $y_{j}=\langle\bm{y},\bm{e}_{j}^{*}\rangle$, with inner product defined
in (3.12).

Figure 2. Eigenvalues dependence as $\tau$ crosses its critical value
$\tau_{c}$. Are shown here, the first 10 pairs
$(\lambda_{j}(\tau),\overline{\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$.

We observe that a dominant pair of complex eigenvalues crosses the imaginary
axis as $\tau$ crosses from below the critical value
$\tau_{\mathfrak{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_{\mathfrak{c}}$ and $H_{\mathfrak{s}}$ as
follows:

(4.13) $H_{\mathfrak{c}}=\mathrm{span}\\{\bm{e}_{1},\bm{e}_{2}\\},\qquad
H_{\mathfrak{s}}=\mathrm{span}\\{\bm{e}_{3},\cdots,\bm{e}_{N}\\},$

in which $(\bm{e}_{1},\bm{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) $\displaystyle\dot{x}_{j}=\lambda_{j}(\tau)x_{j}+\langle
G(\bm{x}+\Psi_{\tau}(\bm{x})),\bm{e}_{j}^{\ast}\rangle,\qquad j=1,2,$

with $\bm{x}=x_{1}\bm{e}_{1}+x_{2}\bm{e}_{2}$, and where the parameterization
$\Psi_{\tau}$ of the neglected modes in $H_{\mathfrak{s}}$ is given by (3.27).

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) $\displaystyle G_{N}^{(2)}(\bm{p},\bm{q})$
$\displaystyle=-3T_{+}\left(\sum_{n=0}^{N-1}p_{n}\right)\left(\sum_{n=0}^{N-1}q_{n}\right)\bm{\nu}_{N}$
$\displaystyle G_{N}^{(3)}(\bm{p},\bm{q},\bm{r})$
$\displaystyle=-\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)\bm{\nu}_{N}.$

The Lyapunov coefficient $\ell^{N}_{1}(\tau_{c})$ is then given by

(4.16)
$\ell_{1}^{N}(\tau_{c})=\mathrm{Re}\Bigl{(}\frac{a^{N}_{20}a^{N}_{11}\sqrt{-1}}{\mathrm{Im}(\lambda_{1}(\tau_{c}))}+a^{N}_{21}\Bigr{)},$

where

(4.17) $\displaystyle a_{20}^{N}$ $\displaystyle=\langle
G_{N}^{(2)}(\bm{e}_{1},\bm{e}_{1}),\bm{e}_{1}^{\ast}\rangle,$ $\displaystyle
a^{N}_{11}$ $\displaystyle=\langle
G_{N}^{(2)}(\bm{e}_{1},\bm{e}_{2}),\bm{e}_{1}^{\ast}\rangle+\langle
G_{N}^{(2)}(\bm{e}_{2},\bm{e}_{1}),\bm{e}_{1}^{\ast}\rangle,$

and

(4.18) $\displaystyle a_{21}^{N}$ $\displaystyle=\langle
G_{N}^{(3)}(\bm{e}_{1},\bm{e}_{1},\bm{e}_{2}),\bm{e}_{1}^{\ast}\rangle+\langle
G_{N}^{(3)}(\bm{e}_{1},\bm{e}_{2},\bm{e}_{1}),\bm{e}_{1}^{\ast}\rangle+\langle
G_{N}^{(3)}(\bm{e}_{2},\bm{e}_{1},\bm{e}_{1}),\bm{e}_{1}^{\ast}\rangle$
$\displaystyle\qquad+\sum_{n=3}^{Nd}\frac{\langle
G_{N}^{(2)}(\bm{e}_{1},\bm{e}_{2}),\bm{e}_{n}^{\ast}\rangle+\langle
G_{N}^{(2)}(\bm{e}_{2},\bm{e}_{1}),\bm{e}_{n}^{\ast}\rangle}{2\mathrm{Re}(\lambda_{1}(\tau_{c}))-\lambda_{n}(\tau_{c})}\Bigl{[}\langle
G_{N}^{(2)}(\bm{e}_{1},\bm{e}_{n}),\bm{e}_{1}^{\ast}\rangle+\langle
G_{N}^{(2)}(\bm{e}_{n},\bm{e}_{1}),\bm{e}_{1}^{\ast}\rangle\Bigr{]}$
$\displaystyle\qquad+\sum_{n=3}^{Nd}\frac{\langle
G_{N}^{(2)}(\bm{e}_{1},\bm{e}_{1}),\bm{e}_{n}^{\ast}\rangle}{2\lambda_{1}(\tau_{c})-\lambda_{n}(\tau_{c})}\Bigl{[}\langle
G_{N}^{(2)}(\bm{e}_{2},\bm{e}_{n}),\bm{e}_{1}^{\ast}\rangle+\langle
G_{N}^{(2)}(\bm{e}_{n},\bm{e}_{2}),\bm{e}_{1}^{\ast}\rangle\Bigr{]}.$

Theorem III.1 of [CKL20] ensures then that subcritical Hopf bifurcation occurs
for the $N$-dimensional GK system (4.12) if $\ell^{N}_{1}(\tau_{c})>0$ and a
supercritical Hopf bifurcation occurs if $\ell^{N}_{1}(\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^{N}_{1}(\tau_{c})$, in
terms of the GK dimension $N$. As can be observed, both $\tau_{c}$ and
$\ell^{N}_{1}(\tau_{c})$ converges quickly as the GK dimension $N$ increases.
Since $\ell^{N}_{1}(\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 Laypunov coefficient $\ell^{N}_{1}(\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^{N}_{1}(\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^{\text{true}}_{c}$ given by [BTR07, Eq. (7)] from
the characteristic equation,

(4.19)
$\tau^{\text{true}}_{c}=\frac{\arccos(\frac{3\alpha-2}{\alpha})}{\sqrt{\alpha^{2}-(3\alpha-2)^{2}}}.$

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

Table 2. The estimated $\tau_{c}$ and Lyapunov coefficient $\ell^{N}_{1}(\tau_{c})$ as GK dimension $N$ increases. | $N=4$ | $N=6$ | $N=8$ | $N=10$ | $N=12$
---|---|---|---|---|---
$\tau_{c}$ | 1.7343471 | 1.7408640 | 1.7408394 | 1.7408395 | 1.7408395
$\ell^{N}_{1}(\tau_{c})$ | 2.1884711 | 2.2246906 | 2.2247594 | 2.2247568 | 2.2247568
Figure 3. Lyapunov coefficient $\ell^{N}_{1}(\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).

###  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^{\sharp}$ 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^{\sharp}$; see
cyan curve in Fig. 4A. The presence of this homoclinic orbit is manifested by
the jump displayed by the curve of UPO’s amplitudes222Obtained by
approximating $\,\widehat{T}(t)$ solving Eq. (4.4) by the UPOs from the
reduced equation according to (LABEL:Eq_SS_reconstruct2) below.,
$\tau\mapsto\max(\,\widehat{T}(t))-\min(\,\widehat{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^{\sharp},\tau_{c})$ only from the
branching family of UPOs emanating from $T_{+}$, i.e. located in the lobe
encircling $T_{+}$ of this homoclinic orbit.

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^{\sharp}\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^{\sharp}$ leads to new UPOs
that encompass the homoclinic orbit, with amplitude that grows as $\tau$
approaches $\tau^{\sharp}$; 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.

After this jump, as $\tau$ is further decreased from $\tau^{\sharp}$, the UPOs
encompass the homoclinic orbit (blue curves in Fig. 4B), and their amplitude
keeps increasing, until eventually, loosing their instability into a Saddle-
Node bifurcation of periodic Orbits (SNO bifurcation) for $\tau=\tau^{\ast}$.
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^{\ast}$, 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 $\,\widehat{T}(t)$ to the DDE (4.4)
emanating from $T_{\textrm{init}}(s)$ ($-\tau\leq s\leq 0$), is approximated
by the solution, $\bm{y}(t)=(y_{1}(t),\cdots,y_{N}(t))^{T}$, to the GK system
Eq. (4.10) that emanates from the projection of $T_{\textrm{init}}$ onto the
Koornwinder polynomials $\mathcal{K}_{n}^{\tau}$ (see (2.8)), according to the
formula

(4.20) $\,\widehat{T}(t)\approx\sum_{j=1}^{N}y_{j}(t),$

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 $\,\widehat{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)\bm{e}_{1}+x_{2}(t)\bm{e}_{2})$, for
$n\geq 3$.

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^{\ast}$ given by
(LABEL:Eq_SS_reconstruct2), 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 (LABEL:Eq_SS_reconstruct2), 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 (LABEL:Eq_SS_reconstruct2), in order to achieve such
approximation skills.

We consider thus the following approximation formula of $\,\widehat{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) $\displaystyle
T^{*}(t)=\sum_{j=1}^{N}\Big{(}x_{1}(t)\bm{e}_{1}^{j}+x_{2}(t)\bm{e}_{2}^{j}+\sum_{n=3}^{N}\Psi_{n,\tau}\big{(}\bm{x}(t)\big{)}\bm{e}_{n}^{j}\Big{)},$
$\displaystyle\mbox{ with
}\bm{x}(t)=x_{1}(t)\bm{e}_{1}+x_{2}(t)\bm{e}_{2}\mbox{ solving
Eq.~{}\eqref{Eq_EffectiveReduced_ENSO}},$

and where $\bm{e}_{\ell}^{j}$ denotes the $j^{\mathrm{th}}$ component of the
eigenmodes $\bm{e}_{\ell}$ for $1\leq\ell\leq 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^{\sharp}$
where the homoclinic orbit emerges ($\tau=1.6$, Fig. 5B), and (iii) $\tau$ is
close to $\tau^{\ast}$ 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 (LABEL:Eq_SS_reconstruct2)
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
$\,\widehat{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)
$\varphi:(z_{1},z_{2})\mapsto\|\Psi_{\tau}(\bm{z})\|_{\mathfrak{s}},\;\;\bm{z}=z_{1}\bm{e}_{1}+z_{2}\bm{e}_{2},$

as a function of $\mathrm{Re}(z_{1})$ and $\mathrm{Im}(z_{1})$, where

(4.23)
$\|\bm{w}\|_{\mathfrak{s}}=\sqrt{\sum_{j=m_{c}+1}^{N}\big{|}\langle\bm{w},\bm{e}_{j}^{*}\rangle\big{|}^{2}},\mbox{
for any }\bm{w}\in H_{\mathfrak{s}}.$

The representation of $\varphi$ as a function of $\mathrm{Re}z_{1}$ and
$\mathrm{Im}z_{1}$ is made possible since $H_{\mathfrak{c}}$ is spanned by
eigenvectors that form a complex conjugate pair. The resulting surface,
$(\mathrm{Re}(z_{1}),\mathrm{Im}(z_{1}))\mapsto\varphi(\mathrm{Re}(z_{1}),\mathrm{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
$\mathrm{Re}(z_{1})\equiv 0$, and that corresponding to
$\mathrm{Im}(z_{1})=a\,\mathrm{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 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 $\mathrm{Re}z_{1}\equiv 0$
(cross section A), and that corresponding to
$\mathrm{Im}z_{1}=a\,\mathrm{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.

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_{\mathfrak{s}}(t),$

where $\Gamma_{j}$ denotes the projection of stable periodic solution
$\Gamma(t)$ to the DDE (4.4) onto $\bm{e}_{j}$ ($j=1,2$), while
$\Gamma_{\mathfrak{s}}(t)$ denotes the projection of $P(t)$ onto
$H_{\mathfrak{s}}=\textrm{span}\\{\bm{e}_{3},\cdots,\bm{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_{\mathfrak{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}\colon
H_{\mathfrak{c}}\rightarrow H_{\mathfrak{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)
$\frac{\,\mathrm{d}\bm{y}}{\,\mathrm{d}t}=A(\tau)\bm{y}+G_{2}(\bm{y})+G_{3}(\bm{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_{\mathfrak{c}}$ and
$H_{\mathfrak{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 $\bm{y}(t)$
of (4.24) evolving in $\mathbb{R}^{N}$, denote by
$\bm{y}_{\mathfrak{c}}=\Pi_{\mathfrak{c}}\bm{y}=\sum_{j=1}^{m_{c}}z_{j}\bm{e}_{j}$
the low-mode projection of $\bm{y}$ and by
$\bm{y}_{\mathfrak{s}}=\Pi_{\mathfrak{s}}\bm{y}=\sum_{j=m_{c}+1}^{N}z_{j}\bm{e}_{j}$
the high-mode projection of $\bm{y}$. Note that the vector
$\bm{z}=(z_{1},\ldots,z_{N})^{T}$ is nothing else than the vector $\bm{y}$
under the eigenbasis of $A(\tau)$. Let $P=[\bm{e}_{1},\ldots,\bm{e}_{N}]$ be
the $N\times N$ matrix whose columns consist of the eigenvectors of $A(\tau)$.
We have $\bm{y}=P\bm{z}$.

We introduce also $P_{\mathfrak{c}}$ as the matrix made of the first $m_{c}$
columns of $P$, namely $P_{\mathfrak{c}}=[\bm{e}_{1},\ldots,\bm{e}_{m_{c}}]$.
Similarly, $P_{\mathfrak{s}}=[\bm{e}_{m_{c}+1},\ldots,\bm{e}_{N}]$,
$\bm{z}_{\mathfrak{c}}=(z_{1},\ldots,z_{m_{c}})^{T}$ and
$\bm{z}_{\mathfrak{s}}=(z_{m_{c}+1},\ldots,z_{N})^{T}$. We have then
$\bm{y}_{\mathfrak{c}}=P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}$ and
$\bm{y}_{\mathfrak{s}}=P_{\mathfrak{s}}\bm{z}_{\mathfrak{s}}$. Finally, let
$Q$ be the $N\times N$ matrix whose columns consist of the eigenvectors of
$A^{*}$. We have then

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

where $Q^{*}$ denotes the complex conjugate transpose. Similarly as for
$P_{\mathfrak{c}}$, the matrix $Q_{\mathfrak{c}}$ denotes the matrix made of
the first $m_{c}$ columns of $Q$.

######  Proposition 4.1.

Consider $\bm{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}\colon H_{\mathfrak{c}}\rightarrow H_{\mathfrak{s}}$ be a
parameterization of $\bm{y}_{\mathfrak{s}}$ in terms of
$\bm{y}_{\mathfrak{c}}$. Assume that $\bm{y}(t)$ satisfies

(4.26) $\|\bm{y}(t)\|\leq
R,\qquad\|\bm{y}_{\mathfrak{c}}(t)+\varphi_{\tau}(\bm{y}_{\mathfrak{c}}(t))\|\leq
R,\qquad t\in[0,t_{m}].$

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

(4.27)
$\overline{\left\|\dot{\bm{z}}_{\mathfrak{c}}-\Lambda_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}-Q_{\mathfrak{c}}^{*}G_{2}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}))\Big{)}\right\|^{2}}\leq
C(R)\overline{\|P_{\mathfrak{s}}\bm{z}_{\mathfrak{s}}-\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\|^{2}},$

where $\overline{(\cdot)}$ denotes the time average over $(0,t_{m})$ and
$\Lambda_{\mathfrak{c}}$ denotes the $m_{c}\times m_{c}$ diagonal matrix made
of the eigenvalues $\lambda_{1},\ldots,\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_{\mathfrak{c}}$. Indeed, note that the GK system (4.24) can
be rewritten as:

(4.28)
$\dot{\bm{z}}=\Lambda(\tau)\bm{z}+Q^{*}\left(G_{2}(P\bm{z})+G_{3}(P\bm{z})\right),$

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

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

(4.29)
$\dot{\bm{z}_{\mathfrak{c}}}=\Lambda_{\mathfrak{c}}(\tau)\bm{z}_{\mathfrak{c}}+Q_{\mathfrak{c}}^{*}\left(G_{2}(P\bm{z})+G_{3}(P\bm{z})\right).$

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

(4.30)
$\dot{\bm{w}}=\Lambda_{\mathfrak{c}}(\tau)\bm{w}+Q_{\mathfrak{c}}^{*}\big{(}G_{2}(P_{\mathfrak{c}}\bm{w}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{w}))+G_{3}(P_{\mathfrak{c}}\bm{w}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{w}))\big{)}.$

We define the model error as

(4.31) $\displaystyle
R_{T}(\bm{z}_{\mathfrak{c}},\varphi_{\tau})=\int_{0}^{t_{m}}\left\|\frac{\,\mathrm{d}\bm{z}_{\mathfrak{c}}}{\,\mathrm{d}t}-\Lambda_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}-Q_{\mathfrak{c}}^{*}\Big{(}G_{2}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}+G_{3}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}\Big{)}\right\|^{2}\,\mathrm{d}t$
$\displaystyle=\int_{0}^{t_{m}}\left\|Q_{\mathfrak{c}}^{*}\Big{(}G_{2}(\bm{y})+G_{3}(\bm{y})-G_{2}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}-G_{3}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}\Big{)}\right\|^{2}\,\mathrm{d}t$
$\displaystyle\leq
C\int_{0}^{t_{m}}\left\|G_{2}(\bm{y})+G_{3}(\bm{y})-G_{2}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}-G_{3}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}\right\|^{2}\,\mathrm{d}t.$

Let us introduce the notation

(4.32)
$\bm{w}_{\mathfrak{s}}=\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}),$

and recall also that
$\bm{y}_{\mathfrak{c}}=P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}$ and
$\bm{y}_{\mathfrak{s}}=P_{\mathfrak{s}}\bm{z}_{\mathfrak{s}}$. Then, by the
definitions of $G_{2}$ and $G_{3}$ given in (4.8), we have

(4.33) $\displaystyle
G_{2}(\bm{y})-G_{2}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}$
$\displaystyle=-3T_{+}\Big{(}\Big{(}\sum_{n=1}^{N}y_{\mathfrak{c},n}+\sum_{n=1}^{N}y_{\mathfrak{s},n}\Big{)}^{2}-\Big{(}\sum_{n=1}^{N}y_{\mathfrak{c},n}+\sum_{n=1}^{N}w_{\mathfrak{s},n}\Big{)}^{2}\Big{)}\bm{\nu}_{N},$
$\displaystyle=-3T_{+}\Big{(}\sum_{n=1}^{N}(2y_{\mathfrak{c},n}+y_{\mathfrak{s},n}+w_{\mathfrak{s},n})\Big{)}\Big{(}\sum_{n=1}^{N}\big{(}y_{\mathfrak{s},n}-w_{\mathfrak{s},n}\big{)}\Big{)}\bm{\nu}_{N},$

and

(4.34) $\displaystyle
G_{3}(\bm{y})-G_{3}\Big{(}P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}}+\varphi_{\tau}(P_{\mathfrak{c}}\bm{z}_{\mathfrak{c}})\Big{)}$
$\displaystyle=-\Big{(}\Big{(}\sum_{n=1}^{N}y_{\mathfrak{c},n}+\sum_{n=1}^{N}y_{\mathfrak{s},n}\Big{)}^{3}-\Big{(}\sum_{n=1}^{N}y_{\mathfrak{c},n}+\sum_{n=1}^{N}w_{\mathfrak{s},n}\Big{)}^{3}\Big{)}\bm{\nu}_{N}$
$\displaystyle=-\mathcal{N}(\bm{y},\bm{w}_{\mathfrak{s}})\Big{(}\sum_{n=1}^{N}\big{(}y_{\mathfrak{s},n}-w_{\mathfrak{s},n}\big{)}\Big{)}\bm{\nu}_{N},$

where

$\displaystyle\mathcal{N}(\bm{y},\bm{w}_{\mathfrak{s}})$
$\displaystyle=3\Big{(}\sum_{n=1}^{N}y_{\mathfrak{c},n}\Big{)}^{2}+\Big{(}\sum_{n=1}^{N}y_{\mathfrak{s},n}\Big{)}^{2}+\Big{(}\sum_{n=1}^{N}w_{\mathfrak{s},n}\Big{)}^{2}$
$\displaystyle\qquad+3\Big{(}\sum_{n=1}^{N}y_{\mathfrak{c},n}\Big{)}\Big{(}\sum_{n=1}^{N}(y_{\mathfrak{s},n}+w_{\mathfrak{s},n})\Big{)}+\Big{(}\sum_{n=1}^{N}y_{\mathfrak{s},n}\Big{)}\Big{(}\sum_{n=1}^{N}w_{\mathfrak{s},n}\Big{)}.$

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

$R_{T}(\bm{z}_{\mathfrak{c}},\varphi_{\tau})\leq
C(R)\int_{0}^{t_{m}}\left|\sum_{n=1}^{N}\big{(}y_{\mathfrak{s},n}-w_{\mathfrak{s},n}\big{)}\right|^{2}\,\mathrm{d}t\leq
C(R)\int_{0}^{t_{m}}\|\bm{y}_{\mathfrak{s}}(t)-\bm{w}_{\mathfrak{s}}(t)\|^{2}\,\mathrm{d}t.\hskip
71.13188pt\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^{\ast},\tau^{\sharp}]$, where
$\tau^{\ast}$ (resp. $\tau^{\sharp}$) 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)
$\displaystyle\,\mathrm{d}\theta(t)=\Big{(}a\theta(t)-\alpha\theta(t-\tau(t))-b\theta^{2}(t)-\theta^{3}(t)\Big{)}\,\mathrm{d}t+\frac{\sigma}{1+|\theta^{2}(t)|}\,\mathrm{d}W_{t},$
(5.1b) $\displaystyle\mbox{with }\tau(t)=\tau_{0}+\epsilon
t,\qquad\tau(t)\in[\tau_{0},\tau_{1}],$

where $a=1-3T_{+}^{2}$, $b=3T_{+}$, and $\epsilon,\sigma\geq 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.333The 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. 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^{\ast}$ and the homoclinic critical value
$\tau^{\sharp}$; see caption of Fig. 4 for numerical values of $\tau^{\ast}$
and $\tau^{\sharp}$. 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 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). 444In the Stratonovich sense ([Gar04, Section 4.3.6]),
an additional term $\sigma^{2}\theta(t)(1+\theta^{2}(t))^{-3}\,\mathrm{d}t$ 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}$. 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
“fingerprints” 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 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^{\sharp}$ ($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.

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)555as
calculated from our 2D reduced equation; see Secns. 4.3 and 4.4. encountered
as $\tau$ evolves from $\tau_{0}$ to $\tau_{1}$ according 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 the 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)666To
avoid multiple notations, we still denote by $T_{-}$ the steady state to Eq.
(4.4) that corresponds to $T_{-}$ given by (4.2). 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 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.

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
$\tau_{0}$ in Eq. (5.1b) and $\sigma$ in Eq. (5.1a). Choosing $\tau_{0}$ too
far below $\tau^{*}$ and/or $\sigma$ 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 $\tau$ 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
$\tau$-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 $\delta
t=2\times 10^{-3}$, in which $\tau(t)$ is now allowed to oscillate in a
periodic fashion between $\tau_{0}$ and $\tau_{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\Delta/\tau$ with $\Delta=349,$ and $\tau=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.

Figure 9. DAH power spectrum and decadal variability. The decadal peak located
at $\approx$ 21 yr is marked by the vertical grey dashed line in panel A. The
ENSO peak is located around $\approx 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.

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 $\tau$ close
to $\tau^{\sharp}$ 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 (LABEL:Eq_SS_reconstruct2). 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
$\tau=1.5827$ in the deterministic model) provides a good synchronization with
the temporal (modulated) patterns exhibited by the decadal variability
extracted from the DAHD (black curve in Fig. 10-lower panel). Thus, one can
infer that the homoclinic bifurcation occurring at $\tau=\tau^{\sharp}$ 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 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$.

##  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)\,\mathrm{d}s$ 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 $\mathcal{P}_{N}$ and $\mathcal{Q}_{N}$ in Eq.
(4.10)

The calculations of $\mathcal{P}_{N}$ and $\mathcal{Q}_{N}$ are based on the
RHS of (4.7). $Q_{N}$ depends on the model parameters, while $\mathcal{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) $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 $\mathcal{P}_{N}$ which results from the collection of
the terms
$\sum_{k=0}^{n-1}a_{n,k}\left(\delta_{j,k}\|\mathcal{K}_{j}\|^{2}_{\mathcal{E}}-1\right)$,
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)
$\frac{\,\mathrm{d}K_{n}}{\,\mathrm{d}s}(s)=\sum_{k=0}^{n-1}a_{n,k}K_{k}(s),\quad
s\in(-1,1),$

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

$\bm{a}_{n}=(a_{n,0},\cdots,a_{n,n-1})^{T},$

solves the following upper triangular system

(A.3) $\mathbf{T}\bm{a}_{n}=\bm{b}_{n},$

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

(A.4) $\displaystyle\mathbf{T}_{i,j}$ $\displaystyle=\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 }\quad 0\leq i,j\leq n-1,$ $\displaystyle
b_{n,i}$ $\displaystyle=\begin{cases}-\frac{1}{2}(2i+1)(n+i+1)(n-i),&\text{ if
$n+i$ is even},\vspace{1em}\\\ (n^{2}+n)(2i+1)-\frac{i}{2}(n+i)(n-i+1)&\\\
\hskip 20.00003pt-\frac{1}{2}(i+1)(n-i-1)(n+i+2),&\text{ if $n+i$ is
odd}.\end{cases}$

Finally, the rescaled Koornwinder polynomials satisfy

(A.5)
$\frac{\,\mathrm{d}K^{\tau}_{n}}{\,\mathrm{d}\theta}(\theta)=\frac{2}{\tau}\sum_{k=0}^{n-1}a_{n,k}K^{\tau}_{k}(\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] by same author, _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. Roose, _Numerical bifurcation analysis of delay differential equations using DDE-BIFTOOL_ , ACM Transactions on Mathematical Software 28 (2002), 1–21.
  * [EvP04] T. Eirola 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] by same author, _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] by same author, _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+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+14] J. Sieber, K. Engelborghs, T. Luzyanina, G. Samaey, and D. Roose, _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.