Continuity, Localization, and Cosmology in Warped Geometry

Gregory Gabadadze, Daniel Older and David Pirtskhalava

Center for Cosmology and Particle Physics, Department of Physics

New York University, New York, NY, 10003, USA

###### Abstract

This is the first of two papers studying localization of massive bulk fields
on a bane in 5D anti-de Sitter spacetime, and some of their cosmological
consequences. Here we focus on a massive 5D scalar, which is known to lack a
localized mode, and discuss how a seeming discontinuity between this theory
and the massless theory – known to support a localized zero mode – is resolved
thanks to peculiar analytic properties of the massive two-point amplitude.
Furthermore, we propose a boundary term that leads to the emergence of a
massless localized mode in the massive theory. Last but not least, we consider
the case when the brane world-volume is de Sitter spacetime, and prove the
existence of a localized massive mode. We discuss how these results, taken
collectively, can be used to describe the accelerated expansion due to the
massive 5D scalar field in an early, or in a late-time universe.

##  1 Introduction and Summary

Randall and Sundrum (RS) [1] have shown that a 5D massless graviton yields a
localized massless 4D graviton in 5D AdS spacetime with a $Z_{2}$ symmetric
brane; the rest of the Kaluza-Klein(KK) modes form a gapless continuum, with
the wavefunctions for the lighter KK modes suppressed on the brane, leading to
approximate recovery of 4D physics at low energies.111Note that $AdS_{5}$
refers to the universal cover of 5D AdS space-time [2]. There exists a stable
solitonic domain wall solution with 4D Poincare invariant worldvolume, and
$Z_{2}$ symmetry in the direction transverse to the worldvolume [3]. This
solution, in a certain approximation, can be thought of as infinitely thin
brane with all its fluctuations becoming negligible for a low energy brane
observer. The RS discovery has a numerous physical and theoretical
consequences, which are well known.

It turns out that massless 5D scalar also yields a localized massless 4D mode
on the brane worldvolume [4], with the gapless KK continuum and the KK
wavefunctions identical to those of massless 5D graviton. However, massless
vectors and spinors do not produce localized massless modes in the minimal
setup.

Our goal is to consider whether massive spin-0, spin-1, and spin-2 fields can
be localized or quasilocalized on a brane in $AdS_{5}$. This is the first of
the two papers on this topic; here we focus on a 5D massive scalar field to
emphasize some of the subtleties which are independent of the spin, but also
to delineate cosmological scenarios specific to such a massive 5D scalar.

The question of localization of a 5D massive field was first studied in [5],
where it was found that there is no localized mode irrespective how small the
5D mass is. Yet, there is a resonant mode, that can mimic physics of the zero
mode in a certain approximation [5].

Seemingly, there is discontinuity between the massless and massive theory –
the former has a massless 4D mode in its spectrum while the latter does not,
irrespective how small the 5D scalar mass is. Indeed, this discontinuity is
real as long as the spectra of the modes are concerned – the bulk mass makes
the wold-be localized mode non-normalizable, for any nonzero value of the
mass. However, we will show that there is no discontinuity in the two-point
amplitude of the theory. This involves careful consideration of the analytic
properties of this amplitude, as it is done in Section 2.

In Section 3 we show how one can modify the brane worldvolume theory by adding
a simple 4D term for the scalar to reinstate the massless 4D localized mode
for a massive 5D bulk scalar field. The new term is a "tachyonic" 4D
worldvolume mass for the scalar, which depends on the bulk mass. We should
note, however, that the 4D term is overwhelmed by the positive 5D mass term
and nowhere one encounters any tachyonic instabilities in the theory. The new
term create an additional "attraction" in a potential for the spectral problem
for the KK modes, and this is enough to reinstate the massless mode.

In Section 4 we consider 5D massive scalar in the background geometry with 4D
de Sitter worldvolume. We prove the existence of a localized massive mode in
this case and a KK continuum starting above the gap determined by the
curvature of the 4D de Sitter spacetime. We then outline how such a geometry
can be an approximation to the one on which 4D curvature is provided by a
slowly rolling 5D massive scalar field. We point out differences of such a
scheme from the one with 4D massive scalar slowly rolling and providing
inflation or late time acceleration (quintessence).

##  2 Massive Scalar in 5D

### The setup

Following [5], we consider a massive scalar field in the 5-dimensional anti de
Sitter spacetime endowed with ${\Large\bf Z_{2}}$ symmetry across its
boundary, described by the following interval

$ds^{2}=\Omega^{2}(z)\,\eta_{MN}\,dx^{M}dx^{N}\,,\quad\Omega(z)=\frac{L}{L+|z|}\,,$
(1)

where the mostly minus metric convention is assumed. The coordinate $z$,
parametrizing the fifth dimension ranges from $-\infty$ to $+\infty$ and the
theory is constrained to be ${\Large\bf Z_{2}}$ invariant under the flip of
sign $z\to-z$, under which the scalar is assumed to be even,
$\phi(x,z)=\phi(x,-z)$ . We are interested in the 4D effective theory on the
physical brane, located at $z=0$, where $\phi$ is additionally coupled to a 4D
source $j(x)$ . The complete action thus reads

$\displaystyle S_{\phi}$ $\displaystyle=\int
d^{4}x\int_{-\infty}^{\infty}dz\,\bigg{[}\sqrt{-g}\,\Big{(}\,\frac{1}{2}g^{MN}\partial_{M}\phi\partial_{N}\phi-\frac{1}{2}m_{5}^{2}\phi^{2}\Big{)}+L\phi(x)j(x)\delta(z)\bigg{]}$
(2) $\displaystyle=\int
d^{4}x\int_{-\infty}^{\infty}dz\,\Omega^{3}\Big{(}\frac{1}{2}\partial_{\mu}\phi\,\partial_{\mu}\phi-\frac{1}{2}(\partial_{z}\phi)^{2}-\frac{1}{2}\Omega^{2}m_{5}^{2}\phi^{2}\Big{)}+\int
d^{4}xL\phi(x)j(x).$

All 4D indices on the last line are assumed to be contracted with the flat 4D
metric, and we do not distinguish between upper and lower 4D components.
Varying this action with respect to $\phi$ yields the following equation of
motion

$\displaystyle\Big{(}-\Box+\partial_{z}^{2}-\frac{3\text{sgn}(z)}{|z|+L}\partial_{z}-\frac{(m_{5}L)^{2}}{(|z|+L)^{2}}\Big{)}\phi(x,z)=-Lj(x)\delta(z)\,,$
(3)

which, when integrated across the brane (that is, within the interval
$z\in[-\epsilon,\epsilon]$, with $\epsilon\rightarrow 0$) implies the
following boundary condition for $\phi$:

$\displaystyle\partial_{z}\phi|_{z=0}=-\frac{L}{2}\,j(x)\,.$ (4)

In order to derive the 4D effective action on the brane, we will need to solve
the system (3) and (4). This is the subject of the next subsection.

### Kaluza-Klein modes

Consistently with the equation of motion (3) and boundary condition (4), the
5D field $\phi$ can be decomposed in terms of the 4D KK modes as follows:

$\displaystyle\phi(x,z)=\int_{0}^{\infty}dmL\,\phi^{(m)}(x)\,\chi^{(m)}_{\nu}(z)\,,$
(5)

where we have defined

$\nu=\sqrt{4+(m_{5}L)^{2}}\,.$

The KK wavefunctions $\chi^{(m)}_{\nu}$ satisfy the following bulk equation

$\displaystyle\Big{(}\partial_{z}^{2}-\frac{3\text{sgn}(z)}{|z|+L}\partial_{z}-\frac{(m_{5}L)^{2}}{(|z|+L)^{2}}\Big{)}\chi^{(m)}_{\nu}(z)=-m^{2}\,\chi^{(m)}_{\nu}(z)\,,$
(6)

complemented with the boundary condition at the brane

$\displaystyle\partial_{z}\chi^{(m)}_{\nu}(z)|_{z=0}=0\,.$ (7)

Explicitly, these KK wavefunctions read [5]

$\displaystyle\chi^{(m)}_{\nu}(z)=\sqrt{\frac{mL}{2}}\left(\frac{|z|+L}{L}\right)^{2}\Big{[}a_{m}J_{\nu}\big{(}m\left(|z|+L\right)\big{)}+b_{m}Y_{\nu}\big{(}m\left(|z|+L\right)\big{)}\Big{]}\,,$
(8)

where the two coefficients $a_{m}$ and $b_{m}$ are given by the following
expressions

$\displaystyle a_{m}$
$\displaystyle=-\frac{Y_{\nu-1}(mL)-\frac{\nu-2}{mL}Y_{\nu}(mL)}{\sqrt{\left(Y_{\nu-1}(mL)-\frac{\nu-2}{mL}Y_{\nu}(mL)\right)^{2}+\left(J_{\nu-1}(mL)-\frac{\nu-2}{mL}J_{\nu}(mL)\right)^{2}}}$
(9) $\displaystyle b_{m}$
$\displaystyle=\frac{J_{\nu-1}(mL)-\frac{\nu-2}{mL}J_{\nu}(mL)}{\sqrt{\left(Y_{\nu-1}(mL)-\frac{\nu-2}{mL}Y_{\nu}(mL)\right)^{2}+\left(J_{\nu-1}(mL)-\frac{\nu-2}{mL}J_{\nu}(mL)\right)^{2}}}\,.$

Importantly, for non-zero bulk mass $m_{5}$ there is no normalizable zero
mode, that is, no mode with $m=0$ [5]. (For $m_{5}$ strictly zero, such a mode
does exist and would correspond to the $z$-independent bulk profile
$\chi^{(0)}=\text{const}$ [4].) Nevertheless, even for a non-zero $m_{5}$, the
spectrum features a _gapless_ continuum of massive KK modes, whose
wavefunctions satisfy the usual orthonormality conditions

$\displaystyle\int_{-\infty}^{\infty}dz\,\left(\frac{L}{|z|+L}\right)^{3}\chi^{(m)}_{\nu}(z)\chi^{(m^{\prime})}_{\nu}(z)=\delta(m-m^{\prime})\qquad(m,m^{\prime}>0)\,.$
(10)

Plugging decomposition (5) into the original 5D action (2), one arrives at the
4D effective action for the KK continuum

$\displaystyle S_{\phi}=L$ $\displaystyle\int
d^{4}x\int_{0}^{\infty}dmL\,\Big{[}\frac{1}{2}(\partial\phi^{(m)})^{2}-\frac{m^{2}}{2}(\phi^{(m)})^{2}+\chi^{(m)}(0)\,\phi^{(m)}j\Big{]}\,.$
(11)

At each KK level, the canonically normalized 4D field $\phi^{(m)}/\sqrt{L}$
couples to the brane source $j$ with strength $\chi^{(m)}(0)/\sqrt{L}$. This
coupling weakens for higher KK modes.

The force due to the exchange of massive KK modes between two four-dimensional
brane sources is determined by the following amplitude (transformed to
momentum space along the brane dimensions)

$\displaystyle\Delta_{\nu}(p^{2})=\int_{0}^{\infty}dmL\,\frac{|\chi^{(m)}_{\nu}(0)|^{2}}{p^{2}-m^{2}+i\epsilon}=\frac{L^{2}}{2}\,\bigg{[}(pL)\frac{H^{(1)}_{\nu-1}(pL)}{H^{(1)}_{\nu}(pL)}-(\nu-2)\bigg{]}^{-1}\quad(\nu>2)\,,$
(12)

where $p\equiv\sqrt{p^{2}}$. This formula, _valid strictly for_ $\nu>2$, was
first derived in Ref. [5]. Addressing the case of $\nu=2$, on the other hand,
is subtle and will be our focus in the remainder of this section.

### Continuity in the bulk mass

The main task of the present subsection is demonstrating continuity of the
physical predictions of the theory in the $m_{5}\to 0$ limit. This limit is
non-trivial: as we have remarked above, depending on whether the parameter
$m_{5}$ is strictly zero or not, there is an extra zero mode present in the KK
spectrum, leading to a discontinuity in the degrees of freedom in the 4D
theory on the brane.

Ref. [5] has demonstrated how exactly continuity occurs for non-zero, but
small $m_{5}$ at the level of the Newtonian potential, experienced by brane
sources due to the exchange of the scalar’s KK modes. The crucial role is
played by a special resonant mode, composed of the KK modes of the 4D theory.
This resonance, which only exists for $m_{5}\neq 0$ (or $\nu\neq 2$), is
characterized by mass and width, given by

$\displaystyle m_{0}^{2}\approx\frac{m_{5}^{2}}{2}\,,\qquad\Gamma$
$\displaystyle\approx\frac{\pi}{16\sqrt{2}L}(m_{5}L)^{3}\,.$ (13)

(Notice that for reasonable values of the bulk mass, $m_{5}\lesssim 1/L$, the
resonance’s width is parametrically smaller than its mass.) The presence of
the pole corresponding to the resonance (13) can be readily established by
studying the analytic structure of the amplitude (12), see Ref. [5] for
details. The negative imaginary part of the pole tends to zero as $m_{5}\to
0$, or equivalently as $\nu\to 2$ (of course, in this limit the real part
vanishes as well and the pole asymptotically merges with the origin of the
complex plane). Therefore, in the limit $\nu\to 2$ one can conveniently write
Eq. (12) as

$\Delta_{\nu}(p^{2})|_{\nu\to
2}=\text{P.V.}\left(\frac{L^{2}}{2}\frac{H_{2}^{(1)}(pL)}{(pL)H^{(1)}_{1}(pL)}\right)-i\pi\delta(p^{2})\,,$
(14)

where $\text{P.V.}\left(\dots\right)$ in the first term denotes the principal
value of the expression in the parentheses, while the second term essentially
provides the $i\epsilon$-prescription for the pole at $p^{2}=0$, corresponding
to the $m_{5}\to 0$ limit of the resonance. To summarize, the analytic
structure of the expression in (14) features the pole, corresponding to the
$m_{5}\to 0$ limit of the 4D resonance scalar (responsible for the $r^{-1}$
piece in the static potential), as well as a branch cut, corresponding to the
entire gapless KK continuum (which leads to the $r^{-3}$ piece in the static
potential).

We’d like to compare the above expression for the amplitude with its
counterpart in the theory with strictly zero bulk scalar mass. To that end,
one can study the KK spectrum just like we did for the massive bulk theory. As
remarked above, in addition to the KK continuum there is now the true zero
mode in the spectrum, and the amplitude becomes:

$\displaystyle\Delta_{\nu=2}(p^{2})$
$\displaystyle=\frac{1}{p^{2}+i\epsilon}+\int_{0}^{\infty}dmL\,\frac{|\chi^{(m)}_{2}(0)|^{2}}{p^{2}-m^{2}+i\epsilon}$
$\displaystyle=\frac{1}{p^{2}+i\epsilon}+L^{2}\int_{0}^{\infty}dx\,\,\frac{1}{(pL)^{2}-x^{2}+i\epsilon}\,\,\frac{2}{\pi^{2}x(Y_{1}(x)^{2}+J_{1}(x)^{2})}\,,$
(15)

where the first term on the right hand side of the first line comes from the
$m=0$ zero mode exchange, while the second term — from the exchange of the
$m>0$ KK modes. With the help of a few identities involving Bessel/Hankel
functions, the second line of (2) can be massaged into an integral over
(almost) the entire real axis

$\displaystyle\Delta_{\nu=2}(p^{2})=\frac{1}{p^{2}+i\epsilon}+\frac{L^{2}}{2\pi
i}\,\lim_{\delta\to
0^{+}}\left(\int_{-\infty}^{-\delta}dx+\int_{\delta}^{+\infty}dx\right)\,\frac{1}{(pL)^{2}-x^{2}+i\epsilon}\,\frac{H_{2}^{(1)}(x)}{H_{1}^{(1)}(x)}\,.$
(16)

Figure 1: A contour for computing the second term in Eq. (16).

The integrand on the right hand side of this expression has 3 poles, those at
$x=\pm(pL+i\epsilon)$ and at $x=0$. The last of these poles is avoided by the
$\delta$-prescription in the integral (16) (that is, by the fact that we are
integrating up to $-\delta$, and from $+\delta$). The way we traverse the
former two poles is determined by the standard Feynman $i\epsilon$
prescription. Closing the integration contour by the large semicircle at
infinity, as well as the small semicircle between $-\delta$ and $+\delta$, the
integral can be readily computed. The integrand decays sufficiently fast at
infinity so that the large semicircle contributes nothing. The small
semicircle, on the other hand, does give a non-zero contribution, which has to
be substracted to get the desired term in (16). This situation is
schematically depicted in Fig. 3, where we chose to close the contour in the
upper half of the complex plane — closing it in the lower half would lead to
the same result, as can be straightforwardly checked. The contour on the left
picks up the single pole at $x=pL+i\epsilon$. Computing the contribution of
this pole and substracting from it the contribution of the semicircle on the
right of Fig. 3, we finally get:

$\displaystyle\Delta_{\nu=2}(p^{2})$
$\displaystyle=\frac{1}{p^{2}+i\epsilon}+\frac{L^{2}}{2}\frac{H_{2}^{(1)}(pL)}{(pL)H^{(1)}_{1}(pL)}-\frac{1}{p^{2}}$
$\displaystyle=\text{P.V.}\left(\frac{L^{2}}{2}\frac{H_{2}^{(1)}(pL)}{(pL)H^{(1)}_{1}(pL)}\right)-i\pi\delta(p^{2})\,,$
(17)

where the second and the third terms in the first line stem from the left and
the right contours in Fig. 3 respectively222It is straightforward to see that
the semicircle on the right panel of Fig. 3 provides a finite contribution to
the integral. To that end, note that $H_{2}^{(1)}/H_{1}^{(1)}\propto x^{-1}$
in the limit $x\to 0$, and the integral over the small semicircle $s_{1}$ is
proportional to $\int_{s_{1}}dx/x=-i\pi$., while the first term in the second
line emerges from the three terms in the first line. Eq. (17) gives the final
form of the $m_{5}=0$ amplitude, which exactly agrees with the $m_{5}\to 0$
limit (14) of its counterpart in the theory with non-zero $m_{5}$. Since the
considered scalar has no self interactions, and its interactions with the 5D
graviton fluctuations are neglected, the above-established continuity of a
two-point amplitude establishes continuity of the theory in the bulk scalar
mass.

### The static potential

It is instructive to understand continuity of the theory in the limit
$m_{5}\to 0$ on a simple physical observable – the static Newtonian potential
between brane sources. We will start with looking at the large-distance limit
of the amplitude (12), corresponding to $pL\ll 1$. Furthermore, we will also
assume that the bulk mass is small compared to AdS curvature, $\nu-2\ll 1$.
One can then expand $\Delta_{\nu}(p^{2})$ as follows

$\displaystyle\Delta_{\nu}(p^{2})$
$\displaystyle\approx\frac{L^{2}}{(pL)^{2}\left(1-(\nu-2)\right)-2(\nu-2)-\frac{(pL)^{4}}{2}\,\text{ln}(pL/2i)}\,,$
(18)

where the corrections to the denominator are of order
$(pL)^{6},\,(\nu-2)^{2}\,(pL)^{2}$ and $(\nu-2)\,(pL)^{4}$. Keeping only terms
of order $pL$ and $\nu-2$ in (18), this gives rise to the Yukawa potential

$\displaystyle
V_{1}(r)=\int\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\,e^{i\vec{\mathbf{p}}\cdot\vec{r}}\,\frac{L^{2}}{-(\mathbf{p}L)^{2}-2(\nu-2)}\approx-\frac{e^{-m_{0}r}}{4\pi
r}\,.$ (19)

In general, the correction $V_{2}$ to this expression depends on the relative
magnitude of the two expansion parameters. In what follows, we will focus on
deriving two different limits of this correction.333It is important to note
that the term ‘correction’ for the extra term will not be adequate at
distances $r\gg m_{0}$ at which the Yukawa potential is exponentially
suppressed. In such a regine, the extra term $V_{2}\propto r^{-7}$ provides
the leading contribution, as stressed in the discussion to follow.

#### Large distances: $r\gg m^{-1}_{5}$

At large distances, momentum transfer is the smallest scale in the problem, in
particular $pL\ll\nu-2$, and expanding the amplitude (18) in momentum yields:

$\displaystyle V(r\gg
m^{-1}_{5})=V_{1}+\int\frac{d^{3}\mathbf{p}}{(2\pi)^{3}}\,\,e^{i{\mathbf{p}}\cdot\mathbf{r}}\,\,\frac{L^{2}\,(\mathbf{p}L)^{4}\,\text{ln}(|\mathbf{p}|L/2)}{8(\nu-2)^{2}}=V_{1}-\frac{45L^{5}}{4\pi(\nu-2)^{2}r^{7}}\,\,.$
(20)

At distances under consideration, $r\gg m_{5}^{-1}$, the Yukawa potential
$V_{1}$ is contributed mainly by modes with $m\,\lower
3.22916pt\hbox{$\sim$}\hbox to0.0pt{\hss\raise 1.1625pt\hbox{$>$}}\,m_{5}$ and
is exponentially suppressed. The leading contribution to the potential is thus
given by the second term in (20), as discussed in Ref. [5].

#### Small bulk mass: $m_{5}\to 0$

The previous, large distance limit $r\gg m^{-1}_{5}$ is obviously not
consistent with the massless limit of the bulk scalar. This limit, not
considered in [5], is our main focus in the present section. When $m_{5}$ is
sent to zero (that is, when $m_{5}$ is the smallest scale in the problem, but
not strictly zero), both the mass and the width of the 4D resonance tend to
zero as well, the latter vanishing faster than the former. The Yukawa
potential of Eq. (19) therefore turns into the Newtonian one:

$V_{1}(r)|_{m_{5}\to 0}=-\frac{1}{4\pi r}\,.$ (21)

While this visually resembles the potential due to the exchange of a 4D
massless scalar, it is important to stress that in the setup under
consideration there is no localized massless mode in the problem. The bulk
mass $m_{5}$ is far smaller than any other scale, but it is still not strictly
zero – the KK spectrum therefore does _not_ possess a massless state (zero
mode), and the potential (21) is due to a linear combination of KK modes with
$m>0$. The correction to the Newtonian potential $V_{1}$ can be found by
expanding the $\nu\to 2$ limit of the amplitude (18)

$\displaystyle V_{2}(r)$
$\displaystyle=\frac{L^{2}}{2}\int\frac{d^{3}{\textbf{p}}}{(2\pi)^{3}}\,e^{i\mathbf{p}\cdot\mathbf{r}}\,\text{ln}(|\mathbf{p}|L/2)=-\frac{L^{2}}{8\pi
r^{3}}\,.$ (22)

Summing up the two limiting contributions in (21) and (22) yields the
expression, equivalent to what one would get for the static potential in the
case of a _strictly massless_ bulk scalar. As we have remarked above, in the
latter case there does exist a zero mode in the four-dimensional spectrum, and
it leads to the Newtonian force, equivalent to (21). Moreover, in the
$m_{5}=0$ theory, the gapless continuum of KK modes works exactly as it does
in the case of a small but non-zero $m_{5}$, providing a correction,
equivalent to (22). We have therefore established continuity between the two
theories with strictly vanishing $m_{5}$ and however small, but non-zero
$m_{5}$: the physical spectra of these theories are somewhat different, but
all observables—in particular the gravitational potential between brane
sources—are perfectly continuous.

##  3 Massive Scalar with a Zero Mode

Formally, even for non-vanishing bulk mass $m_{5}$, the bulk equation (6) for
the scalar’s 4D KK modes admits a would-be (normalizable) zero mode solution
with the $z$-profile proportional to $(|z|+L)^{-\nu+2}$. As already remarked
in Sec. 2, however, this solution is incompatible with the _boundary
condition_ (7), which effectively removes it from the physical KK spectrum.
This observation suggests that the entire theory may be made compatible with
the existence of the zero mode—even for non-zero bulk mass—by modifying the
_boundary part_ of the original scalar action (a similar mechanism has been
considered for the case of a massive bulk vector in [6]).

In what follows, we will show that modifying the theory by an extra scalar
mass term, localized at $z=0$ does the job of re-introducing the zero mode on
the brane. Indeed, consider the following theory

$\displaystyle S$
$\displaystyle=S_{\phi}+\frac{\nu-2}{2L}\int_{z=0}d^{4}x\,\phi^{2}\,,$ (23)

where $S_{\phi}$ denotes the action of the “minimal” massive theory (2).
Notice also, that the new boundary term has a “tachyonic’ sign, which,
however, does not lead to inconsistency of the theory – indeed, we will show
that the effects of the “wrong” sign boundary mass will be overwhelmed by the
effects of the “correct” sign bulk mass in $S_{\phi}$. Physically, the
tachyonic mass can be thought of as giving an additional attractive
contribution to the effective (‘volcano’) potential that traps the scalar zero
mode on the brane.

With the additional boundary term, the dynamics is governed by the modified
equation of motion

$\displaystyle\Big{(}-\Box+\partial_{z}^{2}-\frac{3\text{sgn}(z)}{|z|+L}\partial_{z}-\frac{(m_{5}L)^{2}}{(|z|+L)^{2}}+\frac{\nu-2}{L}\,\delta(z)\Big{)}\phi(x,z)=-Lj(x)\delta(z)\,,$
(24)

which yields the following boundary condition on the brane

$\displaystyle\Big{(}\partial_{z}+\frac{\nu-2}{L}\Big{)}\phi|_{z=0}=-\frac{L}{2}j(x)\,.$
(25)

The KK mode wavefunctions can be found by solving the system (24) and (25) in
the absence of sources, which establishes that the 4D spectrum of the theory
indeed consists of the zero mode scalar, in addition to the continuum of
massive KK states! The bulk profiles for these modes read

$\displaystyle\chi^{(0)}_{\nu}(z)$
$\displaystyle=\sqrt{\nu-1}\Big{(}\frac{L}{|z|+L}\Big{)}^{\nu-2}\,,$
$\displaystyle\chi^{(m)}_{\nu}(z)$
$\displaystyle=\sqrt{\frac{mL}{2}}\Big{(}\frac{|z|+L}{L}\Big{)}^{2}\Big{[}\,\frac{J_{\nu-1}(mL)Y_{\nu}(m(|z|+L))-Y_{\nu-1}(mL)J_{\nu}(m(|z|+L))}{Y_{\nu-1}(mL)^{2}+J_{\nu-1}(mL)^{2}}\,\Big{]}\,,$
(26)

where, to avoid notational clutter, we have kept the same notation for the KK
wavefunctions as in the previous section, although the corresponding functions
are of course different in the modified theory at hand. Plugging the KK
decomposition
$\phi(x,z)=\phi^{(0)}\chi^{(0)}_{\nu}(z)+\int_{0}^{\infty}dmL\phi^{(m)}(x)\chi^{(m)}_{\nu}(z)$
into the original action (23) and using the standard orthonormality properties
of the KK wavefunctions, we arrive at the 4D effective theory of the following
form

$\displaystyle S=L\int d^{4}x\Bigg{[}$
$\displaystyle\frac{1}{2}(\partial\phi^{(0)})^{2}+\sqrt{\nu-1}\,\phi^{(0)}j$
$\displaystyle+\int_{0}^{\infty}dmL\,\left(\frac{1}{2}(\partial\phi^{(m)})^{2}-\frac{m^{2}}{2}(\phi^{(m)})^{2}+\chi^{(m)}_{\nu}(0)\,\phi^{(m)}j\,.\right)\Bigg{]}\,,$
(27)

One can see, that all modes are well-behaved despite the “wrong” sign of the
brane mass term in (23). Moreover, the canonically normalized zero mode and
the massive KK modes couple to brane sources with strength, set by the
quantities $\sqrt{(\nu-1)/L}$ and $\chi^{(m)}(0)/\sqrt{L}$ respectively.
Exchange of these modes between brane sources gives rise to the following
amplitude

$\displaystyle\Delta(p^{2})$
$\displaystyle=\frac{\nu-1}{p^{2}+i\epsilon}+\int_{0}^{\infty}dmL\,\frac{|\chi^{(m)}(0)|^{2}}{p^{2}-m^{2}+i\epsilon}$
$\displaystyle=\frac{\nu-1}{p^{2}+i\epsilon}+\frac{2L^{2}}{\pi^{2}}\int_{0}^{\infty}dx\,\frac{1}{(pL)^{2}-x^{2}+i\epsilon}\,\,\frac{1}{x\,(Y_{\nu-1}(x)^{2}+J_{\nu-1}(x)^{2})}\,.$
(28)

By performing manipulations on the second term, very similar to the ones
discussed in the previous section, one can integrate over $x$ thereby arriving
at the final, closed-form expression for the amplitude

$\displaystyle\Delta(p^{2})=\text{P.V.}\,\left(\frac{L^{2}}{2}\frac{H_{\nu}^{(1)}(pL)}{(pL)H^{(1)}_{\nu-1}(pL)}\right)-i\pi\,(\nu-1)\,\delta(p^{2})\,.$
(29)

This expression differs by the order of the Hankel functions involved, as well
as by the strength of the massless pole, from its counterpart (17) in the
massless bulk theory, described in the previous section.

##  4 de Sitter brane and Cosmology

The above discussion has exclusively concerned the case of a flat brane, which
requires tuning the brane tension $\lambda$ against the bulk cosmological
constant $\Lambda$. It is known due to Kaloper [7] and Nihei [8], that upon
detuning these two quantities, one can end up with an inflating (de Sitter)
brane, instead of a flat one. The metric of the corresponding spacetime reads
[7, 8]:

$ds^{2}=\left(\cosh\frac{y}{L}-q\sinh\frac{y}{L}\right)^{2}(-dt^{2}+e^{2Ht}d\vec{x}^{2})+dy^{2}\,,\qquad\left(q\equiv\sqrt{1+(HL)^{2}}\right)\,,$
(30)

where the brane is located at $y=0$ in the given coordinates, $L$ is the
curvature radius of the bulk and $H$ is the (constant) Hubble rate on the
$dS_{4}$ spacetime on the brane worldvolume, determined by $\Lambda$ and
$\lambda$ – the precise relation won’t be important for our purposes. It will
prove convenient to perform a change of coordinates

$e^{y/L}=\frac{q+1}{HL}\,\frac{e^{H(z+z_{0})}-1}{e^{H(z+z_{0})}+1}\,,\qquad
z_{0}\equiv H^{-1}\ln\frac{q+HL+1}{q-HL+1}\,,$ (31)

which removes the coordinate singularity in the line element (30), putting is
into the following form

$ds^{2}=\frac{(HL)^{2}}{\sinh^{2}H(z+z_{0})}\,\left(-dt^{2}+e^{2Ht}d\vec{x}^{2}+dz^{2}\right)\,.$
(32)

In the new coordinates, the brane is located at $z=0$, and we will assume that
the theory is invariant under reflections of this coordinate, $z\to-z$, so
(32) should be viewed as describing the brane’s “positive side”, which we will
exclusively work with in the following discussion. On this side of the brane,
the $z$ coordinate ranges from $0$ to $\infty$ as $y$ ranges from zero to a
finite value, determined by the location of the coordinate singularity in
(30), the Rindler horizon. Hence, the $z$ coordinate does not cover the entire
space covered by the $y$ coordinate, but only its patch ranging from the
origin to the Rindler horizon.

While we could continue the spacetime past the horizon to values of
$y>L\,\text{ln}\Big{(}\frac{q+1}{HL}\Big{)}$, anything happening beyond this
region would not effect the observers on the brane. While a signal sent from
the brane to the horizon would take a finite proper time to arrive, it would
take an infinite amount of time according to an observer on the brane.
Furthermore, the boundary conditions for the KK modes are already completely
determined at the horizon by the condition of normalizability so we set aside
the question of any extension of this spacetime beyond the horizon as
irrelevant for our purposes.

Consider the Euler-Lagrange equation of motion for a massive, $z$-reflection-
even bulk scalar $\phi(x,z)=\sigma(x)\,\chi(z)$ in the background (32):

$\displaystyle\frac{d^{2}\chi}{du^{2}}-3\,\frac{\cosh u}{\sinh
u}\,\frac{d\chi}{du}+\left(\frac{m^{2}}{H^{2}}-\frac{m_{5}^{2}L^{2}}{\sinh^{2}u}\right)$
$\displaystyle=0\,,$ (33)

where, with the Kaluza-Klein decomposition in mind, we have defined
$\Box_{4}\sigma=m^{2}\sigma$, as well as $\bar{z}\equiv z+z_{0}$ and $u\equiv
H\bar{z}$. To further simplify notation, we will also define

$A^{2}=\frac{m^{2}}{H^{2}}\,,\qquad B^{2}=m^{2}_{5}L^{2}\,,$ (34)

so that the equation of interest (33) becomes

$\displaystyle\chi^{\prime\prime}-3\,\frac{\cosh(u)}{\sinh(u)}\,\chi^{\prime}+\left(A^{2}-\frac{B^{2}}{\sinh^{2}(u)}\right)\,\chi=0\,,\qquad\chi^{\prime}|=0\,.$
(35)

Here, the last equation—the boundary condition on the brane—can be obtained by
integrating the bulk equation across the brane and keeping in mind that $\chi$
is a $z$-reflection-even field (again, the vertical stroke denotes evaluation
at $u=Hz_{0}$). A solution to this equation with given $m^{2}$ (that is, given
value of the parameter $A^{2}$) describes a localized, normalizable mode if it
satisfies

$\int_{0}^{\infty}dz\,\sqrt{g}\,g^{00}\,\chi^{2}=\text{finite}\,.$ (36)

This imposes a second boundary condition, this time at $u\to\infty$, on
normalizable modes. In what follows, we will be looking for precisely such
localized modes.

A general solution to (35) reads

$\displaystyle\chi(u)$
$\displaystyle=\cosh(u)^{\frac{1}{2}\left(3-\sqrt{9-4A^{2}}\right)}\,\left(\tanh
u\right)^{2-\sqrt{4+B^{2}}}$
$\displaystyle\times\,\Big{[}C_{2}\,\cdot{}_{2}F_{1}\left(a_{1},b_{1},c_{1},t\right)-C_{1}\cdot\left(\tanh
u\right)^{2\sqrt{4+B^{2}}}\cdot{}_{2}F_{1}\left(a_{2},b_{2},c_{2},t\right)\Big{]}\,,$
(37)

where $C_{1}$ and $C_{2}$ are the two integration constants and the
coefficients $a_{1,2},b_{1,2}$ and $c_{1,2}$, together with the $u$-dependent
quantity $t$ have been defined as follows

$\displaystyle a_{1}$
$\displaystyle=\frac{1}{4}\left(1+\sqrt{9-4A^{2}}-2\sqrt{4+B^{2}}\right),b_{1}=\frac{1}{4}\left(3+\sqrt{9-4A^{2}}-2\sqrt{4+B^{2}}\right),c_{1}=1-\sqrt{4+B^{2}},$
$\displaystyle a_{2}$
$\displaystyle=\frac{1}{4}\left(1+\sqrt{9-4A^{2}}+2\sqrt{4+B^{2}}\right),b_{2}=\frac{1}{4}\left(3+\sqrt{9-4A^{2}}+2\sqrt{4+B^{2}}\right),c_{2}=1+\sqrt{4+B^{2}}\,,$

and $t=\tanh^{2}u$. In order to fix the integration constants $C_{1}$ and
$C_{2}$, we will need to study the behavior of the solution (4) for both large
and small values of $u$ – something we will turn to next.

#### Large $u$

For large $u$ (corresponding to $t=\tanh^{2}u\to 1^{-}$), and for $c-a-b<0$
which is the case for the solution (4),444Note that
$c_{1}-a_{1}-b_{1}=c_{2}-a_{2}-b_{2}=-\sqrt{9-4A^{2}}/2$ . the relevant limit
of the hypergeometric function reads

$\lim_{t\to
1^{-}}{}_{2}F_{1}(a,b,c,t)=(1-t)^{c-a-b}\,\,\frac{\Gamma(c)\Gamma(a+b-c)}{\Gamma(a)\Gamma(b)}\,,$
(38)

which gives in our case:

$\lim_{u\to\infty}{}_{2}F_{1}(a_{1,2},b_{1,2},c_{1,2},\tanh^{2}(u))=\left(\cosh
u\right)^{\sqrt{9-4A^{2}}}\,\frac{\Gamma(c_{1,2})\Gamma(\sqrt{9-4A^{2}}/2)}{\Gamma\left(a_{1,2})\Gamma(b_{1,2}\right)}\,.$
(39)

Plugging this expression into the general solution then yields

$\chi(u\gg 1)\simeq\Gamma(\sqrt{9-4A^{2}}/2)\,\left(\cosh
u\right)^{\frac{1}{2}\left(3+\sqrt{9-4A^{2}}\right)}\,\left(C_{2}\,\frac{\Gamma(c_{1})}{\Gamma(a_{1})\Gamma(b_{1})}-C_{1}\,\frac{\Gamma(c_{2})}{\Gamma(a_{2})\Gamma(b_{2})}\right)\,.$
(40)

This correctly reproduces one of the growing modes at large $u$. To see this,
we note that when $u\to\infty$, the two independent solutions have the form
$\exp(k_{1,2}\,x)$, where $k_{1}$ and $k_{2}$ are the two solutions of the
quadratic equation $k^{2}-3kx+A^{2}=0$; explicitly,
$k_{1,2}=(3\pm\sqrt{9-4A^{2}})/2$ and one can see that (40) reproduces the
$k_{1}$-mode, while the $k_{2}$-mode corresponds to the other, sub-leading
solution. The latter solution describes a localized, normalizable mode (with
our proper definition of normalizability, given in eq. (36)), while the former
corresponds to a non-normalizable mode and should thus be removed. To that
end, we need to tune $C_{1}$ and $C_{2}$ as follows:

$C_{2}=C_{1}\,\frac{\Gamma(a_{1})}{\Gamma(a_{2})}\,\frac{\Gamma(b_{1})}{\Gamma(b_{2})}\,\frac{\Gamma(c_{2})}{\Gamma(c_{1})}\,.$
(41)

With this tuning at hand, one can go ahead and study the behavior of the
solution in the opposite limit – the one corresponding to $u\ll 1$.

#### Small $u$

Let us now try to understand whether one can satisfy the boundary condition in
eq. (35)

$\chi^{\prime}|\equiv\chi^{\prime}(\epsilon)=0\,.$ (42)

Note, importantly, that the value of $u$ on the brane is small

$\epsilon\equiv Hz_{0}\sim HL\ll 1,$

and let us first try to understand the solution analytically. For $x\ll 1$ our
equation becomes

$\chi^{\prime\prime}-\frac{3}{u}\,\chi^{\prime}-\frac{B^{2}}{u^{2}}\,\chi=0\,,$
(43)

which is solved by

$\chi(u\ll
1)\simeq\tilde{C}_{1}\,u^{2-\sqrt{4+B^{2}}}+\tilde{C}_{2}\,u^{2+\sqrt{4+B^{2}}}\,,$
(44)

where $\tilde{C}_{1,2}$ will be expressed in terms of $C_{1,2}$ for our
particular solution of interest – we will give the explicit expressions for
these coefficients below. Importantly, neither of the $\tilde{C}_{1,2}$
automatically vanish for our solution, and the (derivatives of the) two terms
can balance each other to satisfy the boundary condition (42) at $u=\epsilon$.
Taylor-expanding the solution (4) and using (41), we have

$\displaystyle C_{1}^{-1}\,\chi^{\prime}(\epsilon)$
$\displaystyle=(2-\sqrt{4+B^{2}})\,\frac{\Gamma(a_{1})\Gamma(b_{1})\Gamma(c_{2})}{\Gamma(a_{2})\Gamma(b_{2})\Gamma(c_{1})}\,\epsilon^{-\sqrt{4+B^{2}}}\,(\epsilon+\mathcal{O}(\epsilon^{3}))$
$\displaystyle-(2+\sqrt{4+B^{2}})\,\epsilon^{\sqrt{4+B^{2}}}\left(\epsilon+\mathcal{O}(\epsilon^{3})\right)\,.$
(45)

(The precise expression for $C_{1}$—not important for the discussion to
come—can be found from orthonormality of the AdS mode functions.) Examining
this expression, one can see that there certainly exist possibilities for it
to vanish, one of which can be understood as follows: the two terms balance
each other for

$\epsilon^{2\sqrt{4+B^{2}}}=\frac{2-\sqrt{4+B^{2}}}{2+\sqrt{4+B^{2}}}\,\,\frac{\Gamma(a_{1})\Gamma(b_{1})\Gamma(c_{2})}{\Gamma(a_{2})\Gamma(b_{2})\Gamma(c_{1})}\,;$
(46)

under our assumptions, $\epsilon$ is small, and so should be the right hand
side of this expression. There seems to be at least the following possibility
for this: $2-\sqrt{4+B^{2}}$ is small if $B$ is small and on top of that
$\Gamma(c_{1})$ is large, which makes the left hand side small for a generic
choice of $A$. Of course, for this to be true, one has to check that there are
no surprises at higher orders in Taylor expansion in $\epsilon$. We can prove
the absence of such surprises numerically by examining how the exact
(unexpanded) $C_{1}^{-1}\,\chi^{\prime}(x)$ behaves for different $A$ and $B$.

Figure 2: The dependence of $C_{1}^{-1}y^{\prime}(x)$ on the parameter $A$ for
$B=\epsilon=10^{-4}$.

### (Numerical) study of the exact expression for $y^{\prime}(x)$

The closed-form expression for $C_{1}^{-1}\,\chi^{\prime}(x)$ is quite
cumbersome and we will not reproduce it here. Instead, let us choose a
representative value for $\epsilon$, such as $\epsilon=10^{-4}$, and
numerically explore the dependence of this expression on $A$ and $B$. This
dependence is shown, in the form of a 3D plot, on Fig. 1 (where irregularities
arise whenever $a_{1}$ or $b_{1}$ cross negative integers). One can see that
the function of interest comes close to zero for small B, but it is not clear
enough whether it actually crosses zero. To see whether
$C_{1}^{-1}\,\chi^{\prime}(x)$ crosses zero, let us choose a particular value
$B=10^{-4}$ and see if a value for $A$ exist such that this function vanishes.

Figure 3: The dependence of $C_{1}^{-1}y^{\prime}(x)$ on the parameters $A$
and $B$.

The dependence of $C_{1}^{-1}\,\chi^{\prime}(x)$ on $A$ for
$\epsilon=B=10^{-4}$ is shown on Fig. 2. It is clear that a value of $A$
exists, such that the desired function smoothly crosses zero – exactly like we
predicted from the Taylor expansion argument, given above. (We emphasise again
that the present argument is not using any approximation – we are working with
exact expressions here.) This establishes the presence of a localized mode for
the case of a de Sitter brane. It is straightforward to see that continuum of
the KK excitations in this case starts above a gap determined by the 4D
curvature scale:

$\displaystyle m_{KK}^{2}\geq{9H^{2}\over 4}\,,$ (47)

and the localized mode resides within the gap. It is difficult to calculate
the mass of the localized mode. However, one can deduce scaling dependence of
the mass of this mode on other parameters by approximately matching asymptotic
solutions. Doing so, one obtains for the mass of the localized mode,
$m_{*}\sim m/(HL)$, in the approximation when $HL<<1$ and $m_{*}<<3H/2$. Thus,
in the above regime of the parameter space the mass of the localized mode is
significantly larger than the mass of the bulk mode, $m$, but is smaller than
the curvature scale $H$.

### Comments on Cosmology

We will outline here how the 5D massive scalar in $AdS_{5}/Z_{2}$ can be used
to describe an accelerated expansion of the universe. The rate of the
acceleration will depend on the scalar mass, and therefore, choosing the value
of the mass one can either make the expansion suitable for inflation, or for
the late time acceleration.

To begin with consider the 5D massive scalar embedded in the geometry with a
tuned brane and bulk, as in Section 2, so that the 4D brane world-volume
spacetime is flat. It remains flat as long as the scalar is in its vacuum
state, $\phi=0$.

Let us now imagine that at some earlier time the initial value of of the
scalar was nonzero, $\phi=\phi_{0}$, with its time derivative being
negligible. At that initial moment the scalar would add to the bulk energy
density a positive quantity

$\displaystyle\Delta E={1\over 2}m_{5}^{2}\,\phi_{0}^{2}\,.$ (48)

For simplicity, we will assume that this energy density is less in its
magnitude than the bulk AdS negative energy density; thus, the quadratic
scalar potential will lead to a reduction of the magnitude of the negative
energy density in the bulk. As a result, the bulk energy density and the brane
tension will no longer be tuned, and the brane world-volume will at that time
moment acquire a positive space-time curvature [7, 8]

$\displaystyle H=\sqrt{{\Delta E\over M_{5}^{3}}}\,,$ (49)

where $M_{5}$ denotes the (specifically normalized) Planck mass of the 5D
theory. Let us furthermore assume that after the initial moment the 5D scalar
field is going to roll down its quadratic potential slowly during some
reasonable classical interval of time. As a result of the slow, roll the
curvature of the 5D spacetime will change slowly too, and so will the
curvature of the 4D (quasi) de Sitter universe. This expansion can potentially
describe either inflation in the early universe, or the dark energy driven
acceleration, depending on the values of the scalar mass and $\phi_{0}$.

While we will not pursue the detailed studies of these cosmologies here, we
point out two peculiarities of the proposed scheme: first the 4D Planck
constant will also be changing in time as the 5D scalar rolls down its slope.
Second, the 4D fluctuations, as shown in Section 4, will consist of one
localized massive mode below the gap and and the KK continuum above a gap. If
applied to inflation, The localized mode will then be responsible for density
perturbations; its mass is parametrically different from the mass of the 5D
field that’s rolling down. Thus, it might be interesting to work out the
details of such a cosmological model.

### Acknowledgments

We thanks Nemanja Kaloper for useful correspondence. GG is supported by the
NSF grant PHY-1915219. Work of DP and DO at NYU was supported by the Simons
Foundation under the "Origins of The Universe" program.

## References

  * [1] L. Randall and R. Sundrum, “A Large mass hierarchy from a small extra dimension,” Phys. Rev. Lett. 83 (1999) 3370–3373, arXiv:hep-ph/9905221.
  * [2] S. J. Avis, C. J. Isham, and D. Storey, “Quantum Field Theory in anti-De Sitter Space-Time,” Phys. Rev. D 18 (1978) 3565.
  * [3] K. Behrndt and M. Cvetic, “Supersymmetric domain wall world from D = 5 simple gauged supergravity,” Phys. Lett. B 475 (2000) 253–260, arXiv:hep-th/9909058.
  * [4] B. Bajc and G. Gabadadze, “Localization of matter and cosmological constant on a brane in anti-de Sitter space,” Phys. Lett. B 474 (2000) 282–291, arXiv:hep-th/9912232.
  * [5] S. L. Dubovsky, V. A. Rubakov, and P. G. Tinyakov, “Brane world: Disappearing massive matter,” Phys. Rev. D 62 (2000) 105011, arXiv:hep-th/0006046.
  * [6] B. Batell and T. Gherghetta, “Localized U(1) gauge fields, millicharged particles, and holography,” Phys. Rev. D 73 (2006) 045016, arXiv:hep-ph/0512356.
  * [7] N. Kaloper, “Bent domain walls as brane worlds,” Phys. Rev. D 60 (1999) 123506, arXiv:hep-th/9905210.
  * [8] T. Nihei, “Inflation in the five-dimensional universe with an orbifold extra dimension,” Phys. Lett. B 465 (1999) 81–85, arXiv:hep-ph/9905487.


# A Non-Gaussian Stochastic Model for the El Niño Southern Oscillation

L. T. Giorgini<EMAIL_ADDRESS>Nordita, Royal Institute of Technology
and Stockholm University, Stockholm 106 91, Sweden  W. Moon<EMAIL_ADDRESS>Department of Mathematics, Stockholm University 106 91 Stockholm, Sweden
Nordita, Royal Institute of Technology and Stockholm University, Stockholm 106
91, Sweden  N. Chen<EMAIL_ADDRESS>Department of Mathematics,
University of Wisconsin-Madison, Madison, Wisconsin 53706, USA  J. S.
Wettlaufer<EMAIL_ADDRESS>Yale University, New Haven, Connecticut
06520, USA  Nordita, Royal Institute of Technology and Stockholm University,
Stockholm 106 91, Sweden

###### Summary

A non-autonomous stochastic dynamical model approach is developed to describe
the seasonal to interannual variability of the El Niño-Southern Oscillation
(ENSO). We determine the model coefficients by systematic statistical
estimations using partial observations involving only sea surface temperature
data. Our approach reproduces the observed seasonal phase locking and its
uncertainty, as well as the highly non-Gaussian statistics of ENSO. Finally,
we recover the intermittent time series of the hidden processes, including the
thermocline depth and the wind bursts.

Air-sea interactions in the tropical Pacific drive the largest interannual
variability process in climate called the El Niño-Southern Oscillation (ENSO),
with an influence that reaches the higher latitudes via atmospheric and
oceanic teleconnections McPhaden2006 ; GhilRMP . Due to its spatiotemporal
impacts, estimating the seasonal to interannual state of ENSO is essential for
predicting a wide range of regional and global climate phenomena Battisti1995
; Havlin2015 . However, the complexity of ENSO, involving as it does
stochastic forcing from atmospheric transients (e.g., Kleeman1997, ) and
nonlinear air-sea interactions (e.g., An2004, ), poses great challenges in
both modeling and prediction.

There are two traditional methods of modeling ENSO (see e.g., Tang2018, , for
a review). One is to use the state-of-art coupled general circulation models
(GCMs), which solve the governing conservation (energy, mass, momentum)
equations for the rotating planet with many sub-grid-scale parameterizations
of unresolved small-scale processes. The other is to build low-order
statistical models. The former treat the main physical aspects of ENSO and
thus provide reasonable representations of the spatiotemporal patterns and
air-sea interactions from intraseasonal to interannual time scales Jan2005 .
The latter, commonly designed to describe and predict particular ENSO indices,
characterize the large-scale features of the system and are much cheaper
computationally than GCMs Petrova2020 . Nonetheless, despite great progress in
developing these methods, there remain limitations to advancing analysis and
prediction. On the one hand, GCMs incur systematic errors originating from
inaccurate state estimations and incomplete sub-grid-scale parameterizations
associated with the ocean thermocline, tropical instability waves, and the
Madden-Julian Oscillation (MJO) Chen2008 ; Barnston2015 . On the other hand,
the statistical models based on multivariate regressions cannot provide
detailed physical information regarding the phase and intensity of ENSO
Tippet2019 . One direction for progress is to utilize a hybrid strategy that
simultaneously maximizes the advantages of both physically-oriented and
statistical models.

The basic theoretical understanding of the tropical air-sea interactions
underlying ENSO Cane1985 ; Jin1993  has lead to a hierarchy of simple models
of the main processes that control the sea surface temperature (SST) in the
tropical ocean (e.g., Wang2017, ). A prominent approach is the recharge-
discharge model Jin1 , and its extensions and generalizations, derived from
the forced shallow water equations using the two-strip approximation. Despite
it being a two-dimensional linear model, it (and its generalizations) can
capture inter-annual variability, seasonal variability with time-periodic
coefficients Stein2010 ; Levine2015 , and a range of stochastic forcing
representing weather and intraseasonal processes, such as westerly wind bursts
(WWBs) Chen2017 .

Here we describe a two-stage stochastic model approach that captures the
large-scale dynamical and statistical characteristics of ENSO. Stochastic
recharge-discharge oscillator models have used observations to recover
statistics on inter-annual time-scales Burgers1999 ; Burgers2005 , and have
incorporated the seasonal cycle of the Bjerknes feedback to treat ENSO phase-
locking and the spring predictability barrier Stein2010 ; Stein2011 ;
Levine2015 . However, while our approach builds upon the recharge-discharge
model, we incorporate a slowly-varying low-order deterministic component that
systematically treats time-varying coefficients and multiplicative noise to
accurately describe the central intraseasonal to interannual features of ENSO.
In particular, our framework reproduces ENSO’s observed seasonal phase
locking, uncertainty and highly non-Gaussian statistics. However, we clearly
depart from multivariate regression statistical models that require
observations of all state variables. Rather, our non-Gaussian model
incorporates a subset of observations from which we can infer unobserved
quantities and quantify their uncertainties using an efficient and exact data
assimilation scheme, and can thereby accurately recover the difficult to
observe thermocline depth and wind bursts. Finally, by systematically
determining time-dependent model coefficients, we precisely simulate the
seasonal to inter-annual ENSO statistics.

We begin with the following coupled model,

$\begin{split}&\dot{x}(t)=a(t)x(t)+\omega(t)h(t)+N(t)\xi_{x}(t)~{}~{}\textrm{and}\\\
&\dot{h}(t)=-\omega(t)x(t)+\lambda h(t)+\sigma\xi_{h}(t),\end{split}$ (1)

in which $x(t)$ is the averaged SST anomaly in the equatorial central-eastern
Pacific (the so-called Niño 3.4 region), and $h(t)$ is the thermocline depth,
which is a surrogate of the heat content, averaged over the western tropical
Pacific. The time-dependent function, $a(t)$, represents the seasonal
evolution of the Bjerknes feedback and the noise amplitude, $N(t)$, captures
the role of the relatively short time-scale forcing of the SST anomaly $x(t)$.
Both $\xi_{x}(t)$ and $\xi_{h}(t)$ are independent Gaussian white noise
processes. Finally, the other time-dependent function, $\omega(t)$, controls
the coupling of tropical atmosphere and the upper ocean, which is one of the
key parameters that determines the quasi-oscillatory behavior of the dominant
ENSO modes, and in traditional low-order models is treated as a constant
Stein2010 . Here, however, treating it as time-dependent facilitates a more
accurate description of the intraseasonal and interannual statistics.

We emphasize the role of the seasonality of background fields in the tropical
Pacific by making all the time dependent coefficients in (1) periodic on the
annual cycle. We determine the functions $a(t)$ and $N(t)$ using the monthly
statistics of the Niño 3.4 index Moon2017 . The coefficient $a(t)$ is related
to the relationship between the monthly variance $\langle x^{2}(t)\rangle$ and
the covariance between two adjacent months $\langle x(t)x(t+\Delta t)\rangle$.
We use the derived $a(t)$ to estimate $N(t)$ from the statistics of the
derived data; $y=x(t+\Delta t)-x(t)-a(t)x(t)\Delta t$. It has been shown
previously that $a(t)$ and $N(t)$ can be accurately estimated independent of
the other coefficients (see supplementary material in Moon2017 ). The coupling
function, $\omega(t)$, and the constants $\lambda$ and $\sigma$, are
determined by an expectation-maximization algorithm Dempster1977 ;
Sundberg1974 ; Sundberg1976 , which is a useful and efficient statistical
estimation method involving incomplete data. We note that only the Niño 3.4
index data is used as the partial observation for model calibration, whereas
the thermocline depth, $h(t)$, is assumed to be unknown since it is not
directly available from satellite observations. Therefore, as detailed in
Appendix B, C, we infer $h(t)$ together with unknown coefficients $\omega(t)$,
$\lambda$, and $\sigma$ chen2020learning .

Figure 1: Comparison of the observed SST (blue; bottom x-axis) with the model
Eq. (1) (red; top x-axis). The SST time series are shown in (a), where the
model result is one of the random realizations from Eq. (1) that shares the
same length as the observational data from 1983 to 2016. The comparisons of
the PDFs and ACFs are shown in (b) and (c). The SST standard deviation as a
function of calendar month is shown in (d). In (e) and (f), the diagram of the
persistence forecast from the observational data and that from the model can
be compared.

We evaluate the quality of the model by comparing the model statistics with
those of the observed Niño 3.4 index in Fig. 5. A qualitative comparison
between a random realization of the model and the observed data show
comparable variability in Fig. 5(a) and quantitative comparisons of the
probability density functions (PDFs) and the autocorrelation functions (ACFs)
of the SST anomaly shown in Figs. 5(b) and (c) demonstrate that the model is
able to accurately simulate the interannual quasi-oscillatory behavior of
ENSO. Moreover, as seen in Fig. 5(d), the model reproduces monthly standard
deviations of the Niño 3.4 index data with high precision. Therefore, two
important seasonal features of ENSO, “phase locking” to the seasonal cycle
Tziperman1998  and the “spring predictability barrier” Levine2015 , are
accurately captured in this two dimensional model. Phase locking is related to
the activation of the positive Bjerknes feedback from summer to early winter,
which leads to the concentration of abnormal ENSO events during November and
December, as is clear from the maximum of the standard deviation near the end
of the year. This phase locking feature is clearly seen from the model
trajectory (Fig. 5a), where the peaks of nearly all of the major ENSO events
occur in the boreal winter, consistent with the observations. The spring
predictability barrier is related to the increase of the noise magnitude,
$N(t)$, during the spring, which increases model uncertainties during the
boreal summer thereby causing the loss of predictability Moon2017 . Indeed,
this is reflected by the structural similarity between the model and the
observations in the persistence forecast diagrams shown in Figs. 5(e) and (f).

Despite its simplicity, the model system of Eqs. (1) recovers the monthly
standard deviation and ACFs of the Niño 3.4 index with higher accuracy than
coupled GCMs Bellenger2014 , linear and nonlinear statistical models
Kondrashov2005 , and simpler recharge-discharge based approaches Stein2010 .
Our systematic statistical estimation of coefficients enable this stochastic
model to regenerate the statistics of the Niño 3.4 index. In particular, as
shown in Fig. 5(c), the precise representation of the ACFs beyond 2 years
suggests that the overall dynamical features of the thermocline depth, $h(t)$,
are well captured in this model. Variations in $h(t)$, and thus of the
equatorial warm water volume (WWV) in the western tropical Pacific, influence
the SST through vertical advection of temperature anomalies by mean upwelling,
a process known as the “thermocline feedback”. Thus, the thermocline depth
acts as a precursor for the inter-annual prediction of SST anomalies
Meinen2000 .

Although the linear model described above successfully captures the seasonal
variability and the interannual large-scale dynamical features of ENSO, it
only allows Gaussian PDFs, whereas the observed statistics are highly non-
Gaussian. Importantly, a model must include atmospheric wind burst forcing,
which plays a critical role in generating the extreme ENSO events responsible
for non-Gaussianity Thual2016 ; Chen2017 . Therefore, we include wind bursts
in the stochastic model of Eqs. (1) as follows,

$\begin{split}&\dot{x}(t)=a(t)x(t)+\omega(t)h(t)+\alpha_{1}\tau(t)+N(t)\xi_{x}(t),\\\
&\dot{h}(t)=-\omega(t)x(t)+\lambda
h(t)+\alpha_{2}\tau(t)+\sigma\xi_{h}(t)~{}~{}\textrm{and}\\\
&\dot{\tau}(t)=d_{\tau}\tau(t)+\rho(x)\xi_{\tau}(t),\end{split}$ (2)

where we adopt a simple linear process with multiplicative noise to generate
the stochastic wind bursts $\tau(t)$ that drive the SST and thermocline
dynamics. The parameters in the wind burst equation are chosen to accurately
reproduce the observed wind stress (Fig. 6 (d)-(f)), while the coefficients of
Eqs. (1) were slightly modified to take into account the additional process in
the model. Despite the non-Gaussianity of the model, the thermocline depth and
the wind bursts remain as Gaussian processes conditioned on the observed SST.
Moreover, by using the partial observations that involve only the Niño SST
data, we have closed analytic formulae for solving the conditional
distribution of the unobserved thermocline depth and the wind bursts (see
chen2018conditional, , and Appendix B),

$p(h(t),\tau(t)|x(s),s\leq t)\sim\mathcal{N}(\mu_{f}(t),R_{f}(t)),$ (3)

where $\mu_{f}(t)$ and $R_{f}(t)$ are the so-called filter mean and filter
covariance in data assimilation.

The filter estimate (3) of the non-Gaussian model (2) can be used to carry out
an online time series reconstruction of the unobserved variables and analyze
the associated uncertainty, which can also serve as the forecast
initialization. Note that the observational data for the equatorial heat
content (the integrated WWV above the 20oC isotherm between 5N-5S, 120E to 80W
Meinen2000 , which is equivalent to the thermocline depth $h(t)$ in (2)) is
only available since 1983 111The WWV values are derived from the temperature
analyses of the Bureau National Operations Centre (BNOC) at the Australian
Bureau of Meteorology Smith1995 , whereas SST data is available from the pre-
industrial era. Therefore, the available SST data can be incorporated into the
model (2) to reconstruct the thermocline data over a much longer period.

We first validate the accuracy of the recovered $h(t)$ and $\tau(t)$ in the
period after 1983. Fig. 6 (a) shows the filter mean time series of the
thermocline depth, $h(t)$, (red) from the non-Gaussian model (2), which
compares favorably to the observational value (blue) with a pattern
correlation that is nearly $0.9$. Moreover, the small error between the
observed and the recovered time series lies approximately within one standard
deviation of the filter estimated uncertainty (pink shaded region). In
particular, as expected from the physical role of the WWBs, the extreme events
are accurately recovered. Similar conclusions are reached for the wind burst
time series (Fig. 6e). In addition to a skillful state estimation, the
reconstructed wind burst time series has the further influence of filtering
out the independent white noise in observations. As shown in Figs. 6(b, f),
the model PDFs capture those of the observations for both the thermocline
depth and the wind bursts, particularly the significant non-Gaussian skewness
and kurtosis. Moreover, similar accuracy is found in ACFs and the power
spectra as shown Figs. 6(c, d, g, h). The success in recovering unobserved
processes is due to the confluence of the model containing the essential
physics and properly constrained stochastic forcing and coefficients.

Figure 2: Comparison of the observational value (blue) and the recovered time
series from the filter estimate (red) using Eq. (3), based on the non-Gaussian
model (2) for the thermocline depth (a) and the wind bursts (e). The solid red
curves are the filter mean while the shaded red areas show one standard
deviation in the filter estimate (the square root of the filter variance). The
comparison of the PDFs, the ACFs and the power spectra between the
observations and simulations from the non-Gaussian model (2) are shown in
(b)–(d), for the thermocline depth and (f)–(h), for the wind bursts.

As shown Fig. 6, the accurate nonlinear filter estimates of Eq. (3) provide an
online estimation method of the thermocline depth $h(t)$ and the wind bursts
$\tau(t)$ relevant for model initialization. On the other hand, reconstructing
historical data can be achieved by solving a slightly different conditional
distribution

$p(h(t),\tau(t)|x(s),0\leq s\leq T)\sim\mathcal{N}(\mu_{s}(t),R_{s}(t)),$ (4)

where $t\in[0,T]$. This is the so-called smoother estimate. It differs from
the filter in the sense that it exploits all the available information up to
the current time instant. One desirable feature of the nonlinear smoother
estimate in Eq. (4) associated with the non-Gaussian system in Eqs. (2) is
that it can be solved via closed analytic formulae (see Appendix B). We use
the smoother estimate to recover the historical realizations of the two hidden
variables conditioned on the entire observational period of the Niño 3.4 index
from 1870 to present. The results are shown in Fig. 3. Since the thermocline
depth and the wind bursts are the most important precursors for predicting the
state of ENSO, these recovered variables can be used to evaluate and improve
prediction models and provide a guideline for improving parameterizations in
coupled GCMs.

Figure 3: The recovered time series of the hidden variables (red curves; left
y-axis) for (a) the thermocline depth and (b) the wind bursts, using the non-
Gaussian model, Eqs. (2), with the partial observations involving only SST
data (dashed black curves; right y-axis) from 1870 to 2016. The method used
here is the nonlinear smoother estimate, Eq. (4). The blue curves show the
observational data from 1983 to 2016.

Understanding the complexity of climate dynamics presents a dilemma of
completeness, between using solely statistical regression models lack the
essential physics, to global climate models, which have more variables and
parameterizations than observables. Our approach presents a compromise through
hybridization. We combine a slowly-varying low-order dynamical system with an
observationally faithful statistical representation of short-time processes.
In consequence we find a realistic description of ENSO as measured by the
success in (i) reproducing the observed seasonal phase locking and its
uncertainty; (ii) capturing the spring predictability barrier; (iii)
simulating the observed highly non-Gaussian statistics, and (iv) accurately
recovering the intermittent time series of the hidden processes via a
nonlinear but analytically solvable data assimilation scheme. Furthermore, we
suggest that extending the approach described here to examine climate
variability on a wide range of time scales GhilRMP ; PinkPRL , wherein proxy
data are available, will provide important insights into mechanisms and
feedbacks in the climate system. Finally, given the generality of the
stochastic dynamical systems in physics, it is hoped that our approach can be
adopted not only for other problems in climate, but far more broadly in the
physical sciences.

###### Acknowledgements.

L.T.G., W.M. and J.S.W. gratefully acknowledge support from the Swedish
Research Council Grant No. 638-2013-9243. The research of N.C. is partially
funded by the Office of VCRGE at UW-Madison.

##  Appendix A Observational Data

All the observational data used in the main paper are plotted in Fig. (4) and
described here.

  * (1)

Monthly sea surface temperature (SST) anomalies derived from the temperature
analyses of the Bureau National Operations Centre (BNOC) at the Australian
Bureau of Meteorology Link1 . These data have been averaged over the Niño 3.4
region (5N-5S, 170W-120W) and start from the pre-industrial era spanning the
entire period 1870-2016. The SST data are the only observations used in the
main paper to recover the model coefficients and to reconstruct the hidden
processes involved.

  * (2)

Daily thermocline depth data, defined as the integrated warm water volume
(WWV) above the 200C isotherm, downloaded from the NCEP/GODAS reanalysis Link2
and averaged within 5N-5S and 80W-120W in the tropical Pacific Ocean. This
dataset was used in the main paper to test the performance of both the two and
three dimensional models in estimating the hidden state of the thermocline
depth correctly, together with their long-term (or “equilibrium”) statistical
and dynamical features, evaluated by comparing the probability distribution
functions (PDFs) and the autocorrelation functions (ACFs) obtained from the
models and from observations. The thermocline depth data are available only in
the post-industrial era, from 1983. Therefore, it was not possible to compare
with observations the recovered hidden process over its total length, but only
its realization from 1983.

  * (3)

Daily westerly wind bursts (WWBs) at 850hPa from the NCEP/NCAR reanalysis
Link3 . By introducing the WWBs into our three dimensional model (Eqs. 2 in
the main paper) we reproduced the non-Gaussian features observed in ENSO
statistics that are not possible in a linear model. The relevant observations
were only available from 1983 and thus comparisons of the mean state to the
statistics of the hidden processes recovered through the model to observations
began at that date.

Figure 4: Observational data used (a) the Niño 3.4 index, (b) the thermocline
depth averaged over the west Pacific and (c) the westerly wind bursts.

##  Appendix B State Estimation and Data Assimilation

###  B.1 Online data assimilation

Our two-dimensional model (Eqs. 1 in the main text) is a linear model with
time-varying coefficients. Therefore, given the observations of the SST, the
classical Kalman-Bucy filter Kalman1961  can be directly applied to obtain the
state estimation of the thermocline depth. In our three-dimensional non-
Gaussian model (Eqs. 2 of the main text), the multiplicative noise prevents
the use of the Kalman-Bucy filter for the state estimation. Nevertheless,
conditional on the observed SST, the thermocline depth and the wind bursts
remain as Gaussian processes. For such a conditional Gaussian system, the
conditional distribution $p(h(t),\tau(t)|x(s),s\leq
t)\sim\mathcal{N}(\boldsymbol{\mu}_{f}(t),\textbf{R}_{f}(t))$ is Gaussian.
Note that despite the non-Gaussianity, the conditional Gaussian distribution
can be solved via closed analytic formulae chen2018conditional . This
facilitates an efficient data assimilation scheme to recover the states of
$h(t)$ and $\tau(t)$ and their uncertainties through the filter mean
$\boldsymbol{\mu}_{f}(t)$ and the square root of the filter covariance
$\textbf{R}_{f}(t)$, which for our three dimensional model are given by

$\displaystyle{\,\rm
d}\boldsymbol{\mu}_{f}=(\textbf{c}_{0}+\textbf{c}_{1}\boldsymbol{\mu}_{f}){\,\rm
d}t+\textbf{R}_{f}\textbf{C}_{1}^{T}B^{-2}({\,\rm
d}x-(C_{0}+\textbf{C}_{1}\boldsymbol{\mu}_{f}){\,\rm d}t),$ (5a)
$\displaystyle{\,\rm
d}\textbf{R}_{f}=\big{(}\textbf{c}_{1}\textbf{R}_{f}+\textbf{R}_{f}\textbf{c}_{1}^{T}+\textbf{b}\textbf{b}^{T}-(\textbf{R}_{f}\textbf{C}_{1}^{T})B^{-2}(\textbf{R}_{f}\textbf{C}_{1}^{T})^{T}\big{)}{\,\rm
d}t,$ (5b)

respectively, where

$\begin{split}\textbf{c}_{0}&=\begin{bmatrix}-\omega x\\\
0\end{bmatrix},\;\;\;\textbf{c}_{1}=\begin{bmatrix}\lambda&\alpha_{2}\\\
0&d_{\tau}\end{bmatrix},\\\
C_{0}&=ax,\;\;\;\textbf{C}_{1}=\begin{bmatrix}\omega&\alpha_{1}\end{bmatrix},\\\
B&=N,\;\;\;\textbf{b}=\begin{bmatrix}\sigma&0\\\
0&\rho\end{bmatrix}.\end{split}$ (6)

The filter mean and covariance of the two dimensional model can be obtained
from these expressions by setting $\alpha_{1},\alpha_{2}=0$.

###  B.2 Offline smoothing

The data assimilation process exploits the observations up to the current time
instant for the purpose of initialization. On the other hand, for
reconstructing the unobserved variables, conditioned on the entire
observational signal (such as obtaining reanalysis data), the so-called
smoothing technique is preferable. For example, the smoothing of the three-
dimensional non-Gaussian system aims to solve $p(h(t),\tau(t)|x(s),0\leq s\leq
T)\sim\mathcal{N}(\boldsymbol{\mu}_{s}(t),\textbf{R}_{s}(t))$, where
$t\in[0,T]$. The smoothing technique can be regarded as a two-step data
assimilation process, with a forward pass of filtering and a backward pass of
smoothing. Although there are no closed formulae for the smoother mean
$\boldsymbol{\mu}_{s}(t)$ and covariance $\textbf{R}_{s}(t)$, closed analytic
formulae are available for both the systems studied in the main text with
appeal to the conditional Gaussian structure Chen2019

$\displaystyle{\,\rm
d}(-\boldsymbol{\mu}_{s})=\big{(}-\textbf{c}_{0}-\textbf{c}_{1}\boldsymbol{\mu}_{s}+(\textbf{b}\textbf{b}^{T})\textbf{R}_{f}^{-1}(\boldsymbol{\mu}_{f}-\boldsymbol{\mu}_{s})\big{)}{\,\rm
d}t,$ (7a) $\displaystyle\begin{split}{\,\rm
d}(-\textbf{R}_{s})=&-\big{(}(\textbf{c}_{1}+(\textbf{b}\textbf{b}^{T})\textbf{R}_{f}^{-1})\textbf{R}_{s}\\\
&+\textbf{R}_{s}(\textbf{c}_{1}^{T}+(\textbf{b}\textbf{b}^{T})\textbf{R}_{f}^{-1})-\textbf{b}\textbf{b}^{T}\big{)}{\,\rm
d}t,\end{split}$ (7b)

where $\boldsymbol{\mu}_{f},\textbf{R}_{f}$ are the filter mean and covariance
respectively, defined in Eq. (5), while the coefficients
$\textbf{c}_{0},\textbf{c}_{1},\textbf{b}$ have been specified in Eq. (6).

In (7), the terms of the left hand side are understood as

$\begin{split}{\,\rm d}(-\boldsymbol{\mu}_{s})&=\lim_{\Delta{t}\to
0}\boldsymbol{\mu}_{s}(t)-\boldsymbol{\mu}_{s}(t+\Delta{t}),\\\ {\,\rm
d}(-\textbf{R}_{s})&=\lim_{\Delta{t}\to
0}\textbf{R}_{s}(t)-\textbf{R}_{s}(t+\Delta{t}).\end{split}$

##  Appendix C Parameter Estimation

###  C.1 Parameter estimation of the time-varying parameters $a(t)$ and $N(t)$

In this section, we will use the monthly-averaged SST data to construct the
periodic drift and diffusion coefficients, $a(t)$ and $N(t)$, that appear in
the $x$-equation of Eqs. (1) in the main paper viz.,

$\begin{split}&\dot{x}(t)=a(t)x(t)+\omega(t)h(t)+N(t)\xi_{x}(t)~{}~{}\textrm{and}\\\
&\dot{h}(t)=-\omega(t)x(t)+\lambda h(t)+\sigma\xi_{h}(t).\end{split}$ (8)

Here we consider a discretized version of this two dimensional model with
$\Delta t$ the discrete time unit. For simplicity, we set the yearly frequency
of the coefficients equal to $2\pi$. We define $N_{y}$ as the total number of
years of the observed SST data, each of which contains $\Delta t^{-1}$ data
points. We multiply both sides of the $x$-equation by $x$, and upon taking the
average we obtain

$\langle x_{i}x_{i+1}\rangle_{i}-\langle x_{i}^{2}\rangle_{i}=\Delta
t(a_{i}\langle x_{i}^{2}\rangle_{i}-\omega_{i}\langle
h_{i}\,x_{i}\rangle_{i})\;\;\forall i\in(0,\Delta t^{-1}],$ (9)

where the average has been performed only over the quantities in the $i$-th
time step inside each period $T$ with $T=1$. Using the $h$-equation, the last
term on the right hand side of Eq. (9) can be written as

$\begin{split}\langle
h_{i}\,x_{i}\rangle_{i}&\simeq\frac{1}{N_{y}}\left(\sum_{j=1}^{N_{y}}\sum_{s=0}^{j\Delta
t^{-1}+i\Delta t}e^{-\lambda(j\Delta t^{-1}+i\Delta
t-s)}\omega_{s}x_{s}x_{i}\Delta t\,\right)\\\
&\simeq\frac{1}{N_{y}}\bar{\omega}\left(\sum_{j=1}^{N_{y}}\sum_{s=j\Delta
t^{-1}+i\Delta t-\Lambda(\lambda)}^{j\Delta t^{-1}+i\Delta
t}x_{s}\,x_{i}\Delta t\right)\simeq\\\ &\simeq\bar{\omega}\langle
x_{i+\frac{1}{4\Delta t}}\,x_{i}\rangle_{i}\simeq 0.\end{split}$ (10)

Here we have used the fact that $x$ is periodic and oscillates around a mean
of zero with an average frequency of $\pi/2$ (the natural frequency of ENSO is
approximately 4 years), which is four times smaller than that of $\omega$.
Therefore, this last coefficient has been approximated in the sum with its
average over one period. Moreover, we have exploited the effect of the
exponential inside the sum, which provides a cutoff in the lower limit. Thus,
the average of this periodic function is a phase shift of a quarter of period
and we can approximate this term as zero.

We can then write

$a_{i}=\frac{\langle x_{i}x_{i+1}\rangle_{i}-\langle
x_{i}^{2}\rangle_{i}}{\Delta t\langle x_{i}^{2}\rangle_{i}}+\Delta a_{i},$
(11)

where $\Delta a_{i}$ is the subleading correction to $a_{i}$ arising from the
term $\langle h_{i}\,x_{i}\rangle_{i}$.

Taking the square of the $x$-equation in Eq. (8) and then averaging, we can
isolate $N_{i}$ viz.,

$\begin{split}N_{i}^{2}\Delta t=&\langle\left[x_{i+1}-x_{i}-\Delta
t(a_{i}x_{i}-\omega_{i}h_{i})\right]^{2}\rangle_{i}\simeq\\\
\simeq&\langle\left[x_{i+1}-x_{i}-\Delta
ta_{i}x_{i}\right]^{2}\rangle_{i}+\Delta t^{2}\omega_{i}^{2}\langle
h_{i}^{2}\rangle_{i}\\\ &\forall i\in(0,\Delta t^{-1}],\end{split}$ (12)

where all the terms proportional to $\langle h_{i}\,x_{i}\rangle_{i}$ have
been neglected. Now we define $y_{i}=x_{i+1}-x_{i}-\Delta ta_{i}x_{i}$ and
write

$\begin{split}\frac{\langle y_{i}^{2}\rangle_{i}-\langle
y_{i+1}y_{i}\rangle_{i}}{\Delta t}&=N_{i}^{2}+\Delta t(\omega_{i}^{2}\langle
h_{i+1}h_{i}\rangle_{i}-\omega_{i}^{2}\langle h_{i}^{2}\rangle_{i})\simeq\\\
&\simeq N_{i}^{2}+O(\Delta t^{2})\;\;\forall i\in(0,\Delta
t^{-1}].\end{split}$ (13)

Finally we can write

$N_{i}=\sqrt{\frac{\langle y_{i}^{2}\rangle_{i}-\langle
y_{i+1}y_{i}\rangle_{i}}{\Delta t}}+\Delta N_{i}.$ (14)

We have used Eqs. (11), (14) to reconstruct the coefficients $a(t)$ and $N(t)$
from the data and we have plotted them in Fig. (5). From Fig. (5)(a) can be
noticed a good agreement between our estimation of $a(t)$ (blue) and the
seasonal Bjerknes instability index as given by Stein2010  (red) as
$a(t)=-1-\sin(2\pi t)$.

Figure 5: Shape of $a(t)$ and $N(t)$ obtained from the data (blue) and
comparison with the estimation of $a(t)$ used in Stein2010  (red).

###  C.2 Parameter estimation of the constant coefficients using an
expectation-maximization algorithm

An expectation-maximization (EM) algorithm is used for parameter estimation of
both the two- and three-dimensional systems. We use this algorithm because we
only have the SST data over the entire observational period and hence the
system is only “partially observed”. Therefore, the parameters and the
unobserved states of the thermocline depth and the wind bursts must be
estimated simultaneously. The advantage of the EM algorithm is that it
iterates and updates the parameters and the hidden states alternatively.
Moreover, convergence can be guaranteed. In particular, in the E-step the
smoother is applied to estimate the state while in the M-step a maximum
likelihood approach is used for parameter estimation. In the E-step, we have
used for convenience a normalized (i.e., with unitary variance) smoother
estimate of the hidden process, which is the thermocline depth for the two
dimensional model. The details of the algorithm can be found in Dempster1977 .

Below, the EM algorithm is applied to estimate
$\omega=\omega_{0}+\omega_{1}\sin(2\pi t)+\omega_{2}\cos(2\pi t)$ and $\sigma$
in the linear model. These coefficients have been estimated keeping the drift
term $\lambda$ of the $h$-process constant in order to guarantee the stability
of the algorithm. This procedure has been iterated many times for different
values of $\lambda$ choosing that which minimizes both the relative and
spectral relative entropy (see Hou2019  and references therein).

Figure 6 illustrates the iteration of the coefficients towards their optima;
$\omega_{0}=1.5,\;\omega_{1}=0.6,\;\omega_{2}=-0.5$ and $\sigma=0.9$, with the
optimal value of $\lambda=-0.8$. Note that because of the simplicity of the
model structure, the convergence here is global.

Figure 6: Coefficients of the model extracted from real data as a function of
the number of iterations. The estimated optimal values of the coefficients
obtained using this method are
$\omega_{0}=1.5,\;\omega_{1}=0.6,\;\omega_{2}=-0.5,\;\lambda=-0.8,\;\sigma=0.9.$

###  C.3 Parameter estimation of the wind bursts coefficients

We introduced a three dimensional model (Eqs. 2 in the main paper) to take
into account the effects of wind bursts as

$\begin{split}&\dot{x}(t)=a(t)x(t)+\omega(t)h(t)+\alpha_{1}\tau(t)+N(t)\xi_{x}(t),\\\
&\dot{h}(t)=-\omega(t)x(t)+\lambda
h(t)+\alpha_{2}\tau(t)+\sigma\xi_{h}(t)~{}~{}\textrm{and}\\\
&\dot{\tau}(t)=d_{\tau}\tau(t)+\rho(x)\xi_{\tau}(t),\end{split}$ (15)

with the following parameters

$\begin{split}&a(t)\to 1.5a(t),\;N(t)\to 0.8N(t),\\\
&\omega_{0}=1.5,\;\omega_{1}=0.6,\;\omega_{2}=-0.5,\;\\\
&\lambda=-1.5,\;\sigma=0.8,\\\
&\alpha_{1}=1,\;\alpha_{2}=-0.6,\;d_{\tau}=-1.5,\\\
&\rho(x)=4.5(\tanh(x)+1)+8.\\\ &\,\end{split}$ (16)

The parameters that were part of the two dimensional model have been slightly
modified to take into account the additional process of wind bursts.

The coefficient $\rho(x)$ in the $\tau$ equation has two additive components.
The first is a multiplicative noise effect arising from the fact that the
warmer SST in the eastern tropical Pacific is usually accompanied by strong
wind bursts (e.g., An2020, ). The second is an additive noise effect, which
will increase the variability of the quiescent phases of the SST. These
combined effects reproduce the observed non-Gaussian statistics of the model
variables.

We let $d_{\tau}=\pi/2$, so that the damping time scale of the wind bursts is
$1/4$ year, but we emphasize that this is not the time scale for the wind
bursts themselves. Rather, this is the time scale for the accumulation effect
of the wind bursts, tropical convection and the MJO to trigger ENSO. This
choice ensures that the decorrelation time-scale of the wind bursts in the
model is the same as that of the observations.

## References

  * (1) M. J. McPhaden, S. E. Zebiak, and M. H. Glantz, Science 314, 1740 (2006).
  * (2) M. Ghil and V. Lucarini, Rev. Mod. Phys. 92, 035002 (2020).
  * (3) D. S. Battisti, and E. S. Sarachik, Rev. Geophys. 33, 1367 (1995).
  * (4) D. Zhou, A. Gozolchiani, Y. Ashkenazy, and S. Havlin, Phys. Rev. Lett. 115, 268501 (2015).
  * (5) R. Kleeman, and A. M. Moore, J. Atmos. Sci. 54, 753 (1997).
  * (6) S. I. An, and F. F. Jin, J. Clim. 17, 2399 (2004).
  * (7) Y. Tang, R. H. Zhang, T. Liu, W. Duan, D. Yang, F. Zheng, …, and M. Mu, Natl. Sci. Rev. 5, 826 (2018).
  * (8) G. Jan van Oldenborgh, M. A. Balmaseda, L. Ferranti, T. N. Stockdale, and D. L. Anderson, J. Clim. 18, 3240 (2005).
  * (9) D. Petrova, J. Ballester, S. J. Koopman, and X. Rodo, J. Clim. 33, 163 (2020).
  * (10) Y. G. Ham , J. H. Kim, and J. J. Luo, Nature 573, 568 (2019).
  * (11) D. Chen, and M. A. Cane, J. Comp. Phys. 227, 3625 (2008).
  * (12) A. G. Barnston, “Evolution of ENSO Prediction over the Past 40 Years”, in Science and Technology Infusion Climate Bulletin, 40${}^{\textrm{th}}$ NOAA Annual Climate Diagnostics and Prediction Workshop, pages 94-101 (2015).
  * (13) M. K. Tippett, M. Ranganathan, M. L’Heureux, A. G. Barnston, and T. Del Sole, Clim. Dyn. 53, 7497 (2019).
  * (14) M. A. Cane, and S. E. Zebiak, Science 228, 1085 (1985).
  * (15) F. F. Jin, and J. D. Neelin, J. Atmos. Sci. 50, 3477 (1993).
  * (16) C. Wang, C. Deser, J. Y. Yu, P. Di Nezio, and A. Clement, in Coral Reefs of the Eastern Tropical Pacific, edited by P. Glynn, D. Manzello, I. Enochs (Springer, Dordrecht, 2017), p. 85.
  * (17) F. F. Jin, J. Atmos. Sci. 54, 811 (1997).
  * (18) K. Stein, N. Schneider, A. Timmermann, and F. F. Jin, J. Clim. 23, 5629 (2010).
  * (19) A. F. Levine, and M. J. McPhaden, Geophys. Res. Lett. 42, 5034 (2015).
  * (20) N. Chen, and A. J. Majda, Proc. Natl. Acad. Sci. U.S.A. 114, 1468 (2017).
  * (21) G. Burgers, Clim. Dyn. 15, 521-531 (1999)
  * (22) G. Burgers, F. F. Jin, and G. J. Van Oldenborgh, Geophys. Res. Lett. 32, 13 (2005)
  * (23) K. Stein, A. Timmermann, and N. Schneider, Phys. Rev. Lett. 107, 128501 (2011)
  * (24) W. Moon, and J. S. Wettlaufer, Sci. Rep. 7, 44228 (2017).
  * (25) A. P. Dempster, N. M. Laird, and D. B. Rubin, J. R. Stat. Soc.: Series B (Methodological) 39, 1 (1977).
  * (26) R. Sundberg, Scand. J. Stat. 1, 49 (1974).
  * (27) R. Sundberg, Commun. Stat. Sim. Comp. 5, 55 (1976).
  * (28) N. Chen, J. Comp. Phys. 418, 109635 (2020).
  * (29) E. Tziperman, M. A. Cane, S. E. Zebiak, Y. Xue, and B. Blumenthal, J. Clim. 11, 2191 (1998).
  * (30) H. Bellenger, E. Guilyardi, J. Leloup, M. Lengaigne, and J. Vialard, Clim. Dyn. 42, 1999 (2014).
  * (31) D. Kondrashov, S. Kravtsov, A. W. Robertson, and M. Ghil, J. Clim. 18, 4425 (2005).
  * (32) C. S. Meinen, and M. J. McPhaden, J. Clim. 13, 3551 (2000).
  * (33) N. Chen, and A. J. Majda, Entropy. 20, 644 (2018).
  * (34) S. Thual, A. J. Majda, N. Chen, and S. N. Stechmann, Proc. Natl. Acad. Sci. U.S.A. 113, 10245 (2016).
  * (35) N. R. Smith, Aust. Met. Mag. 44, 93 (1995).
  * (36) W. Moon, S. Agarwal, and J. S. Wettlaufer, Phys. Rev. Lett. 121, 108701 (2018).
  * (37) BNOC ( Bureau National Operations Centre) Operations Bulletin Number 120 "Operational Implementation of ACCESS-S1 Seasonal Prediction System (Dec 2018)"
  * (38) D. W. Behringer, M. Ji, and A. Leetmaa,  Mon. Weather Rev. 126, 1013 (1998).
  * (39) E. Kalnay, M. Kanamitsu, R. Kistler, W. Collins, D. Deaven, L. Gandin, M. Iredell, S. Saha, G. White, J. Woollen, et al.,  Bull. Am. Met. Soc. 77, 437 (1996).
  * (40) R. E. Kalman, and R. S. Bucy, J. Basic Eng. 83, 95 (1961).
  * (41) N. Chen, and A. J. Majda, J. Comp. Phys. 410, 109381 (2020).
  * (42) A. P. Dempster, N. M. Laird, and D. B. Rubin, J. R. Stat. Soci.: Series B (Methodological) 39, 22, (1977).
  * (43) N. Chen, X. Hou, Q. Li, and Y. Li, Atmosphere 10, 248 (2019).
  * (44) S. An, S. Kim, and A. Timmermann, Sci. Rep. 10, 16282 (2020).