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

Conceptual delay models have played a key role in the understanding of El Ni\~no-Southern Oscillation (ENSO) variability. Based on such delay models, we propose 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.

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 τ 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 τ 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 leadingorder (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 τ .These formulas are then applied to the Suarez and Schopf model in Sec. 4 to derive an effective 2D reduced GK system with τ -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.

Galerkin-Koornwinder (GK) approximations of DDEs
We consider nonlinear scalar DDEs of the form (2.1) where a, b and c are real numbers, τ > 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 Koornwinder polynomials are known to form an orthogonal set for the following weighted inner product on [−1, 1] with a point-mass, µ( dx) = 1 2 dx + δ 1 , where δ 1 denotes the Dirac point-mass at the right endpoint x = 1; see [Koo84].In other words, the following orthogonality property holds: It is also worthwhile noting that Koornwinder polynomials augmented by the right endpoint values as follows (2.4) form an orthogonal basis of the product space endowed with the inner product: (2.5) The norm induced by this inner product is denoted by ∥ • ∥ E .The basis function K n has then its ∥ • ∥ E -norm given by [CGLW16, Prop.3.1]: (2.6) This is a useful property for calculating the GK approximations, as it will come apparent below.Finally, since the original Koornwinder basis is given on the interval [−1, 1] and the state space of a DDE such as Eq.(2.1) is made of functions defined on (−τ, 0), we will work mainly with the following rescaled version of Koornwinder polynomials.The rescaled Koornwinder polynomials K τ n are defined as follows (2.7) They form orthogonal polynomials on the interval [−τ, 0] for the L 2 -inner product on (−τ, 0) with a Dirac point-mass at the right endpoint 0.

Space-time representation and GK approximations.
It is well-known that a DDE such as Eq.(2.1) is an infinite-dimensional dynamical system [Hal88] in which a history segment has to be specified over the interval [−τ, 0) to apprehend properly the existence and computation problems of its solutions [HVL93,CZ95,BZ13].Thus, given a function of time, x, solving Eq. (2.1) one distinguishes between the history segment {x(t + θ) : θ ∈ [−τ, 0)} and the current state, x(t).Denoting by u(t, θ) the history segment, one can then rewrite the DDE (2.1) as the transport equation (2.10) subject to the nonlocal and nonlinear boundary condition This reformulation is often called the space-time representation in the literature; see Fig. 1 for an illustration.Using the rescaled Koornwinder polynomials (2.7), one can then seek for approximations of u(t, θ) solving Eqns.(2.10)-(2.11)under the form (2.12) Given the property (2.9), the current state x(t) = u(t, 0) is approximated by 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, θ) provides a rigorous approximation that converges to u(t, θ) as N → ∞.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 Ndimensional 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) Here the Kronecker symbol δ 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 ∥K n ∥ E is given by (2.6).
We rewrite now the GK system (2.14) into the following compact form: (2.15) where y = (y 0 , • • • , y N −1 ) T , and in which the terms A N (τ )y and G 2 (y, τ ) group the linear and nonlinear terms in Eq. (2.14), respectively.The entries of the N × N matrix A N (τ ) are given by [CGLW16, Eq. (5.20)] (2.16) Here, the coefficients a j,k are independent of the DDE model and are determined by [CGLW16,Prop. 5 (2.17) 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.

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 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) dropping the dependence on N .We assume that nonlinear mapping, G : R N → R N , satisfies the following tangency condition The reduction formulas presented below, are derived for the case where G does not depend on τ .This choice is made to simplify the notations.The general case where G depends on τ (see (2.15)) leads to the same type of formulas as long as the tangency condition (3.2) is satisfied for all τ .Assuming that G is sufficiently smooth, then G(y) admits the following Taylor expansion for y near the origin: where denotes a homogenous polynomial of order k ≥ 2. Often, G k (y) is used hereafter as a compact notation for G k (y, • • • , y).We label the elements of the spectrum σ(A(τ )) of A(τ ) according to the lexicographical order.According to this rearrangement, the eigenvalues are labeled by a single positive integer n, so that with, for any 1 ≤ n < n ′ , either 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(τ ) experiences the following change at a critical parameter τ c : for some m c > 0, and for τ in some neighborhood U of τ c .Of course, the PES holds also when one crosses τ c from above with Re λ j (τ ) < 0 for τ > τ c , while Re λ j (τ ) > 0 when τ < τ 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 σ(A(τ )): The PES condition prevents other eigenvalues in σ s (A(τ )) from crossing the imaginary axis as τ varies in U. Hence, no eigenvalues other than those of σ c (A(τ )) change sign in the neighborhood U. Furthermore, the PES condition implies, by a continuity argument, the following uniform spectral gap by possibly reducing U accordingly [CLW15a, Lemma 6.1], (3.10) where k is the leading order of G, and Re(λ j (τ )).
To the modes that lose their stability according to the PES, we associate the reduced state space, H c , given by (3.11) while a mode e n with n ≥ m c + 1 denotes a stable mode.Throughout this article we chose not to make explicit the τ -dependence of the eigenmodes but this dependence should be kept in mind.We denote by H s the subspace spanned by these stable modes.The projector onto the subspace H s (resp. , and y c (resp.y s ) denotes the vector in H c (resp.H s ) of the low-mode amplitudes (resp.stable-mode amplitudes).The inner product in C N is denoted by ⟨•, •⟩ and is defined by In what follows we also denote by e * j the eigenmodes of the conjugate transpose A(τ ) * of A(τ ).Condition (3.10) ensures in particular that the following spectral gap holds for τ in U It is well known that the existence of a (local) exact parameterization or say in other words, of a local m cdimensional invariant manifold, is subject to the following spectral gap condition: where Lip(G| V ) denotes the Lipschitz constant of the nonlinearity G restricted to some neighborhood V of the origin in C N , and C > 0 is typically independent on V. Due to the tangency condition (3.2), condition (3.13) always holds once 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 V.
In the context of e.g.nonlinear oscillations that bifurcate from a steady state, the local invariant manifold provides an exact parameterization1 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.
1As provided for instance by a center manifold or the unstable manifold of the origin.
Theorem 3.1.Assume that G and A(τ ) satisfy the assumptions recalled above, and that the PES condition (3.8) is satisfied.Then for each τ in a neighborhood U of τ c , Eq. (3.1) admits a local invariant manifold, M τ = graph(h τ ), with h τ that maps H c into the stable subspace H s such that h τ (0) = 0 and Dh τ (0) = 0 .
Assume that the following non-resonance condition is satisfied: where G k denotes the leading-order term in the Taylor expansion of G.
ds, is well defined and is a solution to the following homological equation where L A denotes the differential operator acting on smooth mappings ψ from H c into H s , defined as: Moreover, Φ τ (X) provides the leading-order approximation of the invariant manifold function h τ characterizing M τ , in the sense that as long as X lies in a neighborhood N (τ ) of the origin in H c spanned by the m c eigenmodes e j losing their stability (see (3.11)), as τ crosses τ c .Finally, Φ τ (X) possesses the following explicit formula: , and the G n j 1 •••j k denoting the coefficients accounting for the magnitude carried out by e * n , of the nonlinear interactions through the leading-order term G k , between the low modes e j 1 , • • • , e j k (in H c ), namely: Conditions similar to (NR) arise in the smooth linearization of dynamical systems near an equilibrium [Sel85].Here, condition (NR) implies that the eigenvalues of the stable part satisfy a Sternberg condition of order k [Sel85] with respect to the eigenvalues associated with the modes spanning the reduced state space H c .This theorem is essentially a consequence of [CLM20, Theorems 1 and 2], in which condition (NR) is a stronger version of that used for [CLM20, Theorem 2]; see also [CLM20,Remark 1 (iv)].This condition is necessary and sufficient here for 0 −∞ e −sAs(τ ) Π s G k (e sAc(τ ) X) ds 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 τ in U, the following reduction principle holds: (i) Any solution y(t) to Eq. (3.1) such that y(t 0 ) belongs to M τ for some t 0 , stays on M τ over an interval of time [t 0 , t 0 + σ), σ > 0, i.e. (3.21) where y c (t) denotes the projection of y(t) onto the subspace H c .(ii) If there exists a trajectory t → y(t) such that y c (t) belongs, for all −∞ < t < ∞, to the neighborhood B(τ ) (in H c ) over which M τ is well defined, then the trajectory must lie on M τ .
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 Σ to Eq. (3.1) of any type, e.g., equilibria, periodic orbits, invariant tori, must lie in which in turn characterizes the solution y(t) in Σ, since the slaving relationship y s (t) = h τ (y c (t)) holds for any solution y(t) that belongs to an invariant set Σ for which Π c Σ ⊂ B(τ ).
Based on the approximation Theorem 3.1 and the reduction principle, it is thus reasonable to anticipate that the following effective reduced equation, provides an approximation of the invariant set Σ for which Π c Σ ⊂ B(τ ).Since Eq. (3.23) is derived in practice from a GK system, it will be referred to as an effective reduced GK system or mD reduced GK system, where m = dim(H c ). Eq. (3.23) is thus particularly suited for effective approximations of invariant sets to Eq. (3.1), made of solutions of sufficiently small amplitudes, and that bifurcate when one crosses a critical parameter τ c that destabilizes an equilibrium point.Due to the rigorous convergence results as N → ∞ 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 y(t) (and thus of 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 τ close to a critical parameter value τ 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 τ 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 higherorder approximations of an invariant manifold.The invariance equation is insightful in that respect.It is the equation satisfied by the local invariant manifold, h τ , to Eq. (3.1) in a neighborhood of the origin, namely in which we have dropped the dependence on τ 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 c contains typically a large fraction of the solution's total energy.We illustrate this idea to the case where G(y) = G 2 (y, y) + G 3 (y, y, 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 Φ τ + ψ, where ψ 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 Φ τ .By identifying the terms of order three, we obtain the following equation for ψ: where Note that the RHS is a homogeneous cubic polynomial in the X-variable.As recalled above, close to criticality, a large fraction of the energy is contained in the low modes and therefore the energy carried by y s is typically much smaller than ∥y c ∥ 2 .It is then reasonable to expect that the energy carried by Φ τ (X) is much smaller than ∥X∥ 2 for X = y c .Thus, one may assume that in the RHS of (3.25), the term Π s G 3 (X) dominates the other three terms provided that ∥G 2 (y)∥/∥y∥ 2 is on the same order of magnitude as ∥G 3 (y)∥/∥y∥ 3 .As a consequence, it is reasonable to hope for reasonably accurate approximations of ψ with h 3 by simply solving the equation: 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) 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.
where the unknown T represents the non-dimensionalized sea surface temperature (SST) anomalies at the eastern equatorial Pacific Ocean, τ and α are positive constants, and the physically relevant range of α used in [SS88] is (0, 1).For this given range of α, Eq. (4.1) admits three fixed points: 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) The above equation fits into Eq.(2.1) with By applying Eqns.(2.15)-(2.17) to Eq. (4.4), we obtain the following N -dimensional GK system (4.6) where the entries of the matrix A(τ ) are given here by (4.7) As in Sec.2.3, the a n,k are determined thanks to Proposition 1 recalled in Appendix A, and the ∥K j ∥ 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 where ν N denotes the N -dimensional column vector (4.9) We arrive finally at the following N -dimensional GK approximation for the DDE (4.4) (4.10) We refer to Appendix A for the computation of P N and 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 α = 0.75.For this value of α, the steady state is given by T + = √ 1 − α = 0.5.We conclude this section with Table 1 that shows the L ∞ -error achieved by various N -dimensional GK systems in approximating stable limit cycles to the DDE (4.4) for τ -values that will be considered later on.
Table 1.L ∞ -error achieved by the N -dimensional GK approximations to the DDE (4.4).Here, the L ∞ -errors are computed over one period of the stable limit cycle for each tested τ -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 δt = τ /2 18 .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 Ψ τ 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 Ψ τ writes in compact form

GK dimension
for which m c is determined by analyzing the modes that destabilize once a critical delay parameter τ is crossed.
In the eigenbasis coordinate system of A(τ ), the GK system (4.10)writes (4.12) where y j = ⟨y, e * j ⟩, with inner product defined in (3.12).We observe that a dominant pair of complex eigenvalues crosses the imaginary axis as τ crosses from below the critical value τ c ≈ 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 τ above the critical value; see Fig. 2. We choose thus m c = 2 and the subspaces H c and H s as follows: (4.13) in which (e 1 , e 2 ) denotes the conjugate pair of modes that destabilize as τ crosses τ c .Setting m c = 2 in (3.27), the 2D reduced GK system (4.11)(in C 2 ) writes in the eigenbasis coordinates: (4.14) ẋj = λ j (τ )x j + ⟨G(x + Ψ τ (x)), e * j ⟩, j = 1, 2, with x = x 1 e 1 + x 2 e 2 , and where the parameterization Ψ τ of the neglected modes in H 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 Ndimensional 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.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 τ = τ c , but also bifurcations taking place at other critical values of τ .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 τ is further decreased to a value τ ♯ for which the UPOs (dashed black curves in Fig. 4A) hit a saddle equilibrium from both sides, resulting into a homoclinic orbit when τ reaches a critical value τ ♯ ; see cyan curve in Fig. 4A.The presence of this homoclinic orbit is manifested by the jump displayed by the curve of UPO's amplitudes2, τ → max( T (t)) − min( 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 (τ ♯ , τ c ) only from the branching family of UPOs emanating from T + , i.e. located in the lobe encircling T + of this homoclinic orbit.
After this jump, as τ is further decreased from τ ♯ , 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 τ = τ * .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 τ is increased from τ * , 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.

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 τ 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 T (t) to the DDE (4.4) emanating from T init (s) (−τ ≤ s ≤ 0), is approximated by the solution, y(t) = (y 1 (t), • • • , y N (t)) T , to the GK system Eq.(4.10) that emanates from the projection of T init onto the Koornwinder polynomials K τ n (see (2.8)), according to the formula (4.20) 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 ∞ is small in (4.20) (see Table 1), we can expect to reach a good approximation of T (t), at least for τ sufficiently close to τ 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 Ψ n,τ (x 1 (t)e 1 + x 2 (t)e 2 ), for n ≥ 3.
We consider thus the following approximation formula of T (t) that is built solely from the solutions of the 2D reduced GK system Eq.(4.14) and the high-mode parameterization Ψ τ given by (3.27): (4.21) with x(t) = x 1 (t)e 1 + x 2 (t)e 2 solving Eq. (4.14), and where e j ℓ denotes the j th component of the eigenmodes e ℓ for 1 ≤ ℓ ≤ N .Recall that the latter depend on τ and N as eigenvectors of A(τ ) 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) τ is close to τ c from below (τ = 1.7,Fig. 5A), (ii) τ is close (from below) to τ ♯ where the homoclinic orbit emerges (τ = 1.6, Fig. 5B), and (iii) τ is close to τ * at which the SNO bifurcation occurs (τ = 1.562,Fig. 5C).The last case shown corresponds τ away from τ c from above (τ = 1.9, Fig. 5D) for which only a stable limit cycle exists as periodic orbit.In each case, the approximation formula (4.21) based on the 2D reduced GK system (4.14) and the parameterization Ψ τ 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 τ 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 Ψ τ built from the GK system (4.10),already for N = 6.This property remains valid even for solution's amplitude (for T (t)) of order one such as shown for τ = 1.9 in Fig. 6, which corresponds to the furthest value τ to τ c examined in Fig. 5.These comments are made visual in Fig. 6.This figure shows the graph of the norm of Ψ τ , as a function of Re(z 1 ) and Im(z 1 ), where (4.23) ⟨w, e * j ⟩ 2 , for any w ∈ H s .
The representation of φ as a function of Rez 1 and Imz 1 is made possible since H c is spanned by eigenvectors that form a complex conjugate pair.The resulting surface, (Re(z 1 ), Im(z 1 )) → φ(Re(z 1 ), Im(z 1 )), is then intercepted by two vertical planes for a better visual inspection of the distance between the stable periodic trajectory from Eq. (4.10) and this surface.The left and right panels of Fig. 6 show the curves (in red) obtained as cross sections with, respectively, the vertical plane corresponding to Re(z 1 ) ≡ 0, and that corresponding to Im(z 1 ) = a Re(z 1 ) + b, with a = −0.63 and b = 0.21.The intersections between these planes and the stable periodic solution Γ(t) to the DDE (4.4) are shown by the black dots shown in each panel of Fig. 6.The curves (in red) are obtained as intersections between the graph of φ given by (4.22) (for N = 6) with respectively, the vertical plane corresponding to Rez 1 ≡ 0 (cross section A), and that corresponding to Imz 1 = a Rez 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 τ = 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 where Γ j denotes the projection of stable periodic solution Γ(t) to the DDE (4.4) onto e j (j = 1, 2), while Γ s (t) denotes the projection of For the case shown in Fig. 6 (τ = 1.9), the time-averaged over one period of the ratio, ∥R(t)∥ 2 /∥Γ s (t)∥ 2 , is approximately equal to 2.5 × 10 −2 , indicating that Γ(t) is almost slaved to the manifold Ψ τ .This observation explains that the discrepancies observed in Fig. 6 between the DDE solution and the manifold Ψ τ 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.

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 Ψ τ , 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 Let us first recall that the N -dimensional GK system (4.6) of the Suarez-Schopf DDE (4.4) is given by (4.24) with the matrix A(τ ) and the nonlinear terms given by (4.7) and (4.8), respectively.As in Sec. 3, we assume that A(τ ) is diagonalizable over C, which is observed to be true for all the numerical results reported in this study for this DDE model.Recall that H c and H s are the eigen subspaces spanned by the low modes and the high modes, respectively.We will also make use of the following notations.Given a solution y(t) of (4.24) evolving in R N , denote by y c = Π c y = mc j=1 z j e j the low-mode projection of y and by y s = Π s y = N j=mc+1 z j e j the high-mode projection of y.Note that the vector z = (z 1 , . . ., z N ) T is nothing else than the vector y under the eigenbasis of A(τ ).Let P = [e 1 , . . ., e N ] be the N × N matrix whose columns consist of the eigenvectors of A(τ ).We have y = P z.
We introduce also P c as the matrix made of the first m c columns of P , namely P c = [e 1 , . . ., e mc ].Similarly, P s = [e mc+1 , . . ., e N ], z c = (z 1 , . . ., z mc ) T and z s = (z mc+1 , . . ., z N ) T .We have then y c = P c z c and y s = P s z s .Finally, let Q be the N × N matrix whose columns consist of the eigenvectors of A * .We have then (4.25) where Q * denotes the complex conjugate transpose.Similarly as for P c , the matrix Q c denotes the matrix made of the first m c columns of Q. Proposition 4.1.Consider y(t) to be a solution to Eq. (4.24) over a given interval [0, t m ].Assume that A(τ ) is diagonalizable over C. Let φ τ : H c → H s be a parameterization of y s in terms of y c .Assume that y(t) satisfies Then there exists a constant C(R) > 0 such that where (•) denotes the time average over (0, t m ) and Λ c denotes the m c × m c diagonal matrix made of the eigenvalues λ 1 , . . ., λ mc of A(τ ).
Proof.The desired result (4.27) follows from a direct calculation by comparing the vector field of the reduced system associated with φ τ (see (4.30) below) and the projection of the GK vector field in (4.24) onto H c .Indeed, note that the GK system (4.24) can be rewritten as: where Λ(τ ) = Q * A(τ )P .By projecting the above equation onto the low modes, we have ) .On the other hand, the reduced system associated with φ τ is (4.30) We define the model error as Let us introduce the notation (4.32) and recall also that y c = P c z c and y s = P s z s .Then, by the definitions of G 2 and G 3 given in (4.8), we have Finally, by using (4.33) and (4.34) in (4.31), we conclude, in virtue of assumption (4.26), that

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 τ -interval [τ 0 , τ 1 ] containing [τ * , τ ♯ ], where τ * (resp.τ ♯ ) corresponds to the SNO (resp.homoclinic) bifurcation point.Since these invariant sets are τ -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 θ(t) = T − T + : where a = 1 − 3T 2 + , b = 3T + , and ϵ, σ ≥ 0 while W t denotes a Brownian motion.The noise term is a Lorentzian function favouring larger noise when θ ≈ 0 while damping the noise effects when θ 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.3As 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 α = 0.75, σ = 0.2, with τ evolving according to Eq. (5.1b) over an interval [0, t m ], with τ 0 = 1.45, ϵ = 8.4 × 10 −4 , and t m = 237.8 in time unit of the model.With these parameters, we have τ (t m ) = τ 1 = 1.65, and thus the interval [τ 0 , τ 1 ], within which τ (t) varies, contains both the SNO critical value τ * and the homoclinic critical value τ ♯ ; see caption of Fig. 4 for numerical values of τ * and τ ♯ .The rationale behind this numerical setup is guided by a simple intuition.Indeed, as the delay parameter τ slowly drifts across a portion of the bifurcation diagram shown in Fig. 4 and when the noise strength parameter σ 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 3The 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.
intervals preceding the crossing of the SNO critical value τ * , the solution's variability is affected by invariant sets such as limit cycles that emerge only for τ > τ * .
We numerically observed that this ability of the stochastic model to "feel" the deterministic model's dynamics at a given τ -value prior τ (t) crosses this value, is actually independent of how the noise is interpreted (Itô vs. Stratonovich).4 Thus, noise can e.g.excites parameter regimes with oscillatory dynamics as τ (t) approaches τ -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).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)5 encountered as τ evolves from τ 0 to τ 1 according 4In the Stratonovich sense ([Gar04, Section 4.3.6]),an additional term σ 2 θ(t)(1 + θ 2 (t)) −3 dt is added to Eq. (5.1a).In the present parameter setting, we have checked that this extra term has little influence on the dynamics as it scales here with σ 2 .5as calculated from our 2D reduced equation; see Secns.4.3 and 4.4.
to (5.1b): a UPO at τ = 1.5607 encompassing the homoclinic orbit, two UPOs symmetric with respect to the saddle point (for τ = τ 1 = 1.65), and two stable limit cycles, shown as two closed black (plain) curves whose the one with smaller diameter corresponds to τ = 1.65 while the other corresponds to τ = 2.The snippets marked by the turquoise and blue consecutive segments of the TSP as shown in Fig. 7B, correspond to a quiet episode (turquoise) followed by a large excursion (blue) reminiscent with what has been observed for certain El Niño events as recorded by the Niño 3.4 SST index.Such a large excursion following a certain stillness evokes the major El Niño event that occurred in Dec-Jan 2016 [LTW + 17].It is worthwhile noting that the Niño 3.4 SST index during the four years preceding this event exhibited indeed temporal patterns resembling the episode marked in turquoise; cf.e.g. the global reanalysis product of Niño 3.4 index as documented by the E.U.Copernicus Marine Service Information [CME18].The phase portrait reveals that these consecutive events correspond, in our model, to a meandering of the solution path near what corresponds to a La Niña steady state T − (turquoise path in Fig. 7A)6 followed by a nearly cyclic orbit (blue path in Fig. 7A).
The snippets shown in yellow and red in Fig. 7B emphasize a different mechanism leading to an El Niño event.Indeed, by inspecting the corresponding path in Fig. 7A, we observe that while the El Niño excursion corresponds still to a nearly cyclic orbit (red), its ignition is now preceded by a path's wandering around an UPO surrounding itself the La Niña steady state.Due to their behavior starting either from the La Niña steady state or an UPO surrounding it, before engaging into a quasi cyclic excursion, we call such episodes as Transition Snippets.
Such transition snippets provide the generic building blocks at work for the production of El Niño-like patterns displayed by the solution to Eq. (5.1) such as shown in Fig. 7B.Our analysis reveals thus that a solution path wandering around an UPO surrounding the La Niña steady state or having a closer meandering around it can be interpreted as an early warning signal [SBB + 09] for an El Niño-like pattern to occur in this model.
Other episodes displayed by the TSP's temporal patterns correspond to a path connecting the El Niño with the La Niña equilibria.Such an event is shown in Fig. 8 with color coding marking the following three stages: (i) path escape towards the El Niño equilibrium (blue path), (ii) path meandering in a neighborhood of that equilibrium (yellow path), and (iii) path escape from that neighborhood towards the La Niña equilibrium (red path).The mechanisms of fabric of the TSP's temporal patterns reported here are therefore diverse, involving various scenarios such as noise-sustained oscillations associated with the SNO critical transition and its constitutive UPOs and stable limit cycles, or transitions between steady states.Due to the presence of (nonlinear) periodic orbits these mechanisms are characteristic of more intricate interactions between noise and nonlinear effects than classically encountered in noise-induced tipping dynamics between multiple equilibria; see e.g.[GHJM98] and [HL84, Chapter 6].In particular, the good capture of such interactions by stochastic reduced models remains to be explored.We mention that stochastic invariant manifolds techniques [CLW15b,CLMW23] and their extension [CLM23] should be useful in that respect.
Finally, it worth mentioning that such mechanisms are not dependent on the noise path used to drive Eq. (5.1a) but rather conditioned to the choice of τ 0 in Eq. (5.1b) and σ in Eq. (5.1a).Choosing τ 0 too far below τ * and/or σ too small could lead to completely different solution's temporal behavior with for instance more frequent transitions between equilibria and very rare excitation of oscillatory behavior.While allowing τ in Eq. (5.1) to drift is motivated by physical considerations as mentioned above, it is not excluded though, for well-calibrated model's parameters, that stochastic solutions for frozen τ -value could share similar temporal variability attributes than those reported here.Regardless, drifting or not, the nonlinear invariant sets (UPOs, stable limit cycles, steady states) of the deterministic flow play a central role in shaping the temporal patterns exhibited by the TSP.We have been though voluntarily evasive on the role of the homoclinic orbit in the fabric of the TSP variability.The next section clarifies this point.

Decadal variability and homoclinic orbit.
To decipher the role of the homoclinic orbit in the fabric of the TSP variability, we need to quantify this variability in more precise terms.To do so, we analyze the spectral content of a TSP produced over a long-time integration of Eq. (5.1) consisting of 10 6 time steps with δt = 2 × 10 −3 , in which τ (t) is now allowed to oscillate in a periodic fashion between τ 0 and τ 1 in (5.1b).The values of these and other model's parameters are the same as in Sec.5.1.To make a physical sense of the variability displayed by the TSP, we adopt the following conversion formula [BTR07] s = t∆/τ with ∆ = 349, and τ = 1.7, in which s denotes the physical time in year, while t denotes the non-dimensional time unit of the DDE model.To analyse the frequency content, we use the data-adaptive harmonic decomposition (DAHD) method [CK17] that is a signal processing method which has been employed in many applications to analyze the temporal variability and extract coherent patterns of complex (multivariate) time series; see e.g.[KCG18,KCYG18,KC18,KCB18]. Said in simple terms, the method allows for performing a Fourier-like analysis from time-lagged correlations while extracting empirical modes of variability that are naturally ranked per Fourier frequency [CK17, Theorem V.1].A simple projection of the signal onto a selected group of such modes over a given frequency range allows in turn for calculating in the time domain the contribution of that range within the original signal.
We applied DAHD to the aforementioned TSP simulated over a time window of approximately 1, 100 yrs.The DAH power spectrum reveals two distinct broadband peaks: a 21-yr peak associated with decadal variability, and a 5.5-yr peak associated with ENSO interannual variability, corresponding in Fig. 9A, to the small bump made of blue dots for decadal variability, and to the more energetic one located to its right, for the interannual one.The projection of the TSP onto the subspace spanned by the DAH modes (DAHMs) associated with the decadal variability (blue dots in Fig. 9A) is shown by the blue curve in Fig. 9B.
Decadal variability of ENSO has been documented in recent studies [DCP + 21].The origin of this decadal variability in our model (5.1) can be traced back to UPOs that are located close to the homoclinic orbit, i.e. for τ close to τ ♯ from below; see blue curves in inset B of Fig. 4. Recall that the 2D reduced GK system (4.14)allows for a simple computation of the UPOs of the DDE (4.4) by a simple backward integration of Eq. (4.14) followed by a lifting of the orbit via the formula (4.21).Using this approach, UPOs closer to the homoclinic orbit such as shown in Fig. 4B can be easily probed.Among these UPOs, we found that the UPO with periodicity of 18.23-yr (blue curve in Fig. 10-lower panel, corresponding to τ = 1.5827 in the deterministic model) provides a good synchronization with the temporal (modulated) patterns exhibited by the decadal variability extracted from the DAHD (black curve in Fig. 10-lower panel).Thus, one can infer that the homoclinic bifurcation occurring at τ = τ ♯ 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.shown here corresponds to a UPO of period 18.23-yr obtained for τ = 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 τ = 1.7689.

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 τ 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 t t−τ x(s) ds are simply approximated by τ y 0 (t) − τ 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.

Figure 1 .
Figure 1.Panel A shows a solution u(t, θ) to a transport problem of the form (2.10)-(2.11)obtained as a reformulation of the DDE ẋ = ax(t − τ ) − bx(t − τ ) 3 from [CGLW16].Panel B shows the time series u(t, 0) at the right endpoint θ = 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 τ = 0.4.

Figure 2 .
Figure2.Eigenvalues dependence as τ crosses its critical value τ c .Are shown here, the first 10 pairs (λ j (τ ), λ j (τ )) as τ is increased from τ = 1.3 to τ = 2.5, for α = 0.75.These pairs move all from left to right.The pair of eigenvalues that crosses the imaginary axis for the critical delay parameter τ = τ c is shown in red, while the stable pairs are shown in black.These eigenvalues are computed from the GK linear part A(τ ) 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), λ = (1 − 3T 2 + ) − αe −λτ , up to a maximal error of 10 −4 , as τ varies from τ = 1.3 to τ = 2.5.

2ObtainedFigure 4 .
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 Ψ τ 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 τ < τ c , and loses its stability through a subcritical Hopf bifurcation at τ c ≈ 1.7408; see Table2and Eq. (4.19).As τ approaches τ ♯ ≈ 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 τ below τ ♯ leads to new UPOs that encompass the homoclinic orbit, with amplitude that grows as τ approaches τ ♯ ; see blue dashed curves in inset B. This latter branch of UPOs terminates at τ * ≈ 1.5592 through a SNO bifurcation that gives rise to a branch of stable periodic orbits shown by the red curves in inset C.

Figure 5 .
Figure 5. Approximation results: DDE vs Effective Reduced GK system.The solutions are shown in lagged coordinates for α = 0.75.For different values of τ as indicated, the solutions to the DDE (4.4) are compared to those obtained from the formula of T * given by (4.21), built from the solutions of the 2D reduced GK system Eq.(4.14) and the parameterization Ψ τ given by (3.27).The stable (resp.unstable) DDE limit cycles are shown in black (resp.blue).The stable (resp.unstable) limit cycles, obtained from the 2D reduced system after lifting through (4.21), are shown by the orange (resp.cyan) dashed curves.The stable (resp.unstable) equilibria are shown as filled (resp.empty) circles.Here, it is sufficient to use a GK system of dimension N = 6 to build the parameterization Ψ τ and the modes and eigenvalues used in (4.21), in order to achieve such approximation skills.

Figure 6 .
Figure 6.Visualization of the high-mode parameterization Ψ τ used in the effective reduced Eq.(4.14).The curves (in red) are obtained as intersections between the graph of φ given by (4.22) (for N = 6) with respectively, the vertical plane corresponding to Rez 1 ≡ 0 (cross section A), and that corresponding to Imz 1 = a Rez 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 τ = 1.9), with these vertical planes.

Figure 7 .
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 θ(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 τ (t) = τ * (t = 130) and τ (t) = τ ♯ (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 τ -values as pointed out in the text.

Figure 8 .
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.

Figure 9 .
Figure 9. DAH power spectrum and decadal variability.The decadal peak located at ≈ 21 yr is marked by the vertical grey dashed line in panel A. The ENSO peak is located around ≈ 5.5 yr.The projection of the TSP onto the subspace spanned by the DAH modes (DAHMs) associated with the decadal variability (blue dots in panel A) is shown by the blue curve in panel B.

Figure 10 .
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 curveshown here corresponds to a UPO of period 18.23-yr obtained for τ = 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 τ = 1.7689.
[SY02,perty (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, y c (t), is sufficient to entirely determine any solution y(t) that belongs to such an invariant set; see[SY02, Theorem 71.4].Furthermore, y c (t) is obtained as the solution of the following reduced m c -dimensional problem (3.22)