# Spin fluctuations in the 112-type iron-based superconductor
Ca0.82La0.18Fe0.96Ni0.04As2

Tao Xie1,2,3∗, Chang Liu1,2∗, Ryoichi Kajimoto4, Kazuhiko Ikeuchi5, Shiliang
Li1,2,6 and Huiqian Luo1,6† 1 Beijing National Laboratory for Condensed Matter
Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190,
China  2 School of Physical Sciences, University of Chinese Academy of
Sciences, Beijing 100190, China  3 Neutron Scattering Division, Oak Ridge
National Laboratory, Oak Ridge, Tennessee 37831, USA  4 Materials and Life
Science Division, J-PARC Center, Japan Atomic Energy Agency, Tokai, Ibaraki
319-1195, Japan  5 Neutron Science and Technology Center, Comprehensive
Research Organization for Science and Society, Tokai, Ibaraki 319-1106, Japan
6 Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China
$^*$These authors contributed equally to this work<EMAIL_ADDRESS>
###### Abstract

We report time-of-flight inelastic neutron scattering (INS) investigations on
the spin fluctuation spectrum in the 112-type iron-based superconductor (FeSC)
Ca0.82La0.18Fe0.96Ni0.04As2 (CaLa-112). In comparison to the 122-type FeSCs
with a centrosymmetric tetragonal lattice structure (space group $I4/mmm$) at
room temperature and an in-plane stripe-type antiferromagnetic (AF) order at
low temperature, the 112 system has a noncentrosymmetric structure (space
group $P2_{1}$) with additional zigzag arsenic chains between Ca/La layers and
a magnetic ground state with similar wavevector $\mathbf{Q}_{\mathrm{AF}}$ but
different orientations of ordered moments in the parent compounds. Our INS
study clearly reveals that the in-plane dispersions and the bandwidth of spin
excitations in the superconducting CaLa-112 closely resemble to those in 122
systems. While the total fluctuating moments $\langle m^{2}\rangle\approx
4.6\pm 0.2\mu_{B}^{2}$/Fe are larger than 122 system, the dynamic correlation
lengths are similar ($\xi\approx 10$ Å). These results suggest that
superconductivity in iron arsenides may have a common magnetic origin under
similar magnetic exchange couplings with a dual nature from local moments and
itinerant electrons, despite their different magnetic patterns and lattice
symmetries.

###### pacs:

74.70.Xa, 75.30.Gw, 78.70.Nx, 75.40.Gb

Keywords: iron-based superconductors, spin excitations, inelastic neutron
scattering

##  1 Introduction

Unconventional superconductivity usually emerges from the antiferromagentic
(AF) parent compounds such as cuprates, nickelates, iron pnictides and
chalcogenides [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. By introducing charge carrier
doping or applying external pressure to induce superconductivity, the long-
range static AF orders in the parent compounds are gradually suppressed, but
short-range dynamic AF fluctuations persist throughout the superconducting
region and are coupled directly with the occurrence of superconductivity [11,
12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. Thus AF fluctuations are commonly
considered as the pairing glue of the Cooper pairs in unconventional
superconductors [13, 14, 16, 17, 15]. Within this picture, the change of the
magnetic exchange energy ($\Delta E_{ex}$) below and above the superconducting
transition temperature $T_{c}$ should account for the superconducting
condensation energy ($U_{c}$) [23, 24, 25, 26]. Theoretically, it is defined
by: $\Delta
E_{ex}=\Delta\sum_{ij}J_{ij}\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle$,
where $J_{ij}$ is the magnetic exchange coupling and
$\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle$ is the dynamical spin
susceptibility in absolute units [25, 27]. Therefore, it is essential to
experimentally determine both the local magnetic interactions and the absolute
strength of magnetic fluctuations in unconventional superconductors, which can
be measured by inelastic neutron scattering (INS) [13, 14, 15, 16].

While it is a great challenge for INS to map the complete spectrum of spin
fluctuations in high-$T_{c}$ cuprates due to their large energy scale of
bandwidth over than 300 meV [13, 28, 29, 30, 31], the iron-based
superconductors (FeSCs) provide excellent opportunities to extensively compare
the spin fluctuations among different systems due to the relatively lower
bandwidth around 200 meV and available large crystals for different families
with various kinds of doping [14, 15, 16, 17, 32, 33, 34, 35, 36, 37, 38, 39,
26, 40, 41, 42]. Among the FeSCs studied by INS, the 122 systems including
(Ba, Sr, Ca)Fe2As2 and the related doped compounds, have a well established
picture for the doping dependence of the spin fluctuations [32, 33, 34, 35,
36, 37, 38, 39]. Taking BaFe2As2 as an example, while the electron dopings Ni
to Fe sites only suppress the spin excitations below 100 meV [34, 35, 36, 37,
38, 39], the hole dopings K to Ba sites strongly suppress the high energy part
instead [25]. On the other hand, the isovalent dopings P to As sites increase
the effective bandwidth of spin excitations [26]. In spite of the doping
dependent strength and the dramatically different dispersions of spin
fluctuations for various iron-based materials, the next-nearest-neighbor (NNN)
effective coupling in their parent compounds are generally AF and rather
similar (about 20 meV). Additionally, the total fluctuating moments $\langle
m^{2}\rangle\approx 3\mu_{B}^{2}$ per Fe are similar to cuprates [25, 14, 15,
35]. More importantly, a neutron spin resonance mode served as the hallmark of
magnetically driven electron pairing is extensively observed in many FeSCs
[43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57], and the spin
fluctuations in the superconducting compounds are strong enough to make
$\Delta E_{ex}>>U_{c}$ [25, 26, 14, 15].

Figure 1: (a) Crystal structure (above) and magnetic structure (below) of
Ca1-yLayFeAs2 system. The magnetic unit cell marked by
$\mathbf{a}_{\mathrm{M}}$ and $\mathbf{b}_{\mathrm{M}}$ axes is used to define
the reciprocal space in our experiment. (b) Energy dependence of the local
dynamic susceptibility $\chi^{\prime\prime}(E)$ for CaLa-112 in the units of
$\mu_{B}^{2}$/eV/Fe. The solid line is a guide to eyes, the dashed line is the
$\chi^{\prime\prime}(E)$ for BaFe2As2, and the dashed dot line is the DFT+DMFT
result for BaFe2As2 [34, 35, 36]. (c) Energy dependence of the dynamic spin-
spin correlation lengths $\xi$ for CaLa-112. (d) The in-plane dispersions of
spin excitations along high symmetric directions: [0, 0]-[1, 0]-[1, 1]-[0, 0].
The solid dots are obtained from two-gaussian peak fitting of the $H$ or
$K$-cuts at constant-energy window, and the open circles are obtained from the
peak fitting of $E$-cuts at fixed $K$ points. The dashed lines are spin wave
dispersions in the 122-type parent compound CaFe2As2 and BaFe2As2, and the
gradient colors are guides to eyes.

Here, we report INS investigations on the spin fluctuation spectrum in the
112-type FeSC Ca0.82La0.18Fe0.96Ni0.04As2 (CaLa-112) with $T_{c}=22$ K.
Different from those 122/1111/111 systems with a tetragonal lattice structure
(space group I4/mmm or P4/mmm ) at room temperature and an in-plane stripe-
type AF order at low temperature, the 112 system has a unique
noncentrosymmetric structure (space group $P2_{1}$) with additional zigzag
arsenic chains between Ca/La layers and a magnetic ground state with a similar
wavevector $\mathbf{Q}_{\mathrm{AF}}$ but different orientation of magnetic
moments in the parent compounds, which is 45∘ away from the stripe direction
[Fig. 1(a)] [58, 59, 60, 61, 62, 63]. Our previous study on the low energy
spin excitations in the same compound revealed a two-dimensional spin
resonance mode under weak spin-orbit coupling [53]. Using time-of-flight INS
technique, we are able to further map out the high-energy spin excitations in
CaLa-112 and compare with the 122 systems. As shown in Fig. 1 (b)-(d), the
fluctuating strength $\chi^{\prime\prime}(E)$ is initially enhanced at low
energy then decreases to a minimum around 70 meV before increasing again to a
maximum around 180 meV. The in-plane spin-spin correlation lengths at all
energies are similar to that in BaFe2-xNixAs2 ($\xi\approx 10$ Å) [34, 35,
36]. Moreover, the overall dispersions within the FeAs-plane and the energy
bandwidth are also similar to BaFe2As2 and CaFe2As2 [32, 34], suggesting that
the magnetic exchange couplings are similar in these two different systems.

Figure 2: Two dimensional constant-energy slices through the magnetic
excitations of CaLa-112 at energies of $E=11\pm 1$, $21\pm 2$, $40\pm 4$,
$60\pm 4$, $80\pm 10$, $100\pm 10$, $125\pm 10$, $140\pm 10$ and $210\pm 20$
meV. The data in (a), (b), (c,d), (e,f,g,h) and (i) are collected using
$E_{i}=42,73,99,250,449$ meV, respectively. The data sets below 75 meV are
subtracted by a radially symmetric $Q$-dependent background integrated from
the diagonal line of the entire zone $-2<H<2$ and $-2<K<2$, which is mainly
from the phonon scattering of the aluminum sample holders. The data sets above
75 meV are subtracted by a background integrated from $1.8<H<2.2$ and
$-0.2<K<0.2$, which is from the incoherent scattering. The color bars
represent the vanadium normalized absolute spin excitation intensity in the
units of mbarn/sr/meV/f.u. and the dashed red boxes indicate the integration
region for calculating energy-dependent local dynamic susceptibility.

##  2 Experimental details

We used the same sample set of CaLa-112 crystals in our previous INS
experiments [53]. These crystals were grown by self-flux method, and its
composition (Ca0.82La0.18Fe0.96Ni0.04As2) were determined by the inductively
coupled plasma analysis [61]. About 2.3 grams of single crystals ($\sim$1500
pieces) were co-aligned in the scattering plane $[H,0,0]\times[0,0,L]$ defined
by ${\bf Q}=(H,K,L)=(q_{x}a/2\pi,q_{y}b/2\pi,q_{z}c/2\pi)$ in reciprocal
lattice unit (r.l.u.) using the pseudo-orthorhombic magnetic unit cell with
$a_{M}\approx b_{M}\approx 5.54$ Å, $c_{M}=10.27$ Å [Fig. 1(a)]. The mosaics
spread of the co-aligned assembly were about 3.7 degrees in $ab$-plane and 2.8
degrees for out-of-plane case [53]. Time-of-flight INS experiments were
carried out at 4SEASONS chopper spectrometer (BL-01) at J-PARC, Tokai, Japan,
with multiple incident energies $E_{i}$ = 449, 250, 99, 73, 42 meV, $k_{i}$
parallel to the $c$-axis, and chopper frequency $f=$ 250 Hz [64, 65]. All data
were collected at base temperature $T=5$ K and analyzed by the Utsusemi, DAVE
and Horace software packages [66, 67, 68]. To compare easily with spin waves
in BaFe2As2 and CaFe2As2, all data were normalized to the absolute units
(mbarn/sr/meV/f.u.) using a vanadium standard method [32, 34, 35, 36, 69]. The
INS directly measured the differential scattering cross section [25]:
$\frac{d^{2}\sigma}{d\Omega dE}\frac{k_{i}}{k_{f}}=\frac{2(\gamma
r_{e})^{2}}{\pi
g^{2}\mu_{B}^{2}}|F(\mathbf{Q})|^{2}\frac{\chi^{\prime\prime}(Q,E)}{1-\exp(-E/k_{B}T)}.$
Here $F(\mathbf{Q})$ is the magnetic form factor of Fe2+,
$1/[1-\exp(-E/k_{B}T)]$ is the Bose population factor, $k_{B}$ is the
Boltzmann constant, and $\chi^{\prime\prime}(Q,E)$ is the imaginary part of
the dynamic susceptibility. After subtracting the background and correcting
the form factor and Bose factor, the data were multiplied by the constant $\pi
g^{2}\mu_{B}^{2}/(2(\gamma r_{e})^{2})=21.6289$ to covert into the absolute
units of $\mu_{B}^{2}$/eV/f.u.. We calculated the local dynamic susceptibility
using $\chi^{\prime\prime}(E)=\int{\chi^{\prime\prime}({\bf Q},E)d{\bf
Q}}/\int{d{\bf Q}}$ , where $\chi^{\prime\prime}({\bf
Q},E)=(1/3)tr(\chi_{\alpha\beta}^{\prime\prime}({\bf Q},E))$. The total
fluctuating moments can be obtained by further integrating on energy:
$\langle
m^{2}\rangle=(3/\pi)\int{\chi^{\prime\prime}(E)/[1-\exp(-E/k_{B}T)]dE}$
[34, 35, 36]. It should be noticed that there is only one Fe atom in the unit
cell for the single-layer CaLa-112 system but two Fe atoms in double-layer
122-compound, so we compare $\chi^{\prime\prime}(E)$ in the units of
$\mu_{B}^{2}$/eV/Fe for clarity [Fig. 1(b)]. In principle, the absolute
magnitude of magnetic scattering can be also estimated by normalizing it to
phonon intensities [69]. Unfortunately, there were very limited phonon signals
could be found in our measurements. Due to the setup of $k_{i}\parallel c$,
the energy transfer $E$ and momentum transfer along $L$ direction were in
coupled with each other, which made it was difficult to observe the full
phonon signals in the three dimensional reciprocal space. Considering the
large uncertainty in such estimation by phonon intensities, we preferred to
use the vanadium standard method for normalization as we extensively did
before in other iron-based superconductors [34, 35, 36].

Figure 3: Energy dependence of the two-dimensional slices along the $[1,K]$
and $[H,0]$ directions with $E_{i}=450$, 250, 99 and 42 meV for panels (a,c)
and (b,d), respectively. The black and red dashed lines are dispersions of
spin waves in CaFe2As2 and BaFe2As2.

##  3 Results and discussions

In principle, the obtained data sets are four-dimensional (4D) in the three-
axes reciprocal space $[H,K,L]$ plus an energy axis $E$. Since the constant-
energy spin excitations are nearly two-dimensional (2D) in reciprocal space
[53], and the energy transfer $E$ actually couples with $L$ due to
$k_{i}\parallel c$ setup, we present our data details mainly in $[H,K]$ plane
for different energy windows, as shown in Figs. 2 - 4. The 2D constant-energy
slices of spin excitations of CaLa-112 in the $[H,K]$ plane are shown in Fig.
2. The scattering intensity is normalized to absolute units of
mbarn/sr/meV/f.u. using a vanadium standard. The dashed boxes mark the
integration region for calculating energy-dependent local dynamic
susceptibility, which is equivalent to the AF Brillouin zone for the magnetic
unit cell with single Fe2+. As we measured the electron overdoped system,
there is no spin gap due to $J_{c}=J_{s}=0$, which is commonly observed in the
3D spin waves of BaFe2As2 and CaFe2As2 ($E_{g}=10\sim 20$ meV) then suppressed
to zero in the doped compounds [32, 34, 35, 70, 71, 72]. This is consistent
with our previous INS measurements on the same sample at low energies and its
more 2D-like nature in comparison with 122 systems [53, 60, 63, 73, 74]. For
energies below 100 meV [$E=11\pm 1$, $21\pm 2$, $40\pm 4$, $60\pm 4$, $80\pm
10$ meV, Figs. 2 (a)-(e)], the spin excitations form transversely elongated
ellipses centered around the in-plane AF ordering wave vectors
$\mathbf{Q}_{\mathrm{AF}}=(\pm 1,0)$ and (0, $\pm$1), and the peak intensity
decreases with increasing energy. Such features are typically observed in the
low energy spin excitations of those electron doped systems such as
BaFe2-x(Ni, Co)xAs2, which are attributed to the anisotropic contributions
from itinerant electrons due to the mismatched sizes between hole and electron
pockets [12, 38, 39, 35, 36]. When the energy increases to $E=100\pm 10$ meV,
$125\pm 10$ and $140\pm 10$ meV [Figs. 2 (f), (g) and (h)], spin excitations
start to split along the $K$-direction and disperse to the zone boundary
$Q=(\pm 1,\pm 1)$. Finally, the spin excitations become blob-like at the AF
zone boundary ($\pm$1, $\pm$1) when $E=210\pm 20$ meV [Fig. 2 (i)]. As the
high-energy spin excitations are extremely weak, it looks two-fold due to
limited statistics, which suppose to be four-fold like in $[H,K]$ plane in the
pseudo-orthorhombic magnetic unit cell. To improve the statistics, we analysis
the dispersion and intensity based on the four-fold data sets at low energies
and two-fold data sets at high energies.

To firstly determine the rough dispersions of spin excitations in CaLa-112, we
plot the background subtracted scattering for $E_{i}=450$, 250, 99 and 42 meV
projected in the wave vectors $\mathbf{Q}=[1,K]$ and $[H,0]$ and energy space.
As clearly shown in Fig. 3 (a) and (c), the energy top for $K$ dispersion is
about 200 meV, and for $H$ direction is about 150 meV, respectively. Such
behaviors can be attributed to the anisotropic nearest-neighbor (NN) magnetic
exchange coupling, namely $J_{1a}\neq J_{1b}$ in the effective Heisenberg
model [32, 34]. It is reasonable for such system with a monoclinic structure
($\alpha=90^{\circ}$, $\beta=91.4^{\circ}$) [59, 60]. It can be also found
that the low-energy spin excitations ($E<50$ meV) are very strong [Fig. 3 (b)
and (d)]. For comparison, we also plot the spin wave dispersions of BaFe2As2
and CaFe2As2 as dashed lines in Fig. 3, which match overall with the
excitations of CaLa-112.

Figure 4: (a-i) Constant-energy cuts in the spin excitations of CaLa-112 along
the $[1,K]$ and $[H,0]$ directions at different energies, where the wave
vector integration ranges are $0.8<H<1.2$ for the $K$ cuts and $-0.2<K<0.2$
for the $H$ cuts. (j-l) Constant-$Q$ cuts in the spin excitations of CaLa-112
at wave vectors $Q=(1,0.5)$, (1, 0.8), and (1, 1) with thickness $K\pm 0.1$
and $0.8<H<1.2$. The red solid lines are Gaussian fitting results to determine
the dispersions.

To quantitatively illustrate the details of the in-plane dispersion of spin
excitations in CaLa-112, we further show typical constant-energy cuts in Fig.
4 (a)-(i) at different energies along the $\mathbf{Q}=[1,K]$ and $[H,0]$,
respectively. The red solid lines are fitting results with two symmetric
Gaussian functions. The spin excitations become initially incommensurate along
the $[1,K]$ direction for $E=21\pm 2$, $41\pm 4$ and $58\pm 4$ meV [Fig. 4
(a)-(c)], and clearly show a two-peak feature above 80 meV [Fig. 4 (d)-(f)].
The spin excitations seem to be commensurate along the $[H,0]$ direction at
low energies $E=21\pm 2$ and $58\pm 4$ meV [Fig. 4 (g)-(h)], and probably
incommensurate at high energy $E=135\pm 10$ meV with very weak intensity [Fig.
4 (i)]. Constant-$Q$ cuts at $Q=(1,0.5)$, (1, 0.8), and (1, 1) are shown in
Fig. 4 (j)-(l), where the peak positions correspond to the energy of $K$
dispersions approaching the AF zone boundary (open circles in Fig. 1(d)). The
in-plane dispersions both along $[1,K]$ and $[H,0]$ are determined from the
fitting results of above 1D cuts, as shown in Fig. 1 (d) in comparison with
the spin waves of BaFe2As2 and CaFe2As2.

Although it has been argued that the spin waves in BaFe2As2 are more
appropriate to be described by an itinerant model when taking into account
moderate electronic correlation effects [37], a local-moment Heisenberg
Hamiltonian with effective exchange couplings $J_{1a}$, $J_{1b}$ (nearest-
neighbor, NN), $J_{2}$ (NNN) and $J_{c}$ (interlayer) can be also used to fit
the spin waves in twinned samples by considering anisotropic couplings and
dampings [32, 34]. Hence the dispersions are given by:
$E(q)=\sqrt{A_{q}^{2}-B_{q}^{2}}$, where $A_{q}=2S[J_{1b}(\cos(\pi
K)-1)+J_{1a}+2J_{2}+J_{c}+J_{s}]$, $B_{q}=2S[J_{1a}\cos(\pi H)+2J_{2}\cos(\pi
H)\cos(\pi K)+J_{c}\cos(\pi L)]$, $J_{s}$ is the single ion anisotropy
constant, and $q$ is the reduced wave vector away from the AF zone center. For
BaFe2As2, $SJ_{1a}=59.2$ meV, $SJ_{1b}=-9.2$ meV, $SJ_{2}=13.6$ meV,
$SJ_{c}=1.8$ meV, $SJ_{s}=0.084$ meV [34], and other fittings give
$SJ_{c}=0.22$ meV, $SJ_{s}=0.14$ meV [70, 71, 72]. For CaFe2As2, $SJ_{1a}=49$
meV, $SJ_{1b}=-5.7$ meV, $SJ_{2}=18.9$ meV, $SJ_{c}=5.3$ meV, $SJ_{s}=0.063$
meV [32]. As shown by the dashed lines Fig. 1(d), the dispersions of spin
waves in BaFe2As2 and CaFe2As2 are similar under these two sets of exchange
couplings. Except for some slight differences for $K$-dispersion below 50 meV,
which were also observed in BaFe2-xNixAs2 [36], the overall dispersion of
CaLa-112 overlap with the spin waves of BaFe2As2 and CaFe2As2 within error
bar. The in-plane dynamical spin-spin correlation length can be obtained from
the Fourier transforms of the Gaussian fitting results of the constant-energy
cuts in reciprocal space ($\xi=8\ln 2/$FWHM, where FWHM is the full width at
half maximum in the unit of Å-1) [46], which is about 10 Å and similar to
BaFe2-xNixAs2, too [36].

Using the method described in the experimental details, we have calculated the
energy dependence of dynamic local susceptibility $\chi^{\prime\prime}(E)$ by
integrating the intensity of spin excitations in the AF Brillouin zone. Due to
the coupling between $E$ and $L$, we have to carefully choose the energy
windows corresponding to $L\pm 0.5$ where $L=0.5,1,1.5,2,...$, where the first
Brillouin zone in $[H,K]$ plane is marked as the dashed boxes in Fig.2 and the
thickness along $L$ direction should be $L\pm 0.5$ in such single-layered
system. However, due to the $L$-independent signals in CaLa-112 [53], the $L$
integration range $L\pm 0.5$ should be equal to $L\pm 1$ by keeping in mind
that the integration in DAVE and Horace software actually calculates the
average signal over the integrated area. To avoid the interference from the
unclean background at low energies ($E_{i}=$ 42, 73 meV), we perform such
integration within a small area such as $0.7<H<1.3$ and $-0.3<K<0.3$, then
normalize the data to the entire zone marked as the dashed boxes in Fig. 2.
For high energy data ($E_{i}=$ 99, 250, 450 meV), we simply integrate the
intensity in a large area after considering the twinning effect, which is
identical to the dashed boxes in Fig. 2 as illustrated in the case of BaFe2As2
[14]. Similar method was used in the calculation of $\chi^{\prime\prime}(E)$
in 122 system [25, 35, 72]. As shown in Fig. 1(b), the low energy spin
susceptibility below 100 meV in CaLa-112 is indeed stronger than that in
BaFe2As2 system, and the bandwidth, defined as the peak in energy dependence
of $\chi^{\prime\prime}(E)$, is similar for these two compounds ($\sim$ 180
meV). The total fluctuating moments $\langle m^{2}\rangle\approx 4.6\pm
0.2\mu_{B}^{2}$/Fe are slightly stronger than BaFe2As2, which are about
$3.6\mu_{B}^{2}$/Fe. These values are also slightly larger than the results of
cuprates, in which $\langle m^{2}\rangle>1.9\mu_{B}^{2}$/Cu [13, 28, 29, 30,
31]. Using the formula for magnetic moment of a spin $\langle
m^{2}\rangle=(g\mu_{B})^{2}S(S+1)$ (where $g=2$), we can estimate that the
CaLa-112 system has an effective spin $S\approx 0.68$, which likely
corresponds to an $S=1/2$ magnetic ground state. These results are certainly
different from the fully localized case, where $\langle
m^{2}\rangle=24\mu_{B}^{2}$/Fe and $S=2$ under the $3d^{6}$ electronic
configuration [35]. Instead, it may be close to an $S=1$ ground state in the
presence of itinerant electrons [75]. Theoretically, the spin excitations in
FeSCs may be alternatively described by a multi-orbital Hubbard-Hund model
based on the pure itinerant picture as mentioned above [37, 76], in which the
intra- and inter-orbital on-site repulsion $U$ and the Hund’s coupling $J_{H}$
are the effective parameters measuring the electronic correlation strength.
Based on this picture, the density functional theory (DFT) combined with
dynamical mean field theory (DMFT) calculation on the dynamic local
susceptibility can predict the bandwidth of spin excitations in iron
pnictides, as shown by the dashed dot line in Fig. 1(b) [34, 35, 36, 26]. With
increasing the interaction of electrons $U$, the Goldstone mode at low energy
gains additional spectral weight. This could be a possible origin of the
strong peak feature of $\chi^{\prime\prime}(E)$ around 20 meV, since the
electronic correlation in CaLa-112 system may be enhanced by the involvement
of As $4p$ orbitals in hybridization with the Fe $3d$ orbitals [53, 77, 78,
79]. The itinerant nature of the low-energy spin excitations can be further
confirmed by our previous polarized INS measurements, which suggest isotropic
spin excitations in spin space [53]. From the comparison of spin excitations
between CaLa-112 and BaFe2-xNixAs2 systems, we find another fact that the
itinerancy of magnetism only affects the spin excitations below 100 meV in the
electron doped iron pinictides, no matter they are suppressed or enhanced by
dopings [25, 35, 36]. Interestingly, in the electron doped iron chalcogenide
FeSe1-xTex, substitutions on Fe sites by Co, Ni and Cu strongly suppress the
itinerancy, but the localization effects from Cu impurities enhance the low-
energy spin excitations below 100 meV [80]. Therefore, the energy scale
$\sim$100 meV probably is a threshold to separate the dual contributions from
local moments and itinerant electrons in FeSCs.

##  4 Summary

In summary, time-of-flight INS measurements are carried out to map the spin
fluctuation spectrum in the 112-type FeSC Ca0.82La0.18Fe0.96Ni0.04As2. The
obtained results are compared with the parent compounds of 122-type FeSCs,
while the in-plane dispersions, energy bandwidth and spin-spin correlation
lengths are quite similar for both systems, the total fluctuating moments are
stronger than BaFe2As2 probably due to more contributions from itinerant
electrons at low energies. Therefore, the magnetic exchange couplings should
be similar between 112 and 122 systems but the fluctuating effective spin is
not, even though they have different magnetic patterns and lattice symmetries.

## Acknowledgements

This work is supported by the National Key Research and Development Program of
China (Grants No. 2018YFA0704200, No. 2017YFA0303100, and No. 2017YFA0302900),
the National Natural Science Foundation of China (Grants No. 11822411, No.
11961160699, and No. 12061130200), the Strategic Priority Research Program (B)
of the CAS (Grants No. XDB25000000 and No. XDB33000000) and K. C. Wong
Education Foundation (GJTD-2020-01). H. L. is grateful for the support from
the Youth Innovation Promotion Association of CAS (Grant No. Y202001) and
Beijing Natural Science Foundation (Grants No. JQ19002). Work at Oak Ridge
National Laboratory (ORNL) was supported by the U.S. Department of Energy
(DOE), Office of Science, Basic Energy Sciences, Materials Science and
Engineering Division. This work is based on inelastic neutron scattering
experiments performed at the Materials and Life Science Division of J-PARC,
Ibaraki, Japan (Proposal No. 2017B0058).

## References

## References

  * [1] Lee P A, Nagaosa N and Wen X -G 2006 _Rev. Mod. Phys._ 78 17
  * [2] Zhou X, Lee W -S, Imada M, Trivedi N, Phillips P, Kee H -Y, Törmä P and Eremets M 2021 _Nat. Rev. Phys._ 3 462
  * [3] Gu Q, Wen H -H 2022 _The Innovation_ 3 100202
  * [4] Chen X, Dai P, Feng D, Xiang T and Zhang F 2014 _Nat. Sci. Rev._ 1 371
  * [5] Si Q, Yu R and Abrahams E 2016 _Nat. Rev. Mater._ 1 16017
  * [6] White B D, Thompson J D, Maple M B 2015 _Physica C_ 514 246
  * [7] Wu W, Cheng J G, Matsubayashi K, Kong P P, Lin F K, Jin C Q, Wang N L, Uwatoko Y and Luo J L 2014 _Nat. Commun._ 5 5508
  * [8] Cheng J, Matsubayashi K, Wu W, Sun J, Lin F, Luo J and Uwatoko Y 2015 _Phys. Rev. Lett._ 114 117001
  * [9] Yang P, Dong Q, Shan P, Liu Z, Sun J, Dun Z, Uwatoko Y, Chen G, Wang B and Cheng J 2022 _Chin. Phys. Lett._ 39 067401
  * [10] Kamihara Y, Watanabe T, Hirano M, and Hosono H 2008 _J. Am. Chem. Soc._ 130 3296
  * [11] Stewart G R 2011 Rev. Mod. Phys. 83 1589
  * [12] Dai P, Hu J, Dagotto E 2012 Nat. Phys. 8 709
  * [13] Tranquada J M, Xu G and Zaliznyak I A 2014 _J. Magn. Magn. Mater._ 350 148
  * [14] Dai P 2015 Rev. Mod. Phys. 87 855
  * [15] Johnson P D, Xu G, Yin W 2015 _Iron-Based Superconductivity_ (Springer, New York, 2015) Chap. 5
  * [16] Gong D, Luo H 2018 _Acta Phys. Sin._ 67 207407
  * [17] Luo H 2017 _Chin. Sci. Bull._ 62 3955
  * [18] Zhou T, Gao Y, Wang Z D 2020 _Sci. China Phys. Mech. Astron._ 63 287412
  * [19] Luo H, Zhang R, Laver M, Yamani Z, Wang M, Lu X, Wang M, Chen Y, Li S, Chang S, Lynn J W and Dai P 2012 _Phys. Rev. Lett._ 108 247002
  * [20] Lu X, Gretarsson H, Zhang R, Liu X, Luo H, Tian W, Laver M, Yamani Z, Kim Y -J, Nevidomskyy A H, Si Q and Dai P 2013 _Phys. Rev. Lett._ 110 257001
  * [21] Hu D et al 2015 _Phys. Rev. Lett._ 114 157002
  * [22] Gao Q et al 2020 _Chin. Phys. Lett._ 37 087402
  * [23] Stockert O, Arndt J, Faulhaber E, Geibel C, Jeevan H S, Kirchner S, Loewenhaupt M, Schmalzl K, Schmidt W, Si Q and Steglich F 2011 _Nat. Phys._ 7 119
  * [24] Inosov D S, Park J T, Bourges P, Sun D L, Sidis Y, Schneidewind A, Hradil K, Haug D, Lin C T, Keimer B and Hinkov V 2012 _Nat. Phys._ 6 178
  * [25] Wang M, Zhang C, Lu X, Tan G, Luo H, Song Y, Wang M, Zhang X, Goremychkin E A, Perring T G, Maier T A, Yin Z, Haule K, Kotliar G and Dai P 2013 _Nat.Commun._ 4 2874
  * [26] Hu D et al 2016 _Phys. Rev. B_ 94 094504
  * [27] Scalapino D J 2012 _Rev. Mod. Phys._ 84 1383
  * [28] Fujita M, Hiraka H, Matsuda M, Matsuura M, Tranquada J M, Wakimoto S, Xu G and Yamada K 2012 _J. Phys. Soc. Jpn._ 81 011007
  * [29] Headings N S, Hayden S M, Coldea R and Perring T G 2010 _Phys. Rev. Lett._ 105 247001
  * [30] Hayden S M, Aeppli G, Mook H A, Perring T G, Mason T E, Cheong S -W and Fisk Z 1996 _Phys. Rev. Lett._ 76 1344
  * [31] Coldea R, Hayden S M, Aeppli G, Perring T G, Frost C D, Mason T E, Cheong S -W and Fisk Z 2001 _Phys. Rev. Lett._ 86 5377
  * [32] Zhao J, Adroja D T, Yao D, Bewley R, Li S, Wang X F, Wu G, Chen X, Hu J, Dai P 2009 _Nat. Phys._ 5 555
  * [33] Ewings R A, Perring T G, Gillett J, Das S D, Sebastian S E, Taylor A E, Guidi T and Boothroyd A T 2011 _Phys. Rev. B_ 83 214519
  * [34] Harriger L W, Luo H Q, Liu M S, Perring T G, Frost C, Hu J P, Norman M R and Dai P C 2011 _Phys. Rev. B_ 84 054544
  * [35] Liu M, Harriger L W, Luo H, Wang M, Ewings R A, Guidi T, Park H, Haule K, Kotliar G, Hayden S M and Dai P 2012 _Nat. Phys._ 8 376
  * [36] Luo H, Lu X, Zhang R, Wang M, Goremychkin E A, Adroja D T, Danilkin S, Deng G, Yamani Z and Dai P 2013 _Phys. Rev. B._ 88 144516
  * [37] Lu X, Scherer D D, Tam D W, Zhang W, Zhang R, Luo H, Harriger L W, Walker H C, Adroja D T, Andersen B M and Dai P 2018 _Phys. Rev. Lett._ 121 067002
  * [38] Li H F et al 2010 _Phys. Rev. B_ 82 140503(R)
  * [39] Tucker G S, Fernandes R M, Li H -F, Thampy V, Ni N, Abernathy D L, Bud’ko S L, Canfield P C, Vaknin D, Schmalian J and McQueeney R J 2012 _Phys. Rev. B_ 86 024505
  * [40] Lumsden M D, Christianson A D, Goremychkin E A, Nagler S E, Mook H A, Stone M B, Abernathy D L, Guidi T, MacDougall G J, de la Cruz C, Sefat A S, McGuire M A, Sales B C and Mandrus D 2010 _Nat. Phys._ 6 182
  * [41] Lipscombe O J, Chen G, Fang C, Perring T G, Abernathy D L, Christianson A D, Egami T, Wang N, Hu J and Dai P C 2011 _Phys. Rev. Lett._ 106 057004
  * [42] Wang M, Fang C, Yao D, Tan G, Harriger L W, Song Y, Netherton T, Zhang C, Wang M, Stone M B, Tian W, Hu J and Dai P 2011 _Nat. Commun._ 2 580
  * [43] Christianson A D, Goremychkin E A, Osborn R, Rosenkranz S, Lumsden M D, Malliakas C D, Todorov I S, Claus H, Chung D Y, Kanatzidis M G, Bewley R I and Guidi T 2008 _Nature (London)_ 456 930
  * [44] Qiu Y, Bao W, Zhao Y, Broholm C, Stanev V, Tesanovic Z, Gasparovic Y C, Chang S, Hu J, Qian B, Fang M and Mao Z 2009 _Phys. Rev. Lett._ 103 067008
  * [45] Qureshi N, Steffens P, Drees Y, Komarek A C, Lamago D, Sidis Y, Harnagea L, Grafe H -J, Wurmehl S, Büchner B and Braden M 2012 _Phys. Rev. Lett._ 108 117001
  * [46] Zhang C, Yu R, Su Y, Song Y, Wang M, Tan G, Egami T, Fernandez-Baca J A, Faulhaber E, Si Q and Dai P 2013 _Phys. Rev. Lett._ 111 207002
  * [47] Wang M, Luo H, Zhao J, Zhang C, Wang M, Marty K, Chi S, Lynn J W, Schneidewind A, Li S and Dai P 2010 _Phys. Rev. B_ 81 174524
  * [48] Zhang C, Wang M, Luo H, Wang M, Liu M, Zhao J, Abernathy D L, Maier T A, Marty K, Lumsden M D, Chi S, Chang S, Rodriguez-Rivera J A, Lynn J W, Xiang T, Hu J and Dai P 2011 _Sci. Rep._ 1 115
  * [49] Lee C H, Steffens P, Qureshi N, Nakajima M, Kihou K, Iyo A, Eisaki H and Braden M 2011 _Phys. Rev. Lett._ 111 167002
  * [50] Park J T, Friemel G, Li Y, Kim J- H, Tsurkan V, Deisenhofer J, von Nidda H -A K, Loidl A, Ivanov A, Keimer B and Inosov D S 2011 _Phys. Rev. Lett._ 107 177005
  * [51] Wang Q, Shen Y, Pan B, Hao Y, Ma M, Zhou F, Steffens P, Schmalzl K, Forrest T R, Abdel-Hafiez M, Chen X, Chareev D A, Vasiliev A N, Bourges P, Sidis Y, Cao H and Zhao J 2016 _Nat. Mater._ 15 159
  * [52] Ma M, Wang L, Bourges P, Sidis Y, Danilkin S and Li Y 2017 _Phys. Rev. B_ 95 100504(R)
  * [53] Xie T, Gong D, Ghosh H, Ghosh A, Soda M, Masuda T, Itoh S, Bourdarot F, Regnault L -P, Danilkin S, Li S and Luo H 2018 _Phys. Rev. Lett._ 120 137001
  * [54] Xie T, Wei Y, Gong D, Fennell T, Stuhr U, Kajimoto R, Ikeuchi K, Li S, Hu J and Luo H 2018 _Phys. Rev. Lett._ 120 267003
  * [55] Hong W, Song L, Liu B, Li Z, Zeng Z, Li Y, Wu D, Sui Q, Xie T, Danilkin S, Ghosh H, Ghosh A, Hu J, Zhao L, Zhou X, Qiu X, Li S and Luo H 2020 _Phys. Rev. Lett._ 125 117002
  * [56] Xie T, Liu C, Fennell T, Stuhr U, Li S -L and Luo H -Q 2021 _Chin. Phys. B_ 30 127402
  * [57] Liu C, Bourges P, Sidis Y, Xie T, He G, Bourdarot F, Danilkin S, Ghosh H, Ghosh S, Ma X, Li S, Li Y and Luo H 2022 _Phys Rev Lett_ 128 137003
  * [58] Katayama N, Kudo K, Onari S, Mizukami T, Sugawara K, Sugiyama Y, Kitahama Y, Iba K, Fujimura K, Nishimoto N, Nohara M and Sawa H 2013 _J. Phys. Soc. Jpn._ 82 123702
  * [59] Jiang S, Liu C, Cao H, Birol T, Allred J M, Tian W, Liu L, Cho K, Krogstad M J, Ma J, Taddei K M, Tanatar M A, Hoesch M, Prozorov R, Rosenkranz S, Uemura Y J, Kotliar G and Ni N 2016 _Phys. Rev. B_ 93 054522
  * [60] Jiang S, Liu L, Schütt M, Hallas M S, Shen B, Tian W, Emmanouilidou E, Shi A, Luke G M, Uemura Y J, Fernandes R M and Ni N 2016 _Phys. Rev. B_ 93 174513
  * [61] Xie T, Gong D, Zhang W, Gu Y, Huesges Z, Chen D, Liu Y, Hao L, Meng S, Lu Z, Li S, Luo H 2017 _Supercond. Sci. Technol._ 30 095002
  * [62] Yu J, Liu T, Ruan B -B, Zhao K, Yang Q, Zhou M and Ren Z 2021 _Sci. China Phys. Mech. Astron._ 64 267411
  * [63] Jiang C, Song H, Xie T, Liu C, Luo H, Zhao S, Zhang X, Nie X, Meng J, Duan Y, Liu H, Liu S 2020 _Chin. Phys. Lett._ 37 067401
  * [64] Nakamura M, Kajimoto R, Inamura Y, Mizuno F, Fujita M, Yokoo T, Arai M 2009 _J. Phys. Soc. Jpn._ 78 093002\.
  * [65] Kajimoto R _et al._ 2011 _J. Phys. Soc. Jpn._ 80 SB025.
  * [66] Inamura Y, Nakatani T, Suzuki J, Otomo T. 2013 _J. Phys. Soc. Jpn._ 82 SA031
  * [67] Azuah R T, Kneller L R, Qiu Y, Tregenna-Piggott P L W, Brown C M, Copley J R D and Dimeo R M 2009 _J. Res. Natl. Inst. Stan. Technol._ 114 341
  * [68] Ewings R A, Buts A, Le M D, van Duijn J, Bustinduy I, Perring T G 2016 _Nucl. Instrum. Methods Phys. Res. Sect. A_ 834 132
  * [69] Xu G, Xu Z J and Tranquada J M 2013 _Rev. Sci. Instrum._ 84 083906
  * [70] Harriger L W, Schneidewind A, Li S, Zhao J, Li Z, Lu W, Dong X, Zhou F, Zhao Z, Hu J and Dai P 2009 _Phys. Rev.Lett._ 103 087005
  * [71] Park J T, Friemel G, Loew T, Hinkov V, Li Y, Min B H, Sun D L, Ivanov A, Piovano A, Lin C T, Keimer B, Kwon Y S and Inosov D S 2012 _Phys. Rev. B_ 86 024437
  * [72] Luo H, Yamani Z, Chen Y, Lu X, Wang M, Li S, Maier T A, Danilkin S, Adroja D T and Dai P 2012 _Phys. Rev.B_ 86 024508
  * [73] Sóñora D, Carballeira C, Ponte J J, Xie T, Luo H, Li S and Mosqueira J 2017 _Phys. Rev.B_ 96 014516
  * [74] Llovo F, Sóñora D, Mosqueira J, Salem-Sugui Jr S, Sundar S, Alvarenga A D, Xie T, Liu C, Li S and Luo H 2021 _Supercond. Sci. Technol._ 34 115010
  * [75] Wang Q, Shen Y, Pan B, Zhang X, Ikeuchi K, Iida K, Christianson A D, Walker H C, Adroja D T, Abdel-Hafiez M, Chen X, Chareev D A, Vasiliev A N, and Zhao J 2016 _Nat. Commun._ 7 12182
  * [76] Gong D et al 2022 _Front. Phys._ 10 886459
  * [77] Liu X, Liu D, Zhao L, Guo Q, Mu Q, Chen D, Shen B, Yi H, Huang J, He J, Peng Y, Liu Y, He S, Liu G, Dong X, Zhang J, Chen C, Xu Z, Ren Z and Zhou X 2013 _Chin. Phys. Lett._ 30 127402
  * [78] Li M, Liu Z, Zhou W, Yang H, Shen D, Li W, Jiang J, Niu X, Xie B, Sun Y, Fan C, Yao Q, Liu J, Shi Z and Xie X 2015 _Phys. Rev. B_ 91 045112
  * [79] Liu Z, Xing X, Li M, Zhou W, Sun Y, Fan C, Yang H, Liu J, Yao Q, Li W, Shi Z, Shen D and Wang Z 2016 _Appl.Phys. Lett._ 109 042602
  * [80] Wang J, Bao S, Shangguan Y, Cai Z, Gan Y, Li S, Ran K, Ma Z, Winn B L, Christianson A D, Zhong R, Li J, Gu G and Wen J 2022 _Phys. Rev. B_ 105 245129


This version of the paper is submitted to the Quarterly Journal of the Royal
Meteorological Society.
Excitation of mixed Rossby-gravity waves by wave - mean flow interactions on
the sphere

Sándor István Mahó1,*, Sergiy Vasylkevych1, Nedjeljka Žagar1,

1 Meteorological Institute, Center for Earth System Research and
Sustainability, Universität Hamburg, Hamburg, Germany

*<EMAIL_ADDRESS>

## Abstract

The equatorial mixed Rossby-gravity wave (MRGW) is an important contributor to
tropical variability. Its excitation mechanism capable of explaining the
observed MRGW variance peak at synoptic scales remains elusive. This study
investigates wave-mean flow interactions as a generation process for the MRGWs
using the barotropic version of the global Transient Inertia-Gravity And
Rossby wave dynamics model (TIGAR), which employs Hough harmonics as the basis
of spectral expansion, thereby representing MRGWs as prognostic variables.
High accuracy numerical simulations manifest that interactions between waves
emanating from a tropical heat source and zonal mean jets in the subtropics
generate MRGWs with the variance spectra resembling the one observed in the
tropical troposphere. Quantification of spectral tendencies associated with
the MRGW energy growth underscores the significance of wave-mean flow
interactions in comparison to excitation mechanisms driven by external forcing
and wave-wave interactions. The MRGW growth and amplitude depend on the
asymmetry in the zonal mean flow that may explain not only seasonal
variability but also differences between the troposphere and the middle
atmosphere.

##  1 Introduction

Since Yanai and Maruyama (1966) first discovered the mixed Rossby-gravity
waves (MRGWs or the MRG modes) as periodic fluctuations in the meridional wind
in the equatorial stratosphere, their properties and the associated excitation
mechanism have been a subject of numerous studies. Theoretically, the MRGW,
also known as Yanai wave, is a westward-propagating solution of the linearized
shallow water equations on the equatorial $\beta$-plane as derived by Matsuno
(1966).

Initially, the MRGWs were analysed as free modes, however, satellite
irradiance data soon provided evidence that tropospheric MRGWs are coupled to
convection (Chang, 1970; Wallace and Chang, 1972; Hendon and Liebmann, 1991;
Takayabu, 1994; Bergman and Salby, 1994; Wheeler and Kiladis, 1999). The MRGWs
were argued to provide energy for tropical depression type waves over the
equatorial Pacific (Takayabu and Nitta, 1993) that can trigger the formation
of tropical storms (Dickinson and Molinari, 2002). In addition, their vertical
propagation was shown to contribute to the forcing of the easterly phase of
the quasi-biennial oscillation (QBO) (Lindzen and Holton, 1968; Alexander et
al., 2008; Kawatani et al., 2010). Furthermore, MRGWs were shown to contribute
to the MJO initiation Takasuka et al. (2019). In a recent study, tropospheric
MRGWs were associated with the Hadley cell dynamics (Hoskins and Yang, 2021).
Observational studies agree about the dominant zonal wavenumbers of 4-5 of the
MRGWs and periods of about 5 days (e.g. Yanai and Hayashi, 1969; Yanai and
Murakami, 1970; Zangvil and Yanai, 1980).

Early studies of the MRGW excitation mechanism applied the adiabatic
barotropic or baroclinic models and found that MRGWs can be generated by long-
period, large-scale tropospheric heating perturbations asymmetric about the
equator (Holton, 1972). This excitation process is easily demonstrated but can
not explain the observed scale selection. Other excitation mechanisms that
have been studied in depth include the lateral boundary forcing (Mak, 1969),
the wave-CISK (conditional instability of the second kind) mechanism (Hayashi,
1970) and the wave interactions (Raupp and Silva Dias, 2005).

The theory proposed by Mak (1969) and refined by Wilson and Mak (1984) states
that the scale selection of MRGWs is driven by a linear resonance mechanism
that is directly controlled by the properties of the forcing wave exerted at
lateral boundaries. Zhang (1993) later showed that the amplitude of laterally
forced MRGWs changes with the mean flow due to the Doppler-shift effect on the
forcing frequency. Hayashi (1970) argued that instability of the wave-CISK
mechanism, which describes interaction between cumulus convection and
equatorial wave dynamics, is unstable and can produce MRGWs. However, this
mechanism was shown not capable of explaining the scale selection of the MRGWs
(Itoh and Ghil, 1988). As an alternative, Itoh and Ghil (1988) suggested a
non-linear wave-CISK theory in which the MRGWs are amplified by asymmetric
lateral boundary forcing. Finally, a recent study by Raupp and Silva Dias
(2005) suggested wave interactions as a generation mechanism for MRGWs. The
authors numerically solved the non-linear shallow-water equations on the
equatorial $\beta-$plane and argued that MRGWs were excited by a non-linear
mechanism in which an asymmetric heat source produces waves that destabilise
the basic state, which supplies energy for the development of MRGWs.

The present study contributes to the understanding of the MRGWs by showing
that they can be excited by non-linear interactions of wave perturbations and
asymmetric zonal mean flows in the tropics or subtropics. This is demonstrated
by a hierarchy of high-accuracy numerical simulations using idealized and
observed latitudinal profiles of the zonal mean zonal wind.

The new excitation mechanism of MRGWs is presented using the recently
developed barotropic version of the TIGAR model (Transient Inertia-Gravity And
Rossby wave dynamics). TIGAR (Vasylkevych and Žagar, 2021) solves the rotating
shallow water (RSW) equations on the sphere and applies Hough harmonics as
spectral expansion basis. In this framework, the MRGW along with other
eigenmodes of the linearized system (i.e. Kelvin waves, inertia-gravity modes
and Rossby modes on the sphere) is a prognostic variable represented by a
spectral Hough expansion coefficient for each zonal wavenumber. In other
words, the MRG wave is defined by the projection of height and winds onto the
spatial structure of the corresponding mode on the sphere and processes
contributing to its tendencies are quantifiable. MRGW structures have
geopotential height and zonal wind field asymmetric with respect to the
equator and symmetric meridional wind field with equatorial maximum111The
horizontal structure of the zonal wavenumber 1 MRGW for the mean depths from
10 km to about 27 meters is shown at https://modes.cen.uni-
hamburg.de/Hough##part2_2 . Our numerical simulations demonstrate that the
subtropical jet has a prominent role in wave$-$mean flow interactions that
leads to the excitation of the MRG mode.

The remainder of this paper is structured as follows. In Section 2 we describe
the model and its setup for the hierarchy of numerical experiments and
demonstrate the model capabilities. The results are presented in the two
sections focusing on the non-linear excitation in simulations with idealized
(Section 3) and realistic (Section 4) profiles of the background zonal wind.
The discussion and conclusions are presented in Section 5.

##  2 Numerical model formulation and setup

We use the 2D version of the TIGAR model, which was introduced in Vasylkevych
and Žagar (2021). Below we recall its formulation, which is followed by a
description of the linearized version of TIGAR, which is developed for the
purpose of this study. The list of numerical experiments concludes the
section.

###  2.1 The non-linear TIGAR model

The model equations in spherical coordinates
$(\lambda,\varphi)\in[0,2\pi)\times(-\pi/2,\pi/2)$ read as follows:

$\displaystyle\frac{\mathrm{d}u}{\mathrm{d}t}-v\left(2\Omega\sin{\varphi}+\frac{u}{a}\tan{\varphi}\right)+\frac{g}{a\cos{\varphi}}\frac{\partial
h}{\partial\lambda}=-\frac{1}{\tau_{R}}u-F_{u}\,,$ (1a)
$\displaystyle\frac{\mathrm{d}v}{\mathrm{d}t}+u\left(2\Omega\sin{\varphi}+\frac{u}{a}\tan{\varphi}\right)+\frac{g}{a}\frac{\partial
h}{\partial\varphi}=-\frac{1}{\tau_{R}}v-F_{v}\,,$ (1b)
$\displaystyle\frac{\mathrm{d}h}{\mathrm{d}t}+h\nabla\cdot\textbf{V}=Q-\frac{1}{\tau_{N}}\left(h-D\right)-F_{h}\,,$
(1c)

where $\frac{\mathrm{d}}{\mathrm{d}t}$ is the material derivative,
$\textbf{V}=[u,v]$ denotes the horizontal velocity vector composed of the
zonal ($u$) and the meridional component ($v$), $h$ is the total fluid depth,
$D$ is the mean depth, $g$ and $\Omega$ are the gravity constant and the
rotation rate of the Earth respectively, $a$ is the radius of the Earth. The
spectral viscosity denoted as $\textbf{F}=[F_{u},F_{v},F_{h}]$ acts similar to
the one proposed by Gelb and Gleeson (2001) by dampening the smallest resolved
scales in Equation 1 in spectral space. The heat source $Q$ is modelled as a
Gaussian, as explained in Section 2.3 Two other dissipative terms, the so-
called Rayleigh friction in the momentum equations and the Newtonian cooling
in the continuity equation, are defined by their characteristic time scales
$\tau_{R}$ and $\tau_{N}$, respectively.

TIGAR solves Equation 1 by applying the Hough harmonics as basis functions for
the spectral expansion. The Hough harmonics are the eigensolutions of the
linearized rotating shallow-water equations (Longuet-Higgins, 1968), and they
involve three types of solutions: the eastward propagating inertia-gravity
modes (EIG), the westward propagating inertia-gravity modes (WIG), and the
low-frequency Rossby modes. The complete set of eigenmodes involves two
special solutions: the westward propagating equatorial mixed Rossby-gravity
and the eastward propagating equatorial Kelvin mode, which fill the frequency
gap between the slow Rossby and fast inertia-gravity (IG) modes. The solution
of Equation 1 in terms of the Hough harmonics was proposed by Kasahara (1977)
and TIGAR follows the steps therein. First, the non-dimensional variables
$(\tilde{u},\tilde{v},\tilde{h})$ are expanded in terms of the Hough harmonics
as

$[\tilde{u},\tilde{v},\tilde{h}]^{T}(\lambda,\phi,\tilde{t})=\sum_{k,n,l}\
\chi_{n,l}^{k}(\tilde{t})\ \textbf{H}_{n,l}^{k}(\lambda,\phi)\,,$ (2)

where $\chi_{n,l}^{k}$ is the Hough expansion coefficient and
$\textbf{H}_{n,l}^{k}$ is the Hough harmonic with the zonal wavenumber $k$,
meridional mode index $n$ and the type of mode $l$. The latter takes values 1,
2 and 3 for the Rossby, EIG and WIG modes, respectively. Both $k$ and $n$
start with 0. The set of Hough harmonics is uniquely (up to a multiplicative
constant) defined, except for $k=0$ Rossby modes, which are the zonal steady
states of linearized RSW equations. For the latter we use so-called K-modes
derived by Kasahara (1978). The MRG and the Kelvin waves correspond to $n=0$
index for the Rossby and EIG modes, respectively. The relation between the
dimensional and non-dimensional variables is the following:
$\tilde{u}=u/\sqrt{gD}$, $\tilde{v}=v/\sqrt{gD}$ and $\tilde{h}=h/D-1$. The
non-dimensional time $\tilde{t}$ is obtained by multiplying $t$ by $2\Omega$.

The substitution of the expansion (Equation 2) into Equation 1 leads to the
spectral ODE

$\frac{d}{d\tilde{t}}\
\chi_{n,l}^{k}(\tilde{t})+\left(d_{n,l}^{k}+i{\tilde{\omega}}_{n,l}^{k}\right)\
\chi_{n,l}^{k}(\tilde{t})=f_{n,l}^{k}(\tilde{t})\,$ (3)

for the model prognostic variables $\chi_{n,l}^{k}$. In Equation 3
$\tilde{\omega}$ denotes the non-dimensional frequency of the mode $(k,n,l)$.
Term $d_{n,l}^{k}$ stands for the spectral viscosity, which is expressed by
Equations 4-7:

$d_{n,l}^{k}=\left\\{\begin{array}[]{lr}0,&\text{if }N\leq N_{c}\\\ \epsilon
q^{2}N^{2}(N+1)^{2},&\text{if }N>N_{c}\end{array}\right\\}$ (4)

where

$N=k+n,$ (5) $q=\exp\left(\frac{-(N-M)^{2}}{2(N-N_{c})^{2}}\right),$ (6)

and

$\epsilon=\frac{1}{\tau_{SV}M^{3}}\ \ .$ (7)

Note that spectral viscosity only acts beyond a pre-defined cut-off total
wavenumber $N_{c}$. For example, for T170 resolution, $N_{c}$ is set to the
value of 94. The strength of $d_{n,l}^{k}$ is controlled by the dissipation
time $\tau_{SV}$, which is set to 884 hours. In addition, $M$ in Equations 6
and 7 stands for the truncation limit. We employ spectral viscosity to ensure
realistic energy spectra near the truncation limit in the long simulations,
while minimally impacting the dynamics otherwise.

The right-hand side term $f_{n,l}^{k}$ in Equation 3 contains the contribution
of non-linear interactions and forcing as

$f_{n,l}^{k}=\frac{1}{2\pi}\int_{0}^{2\pi}\int_{-\pi/2}^{\pi/2}\left(\textbf{N}+\mathcal{F}\right)\left(\sum_{k^{\prime},\,n^{\prime},\,l^{\prime}}\
\chi_{n^{\prime},\,l^{\prime}}^{k^{\prime}}\
\textbf{H}_{n^{\prime},\,l^{\prime}}^{k^{\prime}}\right)\cdot\left(\textbf{H}_{n,l}^{k}\right)^{*}(\lambda,\varphi)\
d\varphi\ d\lambda\,,$ (8)

where

$\textbf{N}=-\gamma\begin{bmatrix}\tilde{\textbf{V}}\cdot\tilde{\nabla}\tilde{u}-\tilde{u}\tilde{v}\
\tan(\varphi)\\\ \tilde{\textbf{V}}\cdot\tilde{\nabla}\tilde{v}+\tilde{u}^{2}\
\tan(\varphi)\\\
\tilde{\textbf{V}}\cdot\tilde{\nabla}\tilde{h}+\tilde{h}\tilde{\nabla}\cdot\tilde{\textbf{V}}\end{bmatrix}\quad\text{and}\quad\mathcal{F}=\begin{bmatrix}-\frac{1}{\tau_{R}}\tilde{u}\\\
-\frac{1}{\tau_{R}}\tilde{v}\\\ Q-\frac{1}{\tau_{N}}\tilde{h}\end{bmatrix}\,.$
(9)

In the shallow water setting, both the properties of linear waves and the
scaling ratio between linear terms and non-linearity are controlled by a
single non-dimensional quantity $1/\gamma^{2}=(2a\Omega)^{2}/gD$ (Longuet-
Higgins, 1968; Boyd and Zhou, 2008), which is known as Lamb’s parameter. The
forcing and non-linear terms, $\mathcal{F}$ and $\mathbf{N}$ respectively,
appearing in $f_{n,l}^{k}$ are computed pseudo-spectrally, i.e. they are
evaluated in physical space, then transformed to spectral space by applying
FFT and Hough expansion in the zonal and meridional directions, respectively.
To this end Hough functions and their derivatives are precomputed with high
precision from Legendre polynomials; thereby the spectral truncation is the
sole source of numerical errors in the evaluation of spectral tendencies.

The contribution of interactions between different modes to the spectral
tendencies are computed directly in the model by replacing full fields in
$\mathbf{N}$ by the physical space reconstruction of the modes of interest. In
summary, using TIGAR for the present study has a double benefit of combining
high precision computations with the ability to separate the dynamical
contributions of forcing, linear dynamics, and non-linear interactions between
separate modes during the model integration. Note that $\mathbf{N}$ also
includes all interactions of the mean state, which correspond to
$k^{\prime}=0$; thus, in the following, term ”non-linear” will refer to both
wave-wave and wave-mean flow interactions.

###  2.2 Linear and linearized TIGAR

TIGAR can be used in a linear regime that simply omits the evaluation of
$\mathbf{N}$ in Equation 8. We will refer to this setup as ”linear model”. The
linear prognostic equations are

$\displaystyle\frac{\partial u}{\partial
t}-2\Omega\sin{\varphi}v+\frac{g}{a\cos{\varphi}}\frac{\partial
h}{\partial\lambda}=-\frac{1}{\tau_{R}}u-F_{u}\,,$ (10a)
$\displaystyle\frac{\partial v}{\partial
t}+2\Omega\sin{\varphi}u+\frac{g}{a}\frac{\partial
h}{\partial\varphi}=-\frac{1}{\tau_{R}}v-F_{v}\,,$ (10b)
$\displaystyle\frac{\partial h}{\partial
t}+\nabla\cdot\textbf{V}=Q-\frac{1}{\tau_{N}}\left(h-D\right)-F_{h}\,.$ (10c)

”Linearized TIGAR” refers to the setup, which solves the RSW equations
linearized about a prescribed geostrophically balanced mean zonal flow.
Separating the model prognostic variables into the zonal mean part (denoted by
an overbar) and a perturbation (denoted by ’),

$u=\overline{u}(\varphi)+u^{\prime}(\lambda,\varphi,t)\,,\quad
v=v^{\prime}(\lambda,\varphi,t)\,,\quad
h=\overline{h}(\varphi)+h^{\prime}(\lambda,\varphi,t)\,,$

the linearized model equations read

$\displaystyle\frac{\partial u^{\prime}}{\partial t}-v^{\prime}\
2\Omega\sin{\varphi}+\frac{g}{a\cos{\varphi}}\frac{\partial
h^{\prime}}{\partial\lambda}=-\frac{\overline{u}}{a\cos\varphi}\frac{\partial
u^{\prime}}{\partial\lambda}-\frac{v^{\prime}}{a}\frac{\partial\overline{u}}{\partial\varphi}+\frac{\overline{u}v^{\prime}}{a}\tan\varphi-\frac{1}{\tau_{R}}u^{\prime}-F_{u}\,$
(11a) $\displaystyle\frac{\partial v^{\prime}}{\partial t}+u^{\prime}\
2\Omega\sin{\varphi}+\frac{g}{a}\frac{\partial
h^{\prime}}{\partial\varphi}=-\frac{\overline{u}}{a\cos\varphi}\frac{\partial
v^{\prime}}{\partial\lambda}-2\frac{\overline{u}u^{\prime}}{a}\tan\varphi-\frac{1}{\tau_{R}}v^{\prime}-F_{v}\,$
(11b) $\displaystyle\frac{\partial h^{\prime}}{\partial
t}=-\frac{\overline{u}}{a\cos\varphi}\frac{\partial
h^{\prime}}{\partial\lambda}-\frac{v^{\prime}}{a}\frac{\partial\overline{h}}{\partial\varphi}-\frac{\overline{h}}{a\cos\varphi}\frac{\partial
u^{\prime}}{\partial\lambda}+\frac{\overline{h}}{a}\frac{\partial
v^{\prime}}{\partial\varphi}+Q-\frac{1}{\tau_{N}}\left(h^{\prime}-D\right)-F_{h}\,,$
(11c)

which leads to a spectral equation similar to Equation 3, where the new terms
including the means are incorporated into $f_{n,l}^{k}$. Note that unlike
linear TIGAR, the linearized model includes the effects of wave-mean flow
interactions, while neither linear nor linearized model support wave-wave
interactions.

###  2.3 Setup of numerical experiments

Numerical simulations are carried out using geostrophically balanced
background zonal flows, which is prescribed either analytically or derived
from reanalyses. Additional simulations use the background state of rest.
TIGAR generates balanced background flows by projecting the zonal wind field
onto $k=0$ Rossby modes and solving the geostrophic balance equation in
spectral space, thereby obtaining the balanced height field. The formulation
of the background state accounts for the symmetry and strength of the jet. To
compare and quantify the importance of other factors such as linear dynamics,
wave-wave and wave-mean flow interactions, location and symmetry of the wave
source and the mean depth, we perform a number of sensitivity experiments,
which are summarized in Table 1. Note that throughout the article the forcing
is called symmetric whenever $Q$ is centred at the equator, otherwise we refer
to asymmetric forcing.

First, we apply idealized zonal wind profiles with a steady jet (denoted JET).
The JET simulations are further split into SYMJET and ASYJET that include two
jets each located at the subtropics in both hemispheres, which are either
symmetric or asymmetric with respect to the equator. These profiles are shown
in Figure 1 along with an example of the balanced height field. The ASYJET
zonal wind profile has a weaker jet in the Southern Hemisphere (SH) that is
shifted further away from the equator compared to the jet north of the
equator. Simulations with a motionless initial state are denoted NOF.

Figure 1: Latitudinal profiles of the idealized zonal mean background flows
(SYMJET and ASYJET). The balanced height field of SYMJET is also added.

Realistic background zonal wind profiles are derived from ERA5 reanalysis data
(Hersbach et al., 2020). We use seasonal zonally-averaged zonal wind at 500
hPa in years 1993, 1999, 2009, 2012 and 2016 that include both strong El Nino
and La Nina cases. The profiles are shown in Figure 2 for JJA and MAM seasons.

Figure 2: Zonally-averaged zonal wind in ERA5 reanalysis at 500 hPa for (a)
MAM and (b) JJA seasons. Five selected years are 1993, 1999, 2009, 2012 and
2016. The mean profiles are shown by solid black lines.

A visual inspection of Figure 2 suggests that zonally-averaged zonal wind is
more asymmetric in JJA than in MAM. In order to quantify the asymmetry of the
background state, we define a zonal asymmetry measure (ZAM) in the space of
Hough harmonics as the ratio between energies in the asymmetric and symmetric
balanced mean state (i.e., Rossby $k$=0) modes:

$ZAM\equiv\frac{\left|\sum_{n=0}^{M}{\chi_{2n,l=1}^{k=0}}\right|^{2}}{\left|\sum_{n=0}^{M}{\chi_{2n+1,l=1}^{k=0}}\right|^{2}}=\frac{\left|\sum_{n=0}^{M}{\chi_{l=1}^{k=0}(asym.)}\right|^{2}}{\left|\sum_{n=0}^{M}{\chi_{l=1}^{k=0}(sym.)}\right|^{2}}\,.$
(12)

The larger the ZAM, the larger the asymmetry. The mean value of ZAM for JJA
and MAM seasons is 0.321 and 0.053 respectively. The flow asymmetry is greater
for the DJF season compared to SON (not shown), but we find that the
differences are the largest between JJA and MAM and limit our simulations to
these two seasons.

Furthermore, we evaluate the conditions for barotropic instability of the
applied background profiles using the Rayleigh-Kuo criterion (Kuo, 1949) and
find that the jet flanks of ASYJET and SYMJET profiles are unstable. In the
case of the profiles derived from ERA5, the tropics and the midlatitude region
appears barotropically stable, whereas the poles are moderately unstable. The
difference between the stability of the zonally averaged zonal flow in MAM and
JJA is negligible.

Simulations are run with two mean depths: $D=250$ m and $D=400$ m, which
represent two characteristic values of the equivalent depth of the tropical
atmosphere. The former has been commonly used for studying the tropical
response to heating perturbations in adiabatic barotropic models with a lid at
the tropopause (e.g. Kasahara and da Silva Dias, 1986; Raupp and Silva Dias,
2005), whereas the latter is the equivalent depth of vertical mode associated
with deep convection when the stratosphere is included in the computation of
the vertical structure functions (Žagar et al., 2022).

The forcing $Q$ is represented by a Gaussian function with the zonal and
meridional scales $L_{x}$ and $L_{y}$, respectively, central location
($\lambda_{0},\varphi_{0}$) and amplitude $A$:

$Q=A\exp\left(-\frac{(\lambda-\lambda_{0})^{2}}{2L_{x}^{2}}-\frac{(\varphi-\varphi_{0})^{2}}{2L_{y}^{2}}\right)\,.$
(13)

The amplitude of the forcing is $A=1.9$ m/day and 3.05 m/day for $D=250$ m and
$400$ m, respectively, which corresponds to approximately 2.5 K/day heating
rate. This is consistent with typical radiative heating rate of convective
cells in the tropics (McFarlane et al., 2007). The linear damping time scale
is chosen equal for the Rayleigh friction and the Newtonian cooling,
$\tau_{R}=\tau_{N}=5$ days.

Numerical simulations are performed using the triangular truncation T170 that
corresponds to the F128 regular Gaussian grid. The model is integrated in time
using the fourth order exponential Runge-Kutta method with time step of 5
minutes.

Table 1: Summary of the performed numerical experiments. The simulations are defined by the symmetry of the forcing, the shape of the background zonal wind, linearity and and the applied mean depth $D$. NL, LI and LZ refers to non-linear, linear and linearized model setup, respectively. Exp label | Forcing symmetry | Background flow | Linearity | D
---|---|---|---|---
ASYMFOR-NOF | asymmetric | state of rest | NL, LI | 400 m, 250 m
SYMFOR-NOF | symmetric | state of rest | NL | 400 m, 250 m
SYMFOR-SYMJET | symmetric | symmetric jet | NL | 400 m, 250 m
SYMFOR-ASYJET | symmetric | asymmetric jet | NL, LI | 400 m, 250 m
ASYMFOR-ASYJET | asymmetric | asymmetric jet |  | NL
---
LZ, LI
400 m
| SYMFOR-MAM
---
SYMFOR-JJA
symmetric |  | ERA5 mean flows of MAM or
---
JJA
NL | 400 m
| ASYMFOR-MAM
---
ASYMFOR-JJA
asymmetric |  | ERA5 mean flows of MAM or
---
JJA
NL | 400 m

###  2.4 Energy diagnostics in modal space

The energy in $J/kg$ in mode $(k,n,l)$ for $k>0$ is given by

$E_{n,l}^{k}=gD\left|{\chi_{n,l}^{k}}\right|^{2}\,.$ (14)

The sum of energy components $E_{n,l}^{k}$ across all $k$, $n$ and $l$ yields
the sum of the kinetic and available potential energy of the system (Žagar et
al., 2015). Previous studies of the zonal wavenumber energy spectra from
reanalyses typically included the MRGW energy among the Rossby wave energy.
The comparison of energy levels by MRGWs and other modes by Stephan et al.
(2021) showed that the MRGWs contain about 20% of the westward-propagating
non-Rossby wave energy in the troposphere with peak at $k=5-7$. Furthermore,
the subseasonal variability associated with the MRGWs makes about $1/3$ of the
subseasonal variability in non-Rossby modes. While these numbers are
significant, their percentages in the total wave energy (all modes) dominated
by the Rossby waves are small. As an example, the energy spectra from the
SYMFOR-JJA simulation shown in Figure 3 demonstrates a nearly realistic energy
partition. Notably, the Kelvin wave energy has a red spectrum and largest
values at the smallest wavenumber, whereas the MRG mode energy peaks at
synoptic scales in agreement with reanalysis data (Stephan et al., 2021).

Figure 3: Total wave energy spectrum and its partition into contributions from
the IG modes, the Rossby modes, the Kelvin mode and the MRG mode. The spectra
are a result of SYMFOR-JJA simulation in Table 1 run with the JJA background
wind of year 1993. The spectra are averaged over 5 days of simulation.

The second diagnostic is targeted at identifying the role of non-linear
interactions is the evaluation of the contribution of non-linearities and the
forcing to the total energy tendencies. From Equation 14, it follows that for
$k\geq 1$:

$\frac{d}{dt}E_{n,l}^{k}\equiv 2gD\ \Re\left[\left(\frac{d}{d\tilde{t}}\
\chi_{n,l}^{k}\right)\left(\chi_{n,l}^{k}\right)^{*}\right]=2gD\
\Re\Bigg{\\{}\left(f_{n,l}^{k}-q_{n,l}^{k}\right)\left(\chi_{n,l}^{k}\right)^{*}+\left(q_{n,l}^{k}\right)\left(\chi_{n,l}^{k}\right)^{*}\Bigg{\\}}\,.$
(15)

Here, the first and the second term on the right-hand side are the energy
tendencies by non-linear interactions and the forcing respectively, and
$q_{n,l}^{k}$ denotes the projection of the forcing term on the Hough
harmonics. The non-linear tendencies can be further split into interactions
between arbitrary wave components, e.g., IG wave - Rossby wave interactions.
We evaluate wave-wave interactions between Hough modes for $k>0$ and wave-mean
flow interactions between Hough modes $k>0$ and $k=0$ by applying filters in
the non-linear energy tendency in modal space. This diagnostics is
particularly useful for studying wave excitations as it allows to separate
contributions of various factors to the spectral energy tendency of every
mode.

###  2.5 Excitation of MRGWs by external forcing

Placing an asymmetric heat source in the tropics leads to a direct excitation
of the MRGWs, as the forcing associated with heating directly projects onto
the MRG mode, leading to the change in its energy. This experiment is denoted
ASYMFOR-NOF in Table 1 and it is performed both with linear and non-linear
TIGAR using the state of no motion as initial condition. The asymmetric heat
source is centred at 2.5°N.

Figure 4 shows the evolution of the MRG mode energy (panel (a)) and the energy
ratio between MRG mode and the Rossby mode (panel (b)) in the linear and non-
linear models. It reveals that the MRGWs are excited instantaneously in both
models. The energy ratio between the MRG and the Rossby modes is largest
initially and decreases gradually until reaching a quasi-equilibrium state
depending on the damping timescale (at approximately 4-5 days of simulation
time). It is important to notice that there is no qualitative difference
between the linear and non-linear simulations. Thus, wave-wave interactions,
which is a proposed MRGW excitation mechanism by (Raupp and Silva Dias, 2005),
are irrelevant in this setup. While a linear excitation of the MRG wave by an
asymmetric heat source in the tropics has been known for half a century, its
purpose here is to demonstrate how the model and the diagnostics work.

Figure 4: (a) Evolution of the MRG mode energy and (b) the evolution of the
ratio between the MRG mode energy and the Rossby mode energy in the ASYMFOR-
NOF simulations (Table 1). Solid and dashed lines correspond to the
simulations in linear (LI) and non-linear (NL) models, respectively. Purple
and gold curves denote solutions from the experiments with the mean depth
$D=400$ m, whereas blue and silver curves denote the mean depth $D=250$ m.

##  3 Non-linear excitation of MRGWs in the idealized simulations

Now we carry out a number of simulations that provide various details of the
MRGW excitation and evolution in the non-linear framework. First, we provide
evidence that wave-mean flow interactions can excite MRGWs. Then we
concentrate on the scale selection during the excitation process.

Focusing on wave-mean flow interactions, we exclude the direct excitation of
MRGWs by placing the forcing at the equator. The corresponding simulations in
Table 1 are SYMFOR-NOF, SYMFOR-SYMJET and SYMFOR-ASYJET.

Figure 5: Evolution of the MRGW energy in simulations SYMFOR-NOF, SYMFOR-
SYMJET and SYMFOR-ASYJET (see Table 1 for details) during the first 75 hours
of simulation time for mean depths (a) $D=400$ m and (b) $D=250$ m.

Figure 5 shows the time evolution of the MRGW energy in the first 75 hours of
these three simulations. MRGWs grow only in the case of the non-linear
simulation that used an asymmetric background state. This identifies the
interactions between the waves generated by the heat source and the asymmetric
subtropical jet, i.e., wave-mean flow interactions, as the cause of the MRGW
excitation. The fact that at the start of the ASYJET-NL simulation only
symmetric $k>0$ modes are present in the system, which cannot produce
asymmetric modes by triad interactions, rules out the hypothesis that wave-
wave interactions could play a primary role here. Furthermore, the excitation
timescale matches the time when the IG waves excited by the prescribed heat
source reach the jet region. The excitation timescale in the ASYJET-NL
simulations depends on the mean depth $D$ and it is 9 hours and 9.5 hours of
simulation time for $D=400$ m and $D=250$ m, respectively.

Figure 6 shows the mean energy spectra of the MRG mode. The averaging window
is chosen to be 3 days, which is within the typical lifetime of the observed
MRG waves (Zangvil and Yanai, 1980). This time window is also much shorter
than the time required for vortices to develop within the jet region due to
instability. For this reason, instability does not contribute to the MRGW
growth within the analyzed window. Figure 6 demonstrates the dominance of
planetary and synoptic scales in the MRGW response, with a peak in zonal
wavenumber $k=4$ in the ASYJET-NL simulations in case of both mean depths. The
total energy spectrum, with amplitudes 2-3 order of magnitude greater than the
MRGW spectrum, contains no peak at large scales.

Figure 6: MRGW energy spectra of the SYMFOR-NOF, SYMFOR-SYMJET and SYMFOR-
ASYJET (details of simulations in Table 1) for (a) $D=400$ m, and (b) $D=250$
m. Only the non-linear SYMFOR-ASYJET simulation produced significant MRGW
response. Total wave energy is also shown in the SYMFOR-ASYJET simulations.
Figure 7: Spectral energy tendencies of the MRG mode as a function of zonal
wavenumber and simulation time in the non-linear SYMFOR-ASYJET experiment from
Table 1. The MRGW energy tendency is partitioned into the the contributions by
(a) low meridional modes - mean flow interactions and (b) high meridional
modes - mean flow interactions, (c) tendencies due to all non-linear
interactions including wave-mean flow and wave-wave interactions.

To generate MRG waves by wave-mean flow interactions, the westerly jet ought
to be in relative proximity to the tropics. This is verified by carrying out
simulations in which the jets are moved to midlatitudes (i.e., centered at 55°
and 60°, not shown). In such case, the MRGW energy drops by 5 or more orders
of magnitude, and the MRGW energy spectra are white. This is simply explained
by the meridional scale of the MRGW and the need for the wave-mean flow
interactions to take place within the meridional boundaries of MRGWs. This
also means that a crucial role in the MRGW excitation process is played by
interaction of the balanced mean state with waves with low meridional wave
indices that have their strongest signal on the equatorward side of the jet.
This is illustrated in Figure 7 for the non-linear SYMFOR-ASYJET simulation
with $D$ = 400 m by showing separately the MRGW spectral energy tendencies due
to interactions involving low meridional modes $0\leq n\leq 5$ and modes with
high meridional indices $5\leq n\leq 170$. Low meridional modes are comprised
of waves with zonal wind maximum between 30°N and 30°S. Further separation of
tendencies for small $n$ shows that the $n=1$ Rossby mode provides the largest
contribution to the MRGW tendencies although this mode alone can not explain
the energy evolution of MRGWs. While IG modes initiate the interactions by
reaching the jet the fastest, Rossby mode - mean flow interactions become
relevant after about 40 hrs of simulation in case of $D=400$ m, meaning that
their contributions to the spectral energy tendencies become comparable to
that of IG - mean state interactions. Similarly, perturbations from
midlatitudes or locally imposed perturbations on the jet can excite an MRGW
signal depending on the choice of $D$, but the amplitudes of the excited waves
are much smaller than in the case of perturbations emanating from tropical
heating.

The experiments so far showed that the generated MRGWs have the peak signal at
large synoptic scales. Here, we further study the scale selection of the waves
excited by wave-mean flow interactions. The simulations analysed here include
an asymmetric forcing and an asymmetric jet as background flow run with
linear, linearized and non-linear TIGAR (i.e., ASYMFOR-ASYJET simulations from
Table 1). The use of asymmetric forcing ensures that the MRGW is excited in
each case, as evident in the MRGW energy spectra shown in Figure 8.

Figure 8: MRGW energy spectra for the linear (LI), linearized (LZ) and non-
linear (NL) ASYMFOR-ASYJET simulations.

The MRGW energy spectra in Figure 8 are averaged over a 2-day period. Note
that, while the MRGW excitation takes place in all simulations, wave-mean flow
interactions act in the linearized and the non-linear model, whereas wave-wave
interactions are permitted only in the non-linear simulation. No significant
differences can be identified between the MRGW energy spectra of the
linearized and the non-linear simulation in Figure 8 suggesting that wave-wave
interactions do not contribute much to the MRGW signal. Furthermore, all
spectra have a spectral peak at $k=4-6$.

Focusing now on the role of various processes in the growth of MRGW energy, we
analyse the spectra of the absolute energy tendency (Equation 15) in the non-
linear simulation. The tendency is partitioned into the contribution by three
processes: the external forcing, wave-mean flow and wave-wave interactions.
The results are shown in Figure 9. It shows that wave-mean flow interactions
play a dominant role in the MRGW scale selection as they are characterized by
a well-defined spectral peak that is located at $k=3-5$. The forcing, which is
Gaussian, makes the total spectrum more flat at scales $k<15$ and its
amplitude for the selected model setup exceeds that of the wave-mean flow
interactions. The contribution of wave-wave interactions is approximately 3
orders of magnitude smaller than that of the other two mechanisms.

The observed spectral peak in the MRGW energy spectra is somewhat sensitive to
the averaging period and evolves with time (not shown). Nevertheless, we find
that local maxima in the MRGW energy wave spectra are always attributed to the
wave-mean flow interactions. Furthermore, as the mean depth rises, the
relative importance of wave-mean flow interactions increases and can even
exceed the amplitude of energy tendencies due to the forcing at large scales.

Figure 9: Absolute spectral energy tendencies of the MRG mode in the non-
linear ASYMFOR-ASYJET simulation. The contributions by wave-wave, wave-mean
flow interactions and the forcing are denoted WW, WMf and FOR, respectively.
NL+FOR stands for the contribution by the non-linear term and the forcing. The
WW spectrum is multiplied by $10^{2}$. Energy tendencies are averaged over two
days.

##  4 Excitation of MRGWs in simulations with realistic background flow

Now we apply background zonal-mean zonal winds derived from ERA5 for the
boreal spring and summer seasons, as listed in Table 1. Our goal is to verify
to what extent the conclusion based on the idealized flows carry over to a
more realistic setting and to quantify the role of observed asymmetry on the
amplitude of the MRGWs. As discussed in section 2, the JJA season has a higher
asymmetry than the MAM season.

First, we consider simulations in which the excitation by the external forcing
is suppressed, i.e., the forcing is symmetric about the equator. The plots in
Figure 10 show the evolution of the MRG mode energy for the two seasons and
the evolution of the ratio between the MRG and the Rossby mode energy.

Figure 10:  a) MRGW energy evolution and b) evolution of the ratio between the
energy of the MRG mode and the Rossby mode (expressed in %) in the tropical
heat source simulations using realistic background flows. The shades of red
and green correspond to the simulations with the JJA and MAM background flows,
respectively.

Figure 10 reveals large differences in the MRGW energy evolution and its
relative percentage between the simulations with the JJA and MAM seasons, with
substantially larger values for JJA. For example, the maximum of the energy of
the MRG mode on day 4 in the JJA simulations is approximately 5 times larger
compared to the simulations with MAM background flows. The maximum ratio
between the energy of the MRG mode and the Rossby mode is approximately 0.65%
in the JJA simulations, whereas it stays below 0.08% in the MAM simulations.
We note that we did not find any major difference between the stability
profiles of the JJA and MAM flows and we did not observe barotropic
instability development during these simulations, therefore the MRGW
excitation is only a result of wave-mean flow interactions. The 5-day averaged
MRGW energy for every experiment is presented as a function of ZAM, the zonal
asymmetry measure, in Figure 11, which provides a clear evidence that the
larger the zonal asymmetry, the more energetic the MRG waves are. The JJA and
MAM simulations appear as two clusters in Figure 11.

Figure 11: MRGW energy integrated over 5 days of simulation as a function of
the asymmetry measure (ZAM) defined in Section 2.3 The shades of red and green
represent JJA and MAM simulations respectively. Figure 12: Energy spectra of
the MRGWs in simulations with realistic background flows run with a) symmetric
and b) asymmetric forcing. The averaging period is 5 days. Thin and thick
curves denote the individual simulations and the means over JJA and MAM
simulations, respectively.

The mean energy spectra of the MRGWs in the simulations with realistic
background flows is shown in Figure 12 for the case of both symmetric and
asymmetric forcing. The averaging period is 5 days. We find that even with the
excitation by the forcing, simulations with more asymmetric background wind
(JJA) produce more energetic MRG waves at large and at synoptic scales (i.e.
up to $k$ = 9 and 14 in case of the asymmetric and the symmetric forcing
respectively). At the scales $k=3-7$, the JJA simulations yield 10 and 3 times
larger MRGW energy in case of symmetric and asymmetric forcing, respectively.
The differences between the JJA and MAM spectra vanish as the zonal wavenumber
increases towards $k=50$.

The forcing symmetry plays a role for the MRGW scale selection. In case of a
symmetric forcing, the dominant MRGW wavenumber is in range of $k=5-7$ for
both MAM and JJA background flows. When an asymmetric forcing is applied, the
mean MAM energy spectrum flattens without showing any prominent peaks, whereas
the mean JJA spectrum has a well-defined peak in range of $k$=3-7. This result
implies that the spectral tendencies of the MRGW energy due to wave-mean flow
interactions, which are responsible for the MRGW scale selection, are smaller
than tendencies due to forcing when the mean flow is less asymmetric.

##  5 Discussion and Conclusions

The underlying mechanism responsible for the excitation of mixed Rossby-
gravity waves (MRGWs) in the tropical atmosphere is poorly understood,
particularly concerning their synoptic scale selection, which cannot be fully
explained by the excitation associated with the tropical heating (Holton,
1972; Kosaka and Matsuda, 2005). To address this significant knowledge gap, we
employ a new forecasting model TIGAR. This model incorporates the MRG mode as
a prognostic variable alongside other waves. Such an approach enables a direct
quantification of various factors contributing to MRGW development within the
model, presenting a distinct advantage over conventional equatorial wave
filtering techniques.

Using this novel framework, we first demonstrate that MRGWs can be exited by
an asymmetric tropical heat source. Then we show that MRGWs are excited by
wave-mean flow interactions provided that the mean flow is asymmetric with
respect to the equator. In the performed simulations, waves initiating the
interactions are generated by a heat source at the equator. The excited MRGWs
have the maximal energy at the zonal wavenumbers $k=4-5$, in agreement with
observations (e.g. Yang and Hoskins, 2017; Stephan et al., 2021). The non-
linear wave-mean flow interactions as an excitation and scale selection
mechanism is further confirmed through the analysis of the spectral energy
tendencies. Depending on the prescribed mean fluid depth $D$ and the time
stage of the process, the wave-mean flow interactions are mainly due to waves
with low meridional wave index (e.g. for $D=400$ m, the $n=1$ Rossby mode is
the main contributor to the MRGW growth).

Figure 13: Absolute spectral energy tendencies of the MRGW as a function of
zonal wavenumber and time in the non-linear ASYMFOR-JJA simulation with the
zonal mean zonal wind of JJA 1993.

The background flow itself can potentially determine the selected zonal
wavenumber of the excited MRG waves through the instability growth associated
with linearized dynamics. This hypothesis is verified by studying the
eigenmodes of the system linearized about the pertinent background state. The
background eigenmodes are determined following the algorithm suggested by
Kasahara (1980). We find unstable background eigenmodes with substantial
projection onto the MRG modes in range of $k=1-17$ and $k=1-19$ for the
idealized and realistic background zonal flows respectively. However, most of
these unstable eigenmodes have a small growth exponent, which varies only
slightly across nearby scales. For this reason, the selected MRGW scales
generally do not coincide with the fastest growing mode of the linearized
system, rather they are determined by the peak in the wave-mean flow
interactions. This peak fluctuates in time due to the change in phases of
interacting waves, which explains the slight movement of the MRGW spectral
peak. For instance, Figure 13 shows the amplitudes of the contributions of the
wave-mean flow interactions to the MRGW tendencies over 10 days of a
simulation with a realistic background flow in JJA 1993. It shows that the
preferred zonal wavenumber, which is a recipient of most of energy from
interactions, changes over time, but remains at zonal wavenumbers between
$k=3$ and $k=6$, which is consistent with the maximal MRGW signal at $k=5$ in
the reported simulations.

The contribution of the forcing to MRGW tendencies is found to be nearly flat
for the smallest wavenumbers, whereas the MRGW tendencies by wave-mean flow
interactions, which involve waves produced by the forcing, fluctuate over a
small range of $k$ as the phases of these waves change in time. Should the
forcing have a peak zonal wavenumber in the MRGW absolute tendency spectra,
wave-mean flow interactions would also peak at the same $k$, therefore the
scale of the generated MRGWs would be completely predictable. The peak scale
of the MRGW excitation is likely sensitive to waves in the background, which
is a potential subject of a future study.

Two earlier studies of non-linear excitation mechanism of MRGWs require
discussion. First, Itoh and Ghil (1988) showed that asymmetric lateral
boundary forcing of the thermodynamic equation reinforces asymmetric tropical
modes. They suggested that there may be enhanced MRGW activity in the summer
and winter season when the midlatitude westerlies in the winter hemisphere are
intensified. They proposed wave-CISK MRGW excitation mechanism, which
introduces asymmetry through boundary forcing and, therefore, requires certain
properties of extra-tropical flows to maintain. Here we demonstrate a
different mechanism, where MRGW growth is promoted by wave-mean flow
interactions that are local to the tropics and do not require any
extratropical input.

Second, we could not find a strong evidence of the excitation of MRGWs by
wave-wave interactions suggested by Raupp and Silva Dias (2005). In their
simulations with an equatorial $\beta$-plane shallow water model, the initial
state was prepared using an off-equatorial stationary heat source and they
observed strong MRGW development, provided IG waves were suppressed in the
model. In our experiments, an asymmetric heat source leads to the growth of
MRGWs directly regardless of the presence of the IG waves. Moreover,
suppressing IG modes in our simulations does not significantly alter the
contribution of wave-wave interactions to MRGW tendencies (not shown); i.e.
wave-wave interactions are negligible compared to the wave-mean flow
interactions and the excitation process by the external forcing.

Although our simulations are performed in a simplified model, quantitatively
they are consistent with the MRGW energy estimates in reanalysis data. For
instance, an intercomparison of the MRGW energy levels in four re-analyses
datasets by Žagar et al. (2009) showed that they comprise between 4% and 15%
of the non-Rossby wave flow (IG+MRG+Kelvin wave). These percentages roughly
correspond to values obtained in simulations with a realistic background flow
in the present study. In contrast to its energy, the share of the MRGW signal
in intraseasonal variance is large Stephan et al. (2021).

The atmosphere is never in a symmetric zonal mean state. We demonstrated that
the level of asymmetry, defined by a new spectral asymmetry measure (Equation
12), determines the amplitude of the MRGW signal. In particular, the
simulations with the zonal mean state from the JJA season in ERA5 reanalyses
produce a significantly stronger MRGW response in large and synoptic scales
than the ones with less asymmetric zonal wind profiles from the MAM season.
For example, the MRGWs in JJA simulations are roughly 10 times more energetic
compared to MAM simulations when a symmetric forcing is applied.

Similarly to the background flow, the real atmosphere is asymmetric also in
the distribution of convective forcing with the greatest asymmetry during
boreal summer. The quantification of the relative importance of the excitation
by the external forcing (that mimics convection) and by wave-mean flow
interactions in the real atmosphere requires the extension of the present
study to the three-dimensional moist atmosphere.

We speculate that the differences between tropospheric and stratospheric MRGWs
such as the zonal scale with the largest signal could be explained by wave-
mean flow interactions that may peak at different scales compared to the
simulations run with tropospheric flow profiles. The MRGW filtering in the
operational ECMWF forecasts at the MODES webpage, http://modes.cen.uni-
hamburg.de shows much larger scales and greater amplitudes of the MRGWs in the
middle atmosphere compared to the troposphere. This may be a result of the
larger asymmetry of the zonal mean wind, i.e. an enhanced MRGW generation by
wave-mean flow interactions in the middle atmosphere, however, this is a
subject of a follow-on study.

Finally, we note that the presented results can help to better understand the
role of MRGWs in tropical intrinsic and practical predictability and tropical-
extratropical interactions. The non-linear excitation process suggests a
limited intrinsic predictability of tropical synoptic-scale disturbances with
a significant meridional circulation component.

## Acknowledgments

This paper is a contribution to the project W6 (Spectral Energy Fluxes by
Wave-Wave interactions) of the Collaborative Research Centre TRR 181 ”Energy
Transfers in Atmosphere and Ocean” funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer
274762653. We would like to thank Rolando Garcia for the discussion and the
useful comments that aided this manuscript. We further acknowledge the
discussions with Frank Lunkeit, who contributed to the improvement of the
paper.

## References

  * (1)
  * Alexander et al. (2008) Alexander, S., Tsuda, T., Kawatani, Y. and Takahashi, M. (2008), ‘Global distribution of atmospheric waves in the equatorial upper troposphere and lower stratosphere: Cosmic observations of wave mean flow interactions’, J. Geophys. Res.: Atmospheres 113, D24.
  * Bergman and Salby (1994) Bergman, J. W. and Salby, M. L. (1994), ‘Equatorial wave activity derived from fluctuations in observed convection’, J. Atmos. Sci. 51, 3791–3806.
  * Boyd and Zhou (2008) Boyd, J. P. and Zhou, C. (2008), ‘Kelvin waves in the nonlinear shallow water equations on the sphere: nonlinear travelling waves and the corner wave bifurcation’, J. Fluid Mech. 617, 187–205.
  * Chang (1970) Chang, C.-P. (1970), ‘Westward propagating cloud patterns in the tropical pacific as seen from time-composite satellite photographs’, J. Atmos. Sci. 27, 133–138.
  * Dickinson and Molinari (2002) Dickinson, M. and Molinari, J. (2002), ‘Mixed rossby–gravity waves and western pacific tropical cyclogenesis. part i: Synoptic evolution’, J. Atmos. Sci. 59, 2183–2196.
  * Gelb and Gleeson (2001) Gelb, A. and Gleeson, J. P. (2001), ‘Spectral viscosity for shallow water equations in spherical geometry’, Mon. Wea. Rev. 129, 2346–2360.
  * Hayashi (1970) Hayashi, Y. (1970), ‘A theory of large-scale equatorial waves generated by condensation heat and accelerating’, J. Meteor. Soc. Japan. Ser. II 48, 140–160.
  * Hendon and Liebmann (1991) Hendon, H. H. and Liebmann, B. (1991), ‘The structure and annual variation of antisymmetric fluctuations of tropical convection and their association with Rossby–gravity waves’, J. Atmos. Sci. 48, 2127–2140.
  * Hersbach et al. (2020) Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D. et al. (2020), ‘The ERA5 global reanalysis’, Quart. J. Roy. Meteor. Soc. 146, 1999–2049.
  * Holton (1972) Holton, J. R. (1972), ‘Waves in the equatorial stratosphere generated by tropospheric heat sources’, J. Atmos. Sci. 29, 368–375.
  * Hoskins and Yang (2021) Hoskins, B. and Yang, G.-Y. (2021), ‘The detailed dynamics of the hadley cell. part ii: December–february’, J. Climate 34, 805–823.
  * Itoh and Ghil (1988) Itoh, H. and Ghil, M. (1988), ‘The generation mechanism of mixed Rossby-gravity waves in the equatorial troposphere’, J. Atmos. Sci. 45, 585–604.
  * Kasahara (1977) Kasahara, A. (1977), ‘Numerical integration of the global barotropic primitive equations with Hough harmonic expansions’, J. Atmos. Sci. 34, 687–701.
  * Kasahara (1978) Kasahara, A. (1978), ‘Further studies on a spectral model of the global barotropic primitive equations with Hough harmonic expansions’, J. Atmos. Sci. 35, 2043–2051.
  * Kasahara (1980) Kasahara, A. (1980), ‘Effect of zonal flows on the free oscillations of a barotropic atmosphere’, J. Atmos. Sci. 37.
  * Kasahara and da Silva Dias (1986) Kasahara, A. and da Silva Dias, P. L. (1986), ‘Response of planetary waves to stationary tropical heating in a global atmosphere with meridional and vertical shear’, J. Atmos. Sci. 43, 1893–1912.
  * Kawatani et al. (2010) Kawatani, Y., Watanabe, S., Sato, K., Dunkerton, T. J., Miyahara, S. and Takahashi, M. (2010), ‘The roles of equatorial trapped waves and internal inertia–gravity waves in driving the quasi-biennial oscillation. Part I: Zonal mean wave forcing’, J. Atmos. Sci. 67, 963–980.
  * Kosaka and Matsuda (2005) Kosaka, Y. and Matsuda, Y. (2005), ‘Roles of Rossby and gravity waves on circulation associated with tropical and subtropical heating’, J. Meteor. Soc. Japan. Ser. II 83, 481–498.
  * Kuo (1949) Kuo, H.-l. (1949), ‘Dynamic instability of two-dimensional nondivergent flow in a barotropic atmosphere’, J. Atmos. Sci. 6, 105–122.
  * Lindzen and Holton (1968) Lindzen, R. S. and Holton, J. R. (1968), ‘A theory of the quasi-biennial oscillation’, J. Atmos. Sci. 25, 1095–1107.
  * Longuet-Higgins (1968) Longuet-Higgins, M. S. (1968), ‘The eigenfunctions of Laplace’s tidal equation over a sphere’, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 262, 511–607.
  * Mak (1969) Mak, M.-K. (1969), ‘Laterally driven stochastic motions in the tropics’, J. Atmos. Sci. 26, 41–64.
  * Matsuno (1966) Matsuno, T. (1966), ‘Quasi-geostrophic motions in the equatorial area’, J. Meteor. Soc. Japan. Ser. II 44, 25–43.
  * McFarlane et al. (2007) McFarlane, S. A., Mather, J. H. and Ackerman, T. P. (2007), ‘Analysis of tropical radiative heating profiles: A comparison of models and observations’, Journal of Geophysical Research: Atmospheres 112.
  * Raupp and Silva Dias (2005) Raupp, C. F. and Silva Dias, P. L. (2005), ‘Excitation mechanism of mixed Rossby–gravity waves in the equatorial atmosphere: Role of the nonlinear interactions among equatorial waves’, J. Atmos. Sci. 62, 1446–1462.
  * Stephan et al. (2021) Stephan, C. C., Žagar, N. and Shepherd, T. G. (2021), ‘Waves and coherent flows in the tropical atmosphere: New opportunities, old challenges’, Quart. J. Roy. Meteor. Soc. 147, 2597–2624.
  * Takasuka et al. (2019) Takasuka, D., Satoh, M. and Yokoi, S. (2019), ‘Observational evidence of mixed rossby-gravity waves as a driving force for the mjo convective initiation and propagation’, Geophys. Res. Lett. 46, 5546–5555.
  * Takayabu (1994) Takayabu, Y. N. (1994), ‘Large-scale cloud disturbances associated with equatorial waves Part I: Spectral features of the cloud disturbances’, J. Meteor. Soc. Japan. Ser. II 72, 433–449.
  * Takayabu and Nitta (1993) Takayabu, Y. N. and Nitta, T. (1993), ‘3-5 day-period disturbances coupled with convection over the tropical Pacific Ocean’, J. Meteor. Soc. Japan. Ser. II 71, 221–246.
  * Vasylkevych and Žagar (2021) Vasylkevych, S. and Žagar, N. (2021), ‘A high-accuracy global prognostic model for the simulation of Rossby and gravity wave dynamics’, Quart. J. Roy. Meteor. Soc. 147, 1989–2007.
  * Wallace and Chang (1972) Wallace, J. and Chang, L. (1972), ‘On the application of satellite data on cloud brightness to the study of tropical wave disturbances’, J. Atmos. Sci. 29, 1400–1403.
  * Wheeler and Kiladis (1999) Wheeler, M. and Kiladis, G. N. (1999), ‘Convectively coupled equatorial waves: Analysis of clouds and temperature in the wavenumber–frequency domain’, J. Atmos. Sci. 56, 374–399.
  * Wilson and Mak (1984) Wilson, J. D. and Mak, M. (1984), ‘Tropical response to lateral forcing with a latitudinally and zonally nonuniform basic state’, J. Atmos. Sci. 41, 1187–1201.
  * Yanai and Hayashi (1969) Yanai, M. and Hayashi, Y. (1969), ‘Large-scale equatorial waves penetrating from the upper troposphere into the. lower stratosphere’, J. Meteor. Soc. Japan. Ser. II 47, 167–182.
  * Yanai and Maruyama (1966) Yanai, M. and Maruyama, T. (1966), ‘Stratospheric wave disturbances propagating over the equatorial pacific’, J. Meteor. Soc. Japan. Ser. II 44, 291–294.
  * Yanai and Murakami (1970) Yanai, M. and Murakami, M. (1970), ‘A further study of tropical wave disturbances by the use of spectrum analysis’, J. Meteor. Soc. Japan. Ser. II 48, 185–197.
  * Yang and Hoskins (2017) Yang, G.-Y. and Hoskins, B. J. (2017), ‘The equivalent barotropic structure of waves in the tropical atmosphere in the western hemisphere’, J. Atmos. Sci. 74, 1689–1704.
  * Žagar et al. (2015) Žagar, N., Kasahara, A., Terasaki, K., Tribbia, J. and Tanaka, H. (2015), ‘Normal-mode function representation of global 3-D data sets: open-access software for the atmospheric research community’, Geoscientific Model Development 8, 1169–1195.
  * Žagar et al. (2022) Žagar, N., Lunkeit, F., Sielmann, F. and Xiao, W. (2022), ‘Three-dimensional structure of the equatorial kelvin wave: Vertical structure functions, equivalent depths, and frequency and wavenumber spectra’, J. Climate 35, 2209–2230.
  * Žagar et al. (2009) Žagar, N., Tribbia, J., Anderson, J. and Raeder, K. (2009), ‘Uncertainties of estimates of inertia–gravity energy in the atmosphere. Part II: Large-scale equatorial waves’, Mon. Wea. Rev. 137, 3858–3873.
  * Zangvil and Yanai (1980) Zangvil, A. and Yanai, M. (1980), ‘Upper tropospheric waves in the tropics. Part I: Dynamical analysis in the wavenumber-frequency domain’, J. Atmos. Sci. 37, 283–298.
  * Zhang (1993) Zhang, C. (1993), ‘Laterally forced equatorial perturbations in a linear model. Part II: Mobile forcing’, J. Atmos. Sci. 50, 807–821.