# Tidally Forced Planetary Waves in the Tachocline of Solar-like Stars

Gerrit M. Horstmann Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden-
Rossendorf, Bautzner Landstrasse 400, 01328 Dresden, Germany  George
Mamatsashvili Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden-
Rossendorf, Bautzner Landstrasse 400, 01328 Dresden, Germany  Ilia State
University, Cholokashvili ave. 5/3, Tbilisi, Georgia  Abastumani Astrophysical
Observatory, Mount Kanobili, Georgia  André Giesecke Institute of Fluid
Dynamics, Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstrasse 400,
01328 Dresden, Germany  Teimuraz V. Zaqarashvili Institute of Physics, IGAM,
University of Graz, Universitätsplatz 5, 8010, Graz, Austria  Ilia State
University, Cholokashvili ave. 5/3, Tbilisi, Georgia  Abastumani Astrophysical
Observatory, Mount Kanobili, Georgia  Frank Stefani Institute of Fluid
Dynamics, Helmholtz-Zentrum Dresden-Rossendorf, Bautzner Landstrasse 400,
01328 Dresden, Germany

###### Abstract

Can atmospheric waves in planet-hosting solar-like stars substantially
resonate to tidal forcing? Substantially at a level of impacting the space
weather or even of being dynamo-relevant? In particular, low-frequency Rossby
waves, which have been detected in the solar near-surface layers, are
predestined at responding to sunspot cycle-scale perturbations. In this paper,
we seek to address these questions as we formulate a forced wave model for the
tachocline layer, which is widely considered as the birthplace of several
magnetohydrodynamic planetary waves, i.e., Rossby, inertia-gravity (Poincaré),
Kelvin, Alfvén and gravity waves. The tachocline is modeled as a shallow
plasma atmosphere with an effective free surface on top that we describe
within the Cartesian $\beta$-plane approximation. As a novelty to former
studies, we equip the governing equations with a conservative tidal potential
and a linear friction law to account for dissipation. We combine the
linearized governing equations to one decoupled wave equation, which
facilitates an easily approachable analysis. Analytical results are presented
and discussed within several interesting free, damped and forced wave limits
for both mid-latitude and equatorially trapped waves. For the idealized case
of a single tide generating body following a circular orbit, we derive an
explicit analytic solution that we apply to our Sun for estimating leading-
order responses to Jupiter. Our analysis reveals that Rossby waves resonating
to low-frequency perturbations can potentially reach considerable velocity
amplitudes in the order of $10^{1}\,$–$\,10^{2}\,{\rm cm}\,{\rm s}^{-1}$,
which, however, strongly rely on the yet unknown total dissipation.

Sun: interior – Sun: magnetic fields – Sun: oscillations

††software:

##  1 Introduction

For more than 100 years it has been known that the global weather and climate
of our earth is decisively impacted by atmospheric planetary waves, in which
large-scale Rossby waves occupy the most prominent role (Pedlosky, 1987). But
only in the last 20 years strong evidence has accumulated that Rossby waves
play similar vital roles in various astrophysical objects, such as planets,
e.g., Jupiter (Li et al., 2006) or Saturn (Read et al., 2006), accretion discs
(Lyra & Umurhan, 2019) and, most importantly, in the Sun and other solar-like
stars, see Zaqarashvili et al. (2021) for an in-depth review. With regard to
our Sun, it is known today that Rossby waves are promising candidates for
explaining the solar seasons Dikpati et al. (2017, 2018a), inducing angular
momentum transport (Gizon et al., 2020) and impacting or even causing solar
activity cycles, starting from Riega-type periodicities (Zaqarashvili et al.,
2010a) over Schwabe cycle fluctuations (Raphaldini & Raupp, 2015; Raphaldini
et al., 2019) up to long-term modulations in the order of the Gleisberg cycle
(Zaqarashvili, 2018). Rossby waves may also serve a crucial role for the solar
dynamo (Zaqarashvili et al., 2021, section 5.5) and have even been considered
as a key ingredient for dynamo action; an early idea of a self-exciting Rossby
wave dynamo dates back to Gilman (1968).

For a long time, Rossby waves with respect to the Sun were mostly perceived as
a theoretical concept and a variety of explanations for their hypothetical
emergence were proposed. The breakthrough has only very recently been
obtained: different types of Rossby waves have been detected independently by
employing different methodologies. First, McIntosh et al. (2017) have observed
magnetic Rossby-like waves in the solar atmosphere by tracking coronal bright
points from EUV images. Thereupon, Löptien et al. (2018) and Liang et al.
(2019) have further observed classic Rossby modes in solar near-surface layers
from helioseismic measurements. It is important to note that classic and
magnetic Rossby waves are physically very different and have, among other more
nuanced dissimilarities, exactly opposite phase velocities and opposite group
velocities; see Dikpati et al. (2020) for an excellent introduction to the
physics of solar Rossby waves. Both types of waves can therefore play very
different roles in the solar dynamics and it is today of great concern to
understand the different manifestations of Rossby and other
magnetohydrodynamic planetary waves in the Sun (Zaqarashvili et al., 2021).

Today, two of the most important and as yet largely unresolved questions are
where and how solar Rossby waves might originate. Although Rossby waves have
been observed in the outer solar atmosphere, the shallow tachocline layer is
widely believed to be one of the most promising birthplaces of solar planetary
waves. Gilman (2000) has shown in a pioneering work that the solar tachocline
fluid-layer can be treated, in terms of hydrodynamics, fairly analogously to
the lower atmosphere of the Earth, which is why the well-studied geophysical
shallow water equations, governing Rossby and other classical atmospheric
waves, can be transferred almost one-to-one to the solar tachocline. Since
that time, numerous two-dimensional as well as quasi-3D shallow-water models
have been employed to study the global wave instabilities in the tachocline,
see, e.g., Dikpati & Gilman (2001); Schecter et al. (2001); Gilman & Dikpati
(2002); Zaqarashvili et al. (2007, 2009, 2010b); Raphaldini & Raupp (2015);
Klimachkov & Petrosyan (2017); Dikpati et al. (2018a). These studies have
identified a number of different possible causes for the development of Rossby
waves, including different shear instabilities associated to differential
rotation, thermal forcing and nonlinear wave-wave interactions.

One natural creation mechanism, which has gained a special importance in
recent years concerning the solar dynamo, has been disregarded so far: the
possible wave excitation by tidal forcing. It was emphasized in a number of
studies (Stefani et al., 2016, 2018, 2019, 2021) that the combined tidal
action of the planets Earth, Venus and Jupiter might play a significant role
in the synchronization process of the solar dynamo. Although the responding
tidal height is only in the order of $1\,{\rm mm}$, energetically equivalent
velocities can reach up to $1\,{\rm m}\,{\rm s}^{-1}$ due to the high
gravitational acceleration in the tachocline, which could indeed be dynamo-
relevant. While in those works the synchronization mechanism was considered to
rely only on the entrainment of the $\alpha$-effect caused by the Taylor
instability (Weber et al., 2015), the thrilling question has arisen of whether
Rossby and other planetary waves could take on a similar facilitating role.
Can Rossby waves possibly intensify the tidal action, or, in other words, may
they serve as kind of a resonance ground for tidal excitations? Slow magnetic
Rossby waves can have periods in the order of the solar cycle, which is
remarkably close to the 11.07 year period visible in the spring-tide envelope
curve of Jupiter’s, Venus’ and Earth’s tidal potentials (Okhlopkov, 2016). But
also in other respects Rossby waves entail excellent premises on which to act
on the solar dynamo. First, Rossby waves can have a net kinetic helicity
(Dikpati & Gilman, 2001; Gilman & Dikpati, 2002) letting them participate in
the $\alpha$-effect. Second, tachoclinic oscillations can further sensitively
affect magnetic flux tube instabilities since very small variations of the
superadiabacity $\delta$ (stratification of specific entropy) in the order of
$\delta\sim 10^{-4}$ to $\delta\sim 10^{-5}$ considerably alter the magnetic
storage capacity (Ferriz-Mas, 1996; Abreu et al., 2012; Charbonneau, 2022).
Motivated by all these auspicious premises, we devote this study to tidally
forced magnetic planetary waves and present a first theoretical “shallow-
water” model, which can account for arbitrary tidal potentials and, as a
second novelty, further entails Rayleigh friction permitting the study of
damped wave mechanics.

The paper is organized as follows. In Section 2 we formulate the forced wave
model on the local Cartesian $\beta$-plane and rearrange the governing
equations into one decoupled wave equation for the latitude velocity, which
covers the entire wave physics and makes the analysis easily accessible. In
Section 3 we present analytical results for both mid-latitude and equatorially
trapped waves. We start to recover the known free wave problems and gradually
increase in complexity all the way from different freely damped limits up to
the full forced wave problem, which is solved for a single tide generating
body. In Section 4, the solutions are finally applied to the particular
scenario of tachoclinic waves in our Sun forced by Jupiter to estimate
resonant velocity amplitudes.

##  2 Formulation of the shallow-water model

###  2.1 Basic considerations

Shallow water formulations have been widely employed in the last century to
describe several types of planetary waves both in the atmosphere and oceans of
rotating planets (Longuet-Higgins, 1964; Pedlosky, 1987)—with Rossby waves
perhaps being the most prominent member. In its simplest view, an ocean can be
intuitively modeled as a one-layer system of fluid sandwiched between a flat
rigid bottom and a freely movable surface placed on top. If we focus on large-
scale planetary waves, this layer can be considered as ‘shallow’ in the sense
that characteristic wave lengths are much larger than the average water
depths. This has the consequence that the flow is of quasi-two-dimensional
nature, i.e., the velocities do not alter significantly with altitude and the
change of momentum caused by vertical velocities and vertical Coriolis force
components can be neglected in the momentum conservation law. Such shallow
water approaches are very intelligible and simplify the mathematical analysis
enormously. Therefore, it is all the more gratifying that planetary waves
evolving in tachoclines of solar-like stars can be described in very similar
ways. It was pointed out in the pioneering paper by Gilman (2000) that the
geophysical shallow water equations can be transferred almost one-to-one to
the solar tachocline. In a manner of speaking, the tachocline can be treated
as a sort of plasma ocean as schematized in Figure 1.

Figure 1: Schematic drawing of the tachocline modeled as a plasma ocean with
the average depth $H_{0}$ bounded by a rigid bottom and an imaginary free
interface at the top. The interface feels an effective gravity $\bm{g}$
arising from the fact that the convection layer, which is adiabatic or even
superadiabatic, does not offer any buoyancy resistance to fluid volumes coming
from below, whereas fluid volumes within the subadiabatically stratified
tachocline are subject to negative buoyancy.

The tachocline layer with the depth $H_{0}$, which itself consists of a
radiative part and an overshoot layer (not shown here), is the transition
region between the radiative interior and the outer convection zone. The
convection zone is of adiabatic or even superadiabatic nature, whereas the
stably stratified tachocline is subadiabatic. Hence, warmer fluid volumes
entering the convection zone tend to rise (or at least do not experience any
buoyancy resistance in the adiabatic case), while warmer fluid volumes tend to
sink in the tachocline. As a result of this behavior, the transition region
between the tachocline and the convection zone can be replaced by an imaginary
interface, which experiences some effective gravity in the same way as ocean-
atmosphere interfaces are exposed to standard gravity. The effective gravity
force is proportional to the fractional difference between the actual and
adiabatic temperature gradients $\sim|\nabla-\nabla_{\rm ad}|$, taking values
of $10^{-4}-10^{-6}$ in the upper overshoot part of the tachocline and up to
$10^{-1}$ in the lower radiative sublayer (Dikpati & Gilman, 2001). These
gradients yield very different effective gravity constants of $g=0.05-5\,{\rm
cm}\,{\rm s}^{-2}$ (overshoot part) and $g=500-1.5\cdot 10^{4}\,{\rm cm}\,{\rm
s}^{-2}$ (radiative sublayer) (Schecter et al., 2001).

The shallow water equations of the plasma ocean model can be applied
separately to both parts of the tachocline (but not both at the same time)
such that the natural frequencies of planetary waves vary by orders of
magnitude among the sublayers. The density is approximated to be constant and
the lower interface to be stationary since the higher dense radiative interior
appears much more rigid than the tachocline layer above. Note that there are
enhanced two-layer and even continuous descriptions (Hunter, 2015; Petrosyan
et al., 2020; Fedotova et al., 2021), the here adopted one-layer system,
however, is well established and contains the essential wave dynamics,
particularly with regard to Rossby waves, which allowed to make significant
contributions to various astrophysical settings (Zaqarashvili et al., 2021).

###  2.2 Governing equations

We consider a thin layer of an ideally conductive, incompressible and inviscid
fluid on a sphere within the rotating frame of reference. The angular velocity
$\bm{\Omega}_{0}$ as well as the density $\rho$ are set to be constant. The
layer of height $H_{0}$ is further traversed by a magnetic field $\bm{B}$,
which is nearly horizontal in the solar tachocline. As the key new feature, we
also take into account external tidal forces, which can be uniquely expressed
through the gradient of a conservative tidal potential $V$. Since the inviscid
models lack any dissipation source, the inclusion of an external force gives
rise to the problem of singularities appearing under resonance conditions. In
geophysical flows it is common to include linear friction terms, usually
Rayleigh friction and/or Newtonian cooling (see, e.g., Gill (1980); Mofjeld
(1981); Yamagata & Philander (1985); Wu et al. (2001)), which are indicative
in the framework of linear waves. Here, we include Rayleigh’s friction law,
stating that, in the leading order, the effective drag is proportional to the
flow velocity and some a priori unknown damping constant $\lambda$. Although
viscous and thermal dissipation in the tachocline is highly nonlinear on
smaller scales, the linear friction law can nevertheless be an appropriate
approximation for many long-period and large-scale planetary waves—at least
when the wave amplitudes are sufficiently small.

By equipping the classical magnetohydrodynamic shallow water equations
proposed by Gilman (2000) with tidal and friction forces, our governing
equations can be expressed in the following way:

$\displaystyle\frac{\partial\bm{u}}{\partial
t}+(\bm{u}\cdot\bm{\nabla})\bm{u}+2\bm{\Omega}_{0}\times\bm{u}+\lambda\bm{u}$
$\displaystyle=-g\bm{\nabla}H+\frac{1}{4\pi\rho}(\bm{B}\cdot\bm{\nabla})\bm{B}-\bm{\nabla}V,$
(1) $\displaystyle\frac{\partial\bm{B}}{\partial
t}+(\bm{u}\cdot\bm{\nabla})\bm{B}=(\bm{B}\cdot\bm{\nabla})\bm{u},$ (2)
$\displaystyle\frac{\partial H}{\partial t}+\bm{\nabla}\cdot(H\bm{u})=0,$ (3)
$\displaystyle\bm{\nabla}\cdot(H\bm{B})=0,$ (4)

where $\bm{u}$, $\bm{\Omega}_{0}$ and $\bm{B}$ are the horizontal velocity,
angular velocity of the rotating star and magnetic field, $H$ is the total
layer height, $\rho$ is the fluid density, $g$ is effective gravitational
acceleration and $\bm{\nabla}$ denotes here the purely horizontal gradient.
Equation (1) is the Euler equation including the Coriolis force, the Lorentz
force, gravity and Rayleigh friction. Equation (2) is the induction equation
of the magnetic field in the limit of ideally conducting fluids (high magnetic
Reynolds number limit) and (3) is the shallow water version of the continuity
equation. Equation (4) finally ensures that the magnetic field is divergence-
free on the condition that $\bm{B}$ remains parallel to the upper free
surface.

###  2.3 $\beta$-plane approximation

The inclusion of tidal forcing drastically enriches the physical complexity.
For this study, we are mainly interested in finding analytical solutions
rather than conducting simulations to understand the incoming physics in a
fundamental way. Therefore, in order to keep the analysis as simple as
possible, we are going to study the problem in a simpler Cartesian coordinate
system within the framework of the so-called $\beta$-plane approximation, see
Figure 2. The $\beta$-plane can properly

Figure 2: Local Cartesian coordinate system defined at a fixed latitude
$\phi_{0}$ in the non-inertial frame of reference on the tachocline surface of
a sphere constantly rotating with angular velocity $\Omega_{0}$. The $x$-axis
points eastwards (longitude), the $y$-axis northwards (latitude) and the
$z$-axis radially outwards from center (altitude).

describe non-equatorial planetary waves if the wave length is sufficiently
smaller than the size of the sphere, and is perfectly valid for waves trapped
at the equator, where they have exactly the same dispersion relation as in the
analogues spherical coordinate system (Pedlosky, 1987). The local Cartesian
coordinate system is fixed at the latitude $\phi_{0}$, at which we can
evaluate the Coriolis force

$\displaystyle 2\bm{\Omega}\times\bm{u}=2\Omega_{0}\begin{pmatrix}0\\\
\cos\phi_{0}\\\ \sin\phi_{0}\end{pmatrix}\times\begin{pmatrix}u\\\ v\\\
w\end{pmatrix}$ $\displaystyle=\ $ $\displaystyle
2\Omega_{0}\begin{pmatrix}v\sin\phi_{0}-w\cos\phi_{0}\\\ -u\sin\phi_{0}\\\
u\cos\phi_{0}\end{pmatrix}\approx f\begin{pmatrix}v\\\ -u\\\ 0\end{pmatrix},$
(5)

where $u,v,w$ are the longitude, latitude and altitude velocity components and
$f=2\Omega_{0}\sin(\phi_{0})$ is the Coriolis parameter. The shallow water
approximation was applied at the last step, demanding that both the altitude
velocity and Coriolis force components are negligibly small,
$w,u\cos(\phi_{0})\approx 0$. It was the great pioneering contribution of
Carl-Gustaf Arvid Rossby (1939) to locally expand $f$ on the Cartesian plane
for small $\phi_{0}$ variations as

$\displaystyle f\approx f_{0}+\beta y,$ (6)

with

$\displaystyle\beta=\frac{{\rm d}f}{{\rm
d}y}=\frac{2\Omega_{0}}{R_{0}}\cos(\phi_{0}),$ (7)

today referred to as the Rossby parameter. Keeping only the first order
$f=f_{0}=2\Omega_{0}\sin(\phi_{0})$ leads to the $f$-plane approximation,
involving different types of magneto-gravity waves (also magneto-Poincaré
waves) and purely magnetohydrodynamic Alfvén waves (Schecter et al., 2001), or
also magnetostrophic waves evolving under the presence of vertical magnetic
fields (Klimachkov & Petrosyan, 2016). Accordingly, if one further considers
the leading order latitudinal variation of the Coriolis force (second term),
we arrive at the so called $\beta$-plane approximation comprising low-
frequency magneto-Rossby waves in addition to the $f$-plane waves.

As the last step of approximation, we restrict our analysis to the linear
problem, which allows us to derive analytical solutions at the little expense
of being limited to small-amplitude waves. By projecting the governing
equations (1) - (4) onto the $\beta$-plane and linearizing, we get the
equation system

$\displaystyle\frac{\partial u}{\partial t}=-g\frac{\partial\eta}{\partial
x}+fv+\frac{B_{0}}{4\pi\rho}\frac{\partial b_{x}}{\partial x}-\frac{\partial
V}{\partial x}-\lambda u,$ (8) $\displaystyle\frac{\partial v}{\partial
t}=-g\frac{\partial\eta}{\partial y}-fu+\frac{B_{0}}{4\pi\rho}\frac{\partial
b_{y}}{\partial x}-\frac{\partial V}{\partial y}-\lambda v,$ (9)
$\displaystyle\frac{\partial b_{x}}{\partial t}=B_{0}\frac{\partial
u}{\partial x},\ \ \frac{\partial b_{y}}{\partial t}=B_{0}\frac{\partial
v}{\partial x},$ (10) $\displaystyle\frac{\partial\eta}{\partial
t}+H_{0}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial
y}\right)=0,$ (11)

for the individual vector components. Here, $u,v,b_{x}$ and $b_{y}$ are the
perturbational field quantities and $\eta=H-H_{0}$ denotes the perturbed layer
thickness ($\equiv$ wave amplitude). The equations were perturbed with respect
to a uniform toroidal magnetic field $B_{0}$, which is dominant in the solar
tachocline, although generally latitude-dependent. The set of equations is
analogues to the free-wave problem analyzed by Zaqarashvili et al. (2007),
with the sole difference that two additional terms appear in the momentum
equations (8) and (9), addressing tidal forcing and damping.

###  2.4 Decoupled wave equation

Although the linearized $\beta$-plane approximation is by far the most
accessible description of magneto-planetary waves, the associated wave physics
is still very rich. In order to simplify the following analyses and to unify
the free wave physics discussed in preceding studies, it is of interest to
find a generalized wave equation, which allows the straightforward calculation
of dispersion relations in all limiting cases of interest. We found that a
fully decoupled wave equation can be derived only for the latitude velocity
$v$, but not for all other field variables $u,\eta,b_{x}$ or $b_{y}$. After
several algebraic transformations, see Appendix A, the set of Equations (8) -
(11) can be rearranged into the following wave equation:

$\square^{2}_{v_{A}}v-C_{0}^{2}\square_{v_{A}}\Delta
v+f^{2}\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\beta\frac{\partial}{\partial x}\frac{\partial v}{\partial
t}+\frac{\partial}{\partial t}\frac{\partial}{\partial
y}\square_{v_{A}}V-f\frac{\partial}{\partial x}\frac{\partial^{2}V}{\partial
t^{2}}+\lambda\frac{\partial}{\partial y}\frac{\partial^{2}V}{\partial
t^{2}}+2\lambda\frac{\partial}{\partial t}\square_{v_{A}}v-\lambda
C_{0}^{2}\Delta\frac{\partial v}{\partial
t}+\lambda^{2}\frac{\partial^{2}v}{\partial t^{2}}=0,$ (12)

where $C_{0}=\sqrt{gH_{0}}$ and $v_{A}=B_{0}/\sqrt{4\pi\rho}$ are gravity wave
and Alfvén velocities, and
$\square_{v_{A}}:=\partial_{t}^{2}-v_{A}^{2}\partial_{x}^{2}$ denotes the
d’Alembert operator with respect to Alfvén waves. Equation (12) serves as the
only governing equation throughout the rest of the paper.

##  3 Analytical results

We shall now deduce and discuss tachoclinic wave dynamics within several
different wave limits, comprising unbounded ($f\approx f_{0}$) and
equatorially trapped ($f\approx\beta y$), hydrodynamic ($v_{A}=0$) and
magnetic, inviscid ($\lambda=0$) and damped, as well as free ($V=0$) and
forced waves. First, we reiterate some basic result on the free wave problem
in order to then increase progressively in complexity towards the full tidally
forced wave problem.

###  3.1 Free wave dynamics

When neglecting the tidal potential, $V=0$, and damping, $\lambda=0$, Equation
(12) drastically simplifies to the free wave equation

$\square^{2}_{v_{A}}v-C_{0}^{2}\square_{v_{A}}\Delta
v+f^{2}\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\beta\frac{\partial}{\partial x}\frac{\partial v}{\partial
t}=0.$ (13)

If one further considers non-magnetic waves $v_{A}=0$, (13) transforms into

$\frac{\partial^{3}v}{\partial t^{3}}+f^{2}\frac{\partial v}{\partial
t}-C_{0}^{2}\frac{\partial}{\partial t}\Delta v-C_{0}^{2}\beta\frac{\partial
v}{\partial x}=0,$ (14)

which is the classic planetary wave equation describing Rossby, Poincaré,
Kelvin and gravity waves (Pedlosky, 1987). In the following, we want to
shortly recapitulate the wave physics captured in Equation (13) within two
different limits. First, we focus on non-equatorial waves in the latitude
range $30^{\circ}\lesssim\phi_{0}\lesssim 60^{\circ}$, where we readily find
$f_{0}\gg\beta y$ allowing to set $f\approx f_{0}$. In the vicinity of the
equator, we observe $f_{0}\ll\beta y$ in contrast, such that one must keep the
$y$-dependent term of the Coriolis parameter $f\approx\beta y$, which results
in a nonlinear wave equation coming along with fundamentally different wave
dynamics.

####  3.1.1 Non-equatorial waves

Applying $f=f_{0}$ to Equation (13) and inserting a simple Fourier ansatz of
the form

$v=v_{0}\exp(ik_{x}x+ik_{y}y-i\omega t),$ (15)

where $v_{0}$ is an arbitrary wave amplitude, $k_{x}$ and $k_{y}$ are the
longitudinal and latitudinal wave numbers and $\omega$ is the angular
frequency, we can straightforwardly deduce the following fourth-order
dispersion relation

$\displaystyle\omega^{4}-\left[2k_{x}^{2}v_{A}^{2}+f_{0}^{2}+C_{0}^{2}(k_{x}^{2}+k_{y}^{2})\right]\omega^{2}-C_{0}^{2}k_{x}\beta\omega$
$\displaystyle+k_{x}^{2}v_{A}^{2}\left[k_{x}^{2}v_{A}^{2}+C_{0}^{2}(k_{x}^{2}+k_{y}^{2})\right]=0.$
(16)

This dispersion relation was first derived and discussed by Zaqarashvili et
al. (2007). Let us first consider the limit $\beta=0$ corresponding to the
$f$-plane approximation. In this case, only symmetric frequency polynomials
$\omega^{4},\omega^{2},\omega^{0}$ remain in Equation (16), allowing us to
rearrange it into the explicit form

$\displaystyle\omega_{\pm}^{2}=k_{x}^{2}v_{A}^{2}+\frac{f_{0}^{2}}{2}+\frac{C_{0}^{2}}{2}\Bigg{[}\left(k_{x}^{2}+k_{y}^{2}\right)$
$\displaystyle\pm\sqrt{\frac{f_{0}^{4}}{C_{0}^{4}}+\frac{f_{0}^{2}}{C_{0}^{2}}\left[\frac{4k_{x}^{2}v_{A}^{2}}{C_{0}^{2}}+2(k_{x}^{2}+k_{y}^{2})\right]+\left(k_{x}^{2}+k_{y}^{2}\right)^{2}}\Bigg{]}.$
(17)

This dispersion relation was analyzed in depth by Schecter et al. (2001), who
referred the fast $\omega_{+}$ branch to magnetogravity waves and the slow
$\omega_{-}$ modes to Alfvén waves. The $\omega_{+}$ branch can be
approximated as

$\displaystyle\omega^{2}\approx
f_{0}^{2}+C_{0}^{2}(k_{x}^{2}+k_{y}^{2})+2k_{x}^{2}v_{A}^{2},$ (18)

showing that magnetogravity waves are magnetically modified Poincaré waves,
which appear on time scales less than $2\pi/f_{0}$. The Alfvén branch is
intrinsically connected to the magnetic field and disappears for
$v_{A}\longrightarrow 0$. Alfvén waves can become arbitrarily slow at large
length scales ($k_{x},k_{y}\longrightarrow 0$).

Keeping the next order in the expansion (6) ($\beta\neq 0$) gives rise to the
emergence of Rossby waves. They appear in the low frequency limit $\omega\ll
f_{0}$ of Equation (16), where the dispersion relation reads

$\displaystyle\omega^{2}+\frac{R_{R}^{2}\beta
k_{x}}{1+R_{R}^{2}(k_{x}^{2}+k_{y}^{2})}\omega-\frac{R_{R}^{2}v_{A}^{2}k_{x}^{2}(k_{x}^{2}+k_{y}^{2})}{1+R_{R}^{2}(k_{x}^{2}+k_{y}^{2})}=0,$
(19)

here expressed in terms of the Rossby radius $R_{R}=C_{0}/f_{0}$. The Rossby
radius describes the length scale at which Coriolis force-driven inertial
waves become equally significant as buoyancy-driven gravity waves in the
spatiotemporal evolution of linear disturbances. Interestingly, $R_{R}$
strongly varies throughout the tachocline. One can estimate $R_{R}\sim
10^{4}-10^{6}\,{\rm km}$ in the overshoot layer but larger values $R_{R}\sim
10^{7}-10^{8}\,{\rm km}$ in the radiative part of the tachocline at the
latitude $\phi_{0}=30^{\circ}$. The tachocline radius is approximately
$R_{0}\approx 5\cdot 10^{5}\,{\rm km}$. In the radiative sublayer, we can
safely assume $R_{R}\gg R_{0}\gtrsim 1/k_{x},1/k_{y}$, allowing us to simplify
Equation (19) into the pure magneto-Rossby wave dispersion relation

$\displaystyle\omega^{2}+\frac{\beta k_{x}}{k_{x}^{2}+k_{y}^{2}}\omega-
v_{A}^{2}k_{x}^{2}=0.$ (20)

It contains a high frequency solution

$\displaystyle\omega\approx-\frac{k_{x}\beta}{k_{x}^{2}+k_{y}^{2}}$ (21)

as well as a low frequency branch

$\displaystyle\omega\approx\frac{k_{x}v_{A}^{2}(k_{x}^{2}+k_{y}^{2})}{\beta}$
(22)

in the limit of large length scales. Equation (21) is the very classic
dispersion relation of hydrodynamic Rossby waves. The minus sign reveals a
retrograde propagation relative to the rotation of the reference frame.
Equation (22) describes the class of magneto-Rossby waves arising only by the
effect of horizontal magnetic fields $v_{A}\neq 0$. These wave modes are quite
similar to classic Rossby waves, but move in the prograde direction and have
far slower eigen-frequencies, see Zaqarashvili et al. (2007) for more details.
A very illustrative description of both classic and magneto-Rossby waves is
further given in Dikpati et al. (2020).

####  3.1.2 Equatorial waves

In the vicinity of the equator, the Coriolis force approximates to
$f\approx\beta y$ in the leading order. This yields the wave equation

$\square^{2}_{v_{A}}v-C_{0}^{2}\square_{v_{A}}\Delta
v+\beta^{2}y^{2}\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\beta\frac{\partial}{\partial x}\frac{\partial v}{\partial
t}=0$ (23)

as the counterpart to Equation (13). Equation (23) is nonlinear in $y$
therefore it is convenient to first perform a Fourier analysis for the
$x$-part only

$\displaystyle v=v_{y}(y)\exp(ik_{x}x-i\omega t)$ (24)

yielding

$\frac{{\rm d}^{2}v_{y}}{{\rm
d}y^{2}}+\left[\frac{\omega^{2}-k_{x}^{2}(C_{0}^{2}+v_{A}^{2})}{C_{0}^{2}}-\frac{k_{x}\beta\omega}{\omega^{2}-k_{x}^{2}v_{A}^{2}}-\mu^{2}y^{2}\right]v_{y}=0$
(25)

for the latitude-dependent velocity part $v_{y}$. The parameter $\mu$, which
is a reciprocal measure for the wave’s equatorial expansion, is given as

$\displaystyle\mu=\frac{\beta\omega}{C_{0}\sqrt{\omega^{2}-k_{x}^{2}v_{A}^{2}}}.$
(26)

The differential Equation (25) was first derived by Zaqarashvili (2018) and
can be identified as the classic equation of the quantum harmonic oscillator.
It has bounded solutions of the form

$\displaystyle
v_{y}=v_{0}\exp\left(-\frac{|\mu|y^{2}}{2}\right)H_{n}(\sqrt{|\mu|}y),$ (27)

if and only if

$\displaystyle\frac{\omega^{2}-k_{x}^{2}(C_{0}^{2}+v_{A}^{2})}{C_{0}^{2}}-\frac{k_{x}\beta\omega}{\omega^{2}-k_{x}^{2}v_{A}^{2}}=|\mu|(2n+1).$
(28)

Here, $v_{0}$ is some arbitrary amplitude and $H_{n}$ the Hermite polynomial
of order $n$. The solutions are oscillatory inside the latitude interval

$y<\left|\sqrt{\frac{2n+1}{|\mu|}}\right|$ (29)

and exponentially tend to zero outside. Although we have not incorporated any
boundary conditions, we arrive at a natural latitudinal boundary condition
capturing the waves around the equator. The integers $n$ can be identified as
discrete latitudinal wave numbers specifying the number of vortices within the
equatorial band defined by Equation (29).

The constraint (28) defines the dispersion that can be expressed in the
following form

$\displaystyle(\omega^{2}-k_{x}^{2}v_{A}^{2})(\omega^{2}-k_{x}^{2}(C_{0}^{2}+v_{A}^{2}))-k_{x}\beta
C_{0}^{2}\omega$ $\displaystyle=\beta
C_{0}|\omega|\sqrt{\omega^{2}-k_{x}^{2}v_{A}^{2}}(2n+1).$ (30)

Neglecting the magnetic field, $v_{A}=0$, we arrive at the classic geophysical
dispersion relation

$\displaystyle\omega^{2}-k_{x}^{2}C_{0}^{2}-\frac{k_{x}\beta
C_{0}^{2}}{\omega}=\beta C_{0}(2n+1)$ (31)

(Matsuno, 1966) comprising gravity-inertia and Rossby waves. The essential
novelty of Equation (30) compared to the hydrodynamic case (31) is that the
toroidal field creates low frequency cut-off areas at $\omega=\pm k_{x}v_{A}$,
suppressing the low-frequency Rossby modes from the hydrodynamic solution.
This behavior can, however, change drastically when considering inhomogeneous
toroidal magnetic field profiles $B_{x}\sim B_{0}y/R$, giving rise to super-
slow magneto-Rossby waves, which can reach periods up to the order of the
$100\,{\rm yr}$ Gleissberg cycle (Zaqarashvili, 2018).

The dispersion relation (30) can only be solved numerically owing to the
square root on the right-hand side. There are, however, some approximated
solutions for certain interesting limit. Restricting to high frequencies
$\omega\gg k_{x}v_{A}$, Equation (30) approximates to

$\displaystyle\omega^{2}\approx\beta
C_{0}(2n+1)+k_{x}^{2}C_{0}^{2}+2k_{x}^{2}v_{A}^{2},$ (32)

describing magneto-inertia-gravity waves, which propagate always faster than
the hydrodynamic counterparts. By eliminating the high frequency branches from
Equation (30), one finds the low-frequency dispersion relation describing
Rossby waves

$\displaystyle k_{x}^{2}\omega^{2}+\beta k_{x}\omega-
k_{x}^{4}v_{A}^{2}=-\frac{\beta}{C_{0}}|\omega|\sqrt{\omega^{2}-k_{x}^{2}v_{A}^{2}}(2n+1)$
(33)

which also is not amenable to analytical solution. We can only find explicit
solutions in the limit of weak magnetic fields $v_{A}\ll 1$, where Equation
(33) reduces to

$\displaystyle\omega^{2}+\frac{2C_{0}\beta
k_{x}\omega}{C_{0}k_{x}^{2}+\beta(2n+1)}-\frac{v_{A}^{2}k_{x}^{2}[2C_{0}k_{x}^{2}+\beta(2n+1)]}{2[C_{0}k_{x}^{2}+\beta(2n+1)]}=0,$
(34)

which involves again a high-frequency solution

$\displaystyle\omega\approx-\frac{C_{0}\beta
k_{x}}{C_{0}k_{x}^{2}+\beta(2n+1)}$ (35)

and a low-frequency solution

$\displaystyle\omega\approx
k_{x}v_{A}^{2}\left(\frac{k_{x}^{2}}{\beta}+\frac{1}{2}\frac{2n+1}{C_{0}}\right).$
(36)

Equation (35) is the classic dispersion relation of hydrodynamic, equatorially
trapped Rossby waves. Similarly as for the non-equatorial waves, Equation (36)
describes prograde magneto-Rossby waves and reduces exactly to Equation (22)
in the limit of purely longitudinally propagating waves ($k_{y}=0$) and large
gravity velocities $C_{0}$. In case of the solar tachocline, however, these
solutions are truncated by the Alfvén wave branches for magnetic field
strength $>10\,{\rm kG}$ (Zaqarashvili, 2018).

###  3.2 Damped wave dynamics

We have now prepared the groundwork to discuss novel wave solutions under the
effect of damping. Keeping the $\lambda$-dependent terms in Equation (12) but
still neglecting the forcing potential $V$, we arrive at the wave equation

$\displaystyle\square^{2}_{v_{A}}v-C_{0}^{2}\square_{v_{A}}\Delta
v+(f^{2}+\lambda^{2})\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\beta\frac{\partial}{\partial x}\frac{\partial v}{\partial t}$
$\displaystyle+2\lambda\frac{\partial}{\partial t}\square_{v_{A}}v-\lambda
C_{0}^{2}\Delta\frac{\partial v}{\partial t}=0.$ (37)

####  3.2.1 Non-equatorial waves

For non-equatorial waves, we can again perform a simple Fourier analysis
$v\sim\exp(ik_{x}x+ik_{y}y-i\omega t)$, yielding the complex dispersion
relation

$\displaystyle\omega^{4}+2i\lambda\omega^{3}-\left[2k_{x}^{2}v_{A}^{2}+f_{0}^{2}+\lambda^{2}+C_{0}^{2}(k_{x}^{2}+k_{y}^{2})\right]\omega^{2}$
$\displaystyle-[C_{0}^{2}k_{x}\beta+2i\lambda k_{x}^{2}v_{A}^{2}+\lambda
iC_{0}^{2}(k_{x}^{2}+k_{y}^{2})]\omega$
$\displaystyle+k_{x}^{2}v_{A}^{2}\left[k_{x}^{2}v_{A}^{2}+C_{0}^{2}(k_{x}^{2}+k_{y}^{2})\right]=0.$
(38)

The attenuation of the waves is manifested in the imaginary part of the
frequency $\omega$, which translates into an exponential decay of the Fourier
modes. Although Equation (38) appears quite delicate, the damping behavior
turns out to be surprisingly simple for most wave modes. In the limit of high-
frequency magneto-inertia waves, Equation (38) approximates to

$\displaystyle\omega\approx\pm\sqrt{f_{0}^{2}+C_{0}^{2}(k_{x}^{2}+k_{y}^{2})+2k_{x}^{2}v_{A}^{2}}-i\lambda.$
(39)

The real part of Equation (39) is the same as in the inviscid case (18), here
Rayleigh frictions does not do anything except letting the waves decay
exponentially, where the decay rate directly corresponds to the damping rate
$\lambda$. For the slow frequency solutions in Equation (38) in the limit of
large-scale waves $k_{x},k_{y}\rightarrow 0$ we find

$\displaystyle\omega^{2}+\frac{R_{R}^{2}\beta k_{x}+i\lambda
R_{R}^{2}(k_{x}^{2}+k_{y}^{2})}{1+\lambda^{2}/f_{0}^{2}+R_{R}^{2}(k_{x}^{2}+k_{y}^{2})}\omega$
$\displaystyle-\frac{R_{R}^{2}v_{A}^{2}k_{x}^{2}(k_{x}^{2}+k_{y}^{2})}{1+\lambda^{2}/f_{0}^{2}+R_{R}^{2}(k_{x}^{2}+k_{y}^{2})}=0.$
(40)

In contrast to the inertia wave solution, damping is here also affecting the
real part of $\omega$ by the $\lambda^{2}/f_{0}^{2}$ term in the denominator,
effectively reducing the natural frequencies. However, this detuning effect is
negligibly small if we focus on underdamped waves($\lambda<\omega$), where we
readily find $\lambda/f_{0}\ll 1$. If we further, as in section 3.1.1,
consider the radiative part of the tachocline, Equation (40) simplifies to

$\displaystyle\omega^{2}+\frac{\beta
k_{x}}{k_{x}^{2}+k_{y}^{2}}\omega+i\lambda\omega-v_{A}^{2}k_{x}^{2}=0$ (41)

as the analogue to Equation (20). The hydrodynamic Rossby branch

$\displaystyle\omega\approx-\frac{k_{x}\beta}{k_{x}^{2}+k_{y}^{2}}-i\lambda$
(42)

is in the same way affected by damping as magneto-inertia waves (39), the
damping rate $\lambda$ is equal to the exponential decay rate. Intriguingly,
we find a more complex and novel behavior for magneto-Rossby waves, underlying
the dispersion relation

$\displaystyle\omega\approx\frac{k_{x}^{2}v_{A}^{2}(k_{x}^{2}+k_{y}^{2})}{k_{x}\beta+i\lambda(k_{x}^{2}+k_{y}^{2})}$
$\displaystyle=\frac{v_{A}^{2}\beta
k_{x}^{3}(k_{x}^{2}+k_{y}^{2})}{\beta^{2}k_{x}^{2}+\lambda^{2}(k_{x}^{2}+k_{y}^{2})^{2}}-\frac{i\lambda
v_{A}^{2}k_{x}^{2}(k_{x}^{2}+k_{y}^{2})^{2}}{\beta^{2}k_{x}^{2}+\lambda^{2}(k_{x}^{2}+k_{y}^{2})^{2}}.$
(43)

Since this solution is only valid for large scales $k_{x}R_{0}\lesssim 1$, see
Zaqarashvili et al. (2007), we can again neglect the detuning term
$\lambda^{2}(k_{x}^{2}+k_{y}^{2})^{2}$, allowing to write the dispersion
relation as follow

$\displaystyle\omega=\omega_{0}-i\lambda\frac{\omega_{0}^{2}}{\omega_{A}^{2}},$
(44)

where $\omega_{0}$ is the inviscid natural frequency of slow magneto-Rossby
waves (22) and $\omega_{A}=\pm k_{x}v_{A}$ the Alfvén frequency. This time,
the decay rate differs from the damping rate by the factor
$\omega_{0}^{2}/\omega_{A}^{2}$. Zaqarashvili et al. (2007) have shown that
the magneto-Rossby waves are always slower than Alfvén waves (in contrast to
hydrodynamic Rossby waves (42)) such that we always have
$\omega_{0}^{2}/\omega_{A}^{2}<1$ effectively reducing the decay rate. As a
consequence, we can conclude that magneto-Rossby waves are more resilient to
(linear) viscous damping than any other planetary wave modes, giving some more
evidence that magneto-Rossby waves are indeed more relevant for the space
weather dynamics than classic Rossby waves (Dikpati & McIntosh, 2020).

As a last remark, we want to stress the practical benefit of our approach to
include damping. Treating Rayleigh’s linear friction law as the only
dissipation source is surely an oversimplification, but it can also be viewed
as an empirical law since we know how it translates to decay rates in the
different solutions. In principle, one can experimentally measure the
exponential decay rate and translate it back to the damping rate $\lambda$
(e.g., by applying the factor $\omega_{0}^{2}/\omega_{A}^{2}$ in the case of
magneto-Rossby waves). This way, the model can correctly reflect all
participating linear damping mechanisms and not only viscous damping. Such
empirical approaches have been successfully implemented for simple free-
surface wave systems (Horstmann et al., 2020, 2021). In these controllable
mechanical systems it is, however, easy to determine decay rates, while it
seems rather hopeless to extract them from space weather data, especially
since solar Rossby waves are constantly and unpredictably excited.
Nevertheless, one can at least roughly estimate the order of magnitude of
$\lambda$, e.g., by measuring typical Rossby wave amplitudes (Mandal &
Hanasoge (2020) determined amplitudes in the order of $1\,{\rm m}\,{\rm
s}^{-1}$) and estimating the energy input, i.e., the energy stored in the
differential rotation and toroidal magnetic field (Dikpati & Gilman, 2001).
Apart from that, it is worthwhile for future studies to investigate further
dissipation mechanisms, be it Newtonian cooling Chang & Lim (1982); Wu et al.
(2001); Tsai et al. (2014) or magnetic damping, which have not yet been
regarded in the framework of magnetohydrodynamic Rossby waves.

####  3.2.2 Equatorial waves

At the equator, damped waves are more complex and different to non-equatorial
waves since friction does not only alter the frequencies but also the
meridional structure and phase of the wave profiles, as we will see in the
following. In the vicinity of the equator, we have $f\approx\beta y$ and
Equation (12) transforms into

Figure 3: Dispersion curves of different wave modes calculated from equation
(48) for the same parameters used by Zaqarashvili (2018), i.e.,
$\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$, $R_{0}=5\cdot 10^{10}\,{\rm cm}$,
$C_{0}=13\cdot 10^{3}\,{\rm cm}\,{\rm s}^{-1}$, $v_{A}=12.6\cdot 10^{3}\,{\rm
cm}\,{\rm s}^{-1}$, $\beta=1.04\cdot 10^{-16}\,{\rm cm}^{-1}\,{\rm s}^{-1}$
under the action of the damping rate $\lambda=0.01\Omega_{0}$. The blue curves
correspond to magneto-inertia waves, the green curve to a mixed magneto-
Rossby-gravity wave solution and the black curve describes the magneto and
hydrodynamic Rossby waves with $n=1$.

$\displaystyle\square^{2}_{v_{A}}v-C_{0}^{2}\square_{v_{A}}\Delta
v+(\beta^{2}y^{2}+\lambda^{2})\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\beta\frac{\partial}{\partial x}\frac{\partial v}{\partial t}$
$\displaystyle+2\lambda\frac{\partial}{\partial t}\square_{v_{A}}v-\lambda
C_{0}^{2}\Delta\frac{\partial v}{\partial t}=0$ (45)

for $V=0$. By using the same ansatz as before, $v=v_{y}(y)\exp(ik_{x}x-i\omega
t)$, one readily gets the determining equation

$\displaystyle\frac{{\rm d}^{2}v_{y}}{{\rm
d}y^{2}}+\left[\frac{\omega^{2}+i\lambda\omega-
k_{x}^{2}v_{A}^{2}}{C_{0}^{2}}-k_{x}^{2}\right.-$
$\displaystyle\left.-\frac{k_{x}\beta\omega}{\omega^{2}+i\lambda\omega-
k_{x}^{2}v_{A}^{2}}-\mu^{2}y^{2}\right]v_{y}=0,$ (46)

where

$\displaystyle\mu=\frac{\beta\omega}{C_{0}\sqrt{\omega^{2}+i\lambda\omega-
k_{x}^{2}v_{A}^{2}}}.$ (47)

The dispersion relation follows again from the solvability condition (28):

$\displaystyle(\omega^{2}+i\lambda\omega-
k_{x}^{2}v_{A}^{2})(\omega^{2}+i\lambda\omega-
k_{x}^{2}(C_{0}^{2}+v_{A}^{2}))-k_{x}\beta C_{0}^{2}\omega$
$\displaystyle=\beta C_{0}|\omega|\sqrt{\omega^{2}+i\lambda\omega-
k_{x}^{2}v_{A}^{2}}(2n+1).$ (48)

The complex-valued square root on the right-hand side highly complicates
further analysis. It is no longer feasible to derive analytical approximations
for the different wave types. Therefore, we analyze the dispersion relation
numerically, which is a bit delicate since the square root introduced spurious
solutions. In order to get around of this issue we squared both sides of the
dispersion relation—for the price of introducing four invalid solutions—and
determined the eight zeros of the resulting polynomial by calculating the
eigenvalues of the companion matrix. Afterwards, we checked the validity of
the solutions to guarantee that they fulfill the initial dispersion relation
(48) . Figure 3 shows the wave frequencies Re$(\omega)$ as a function of the
wave numbers $k_{x}$ on the example of the solar tachocline with the same
parameters as used by Zaqarashvili (2018) under the effect of a hypothetical
friction parameter of $\lambda=0.01\Omega_{0}$.

Figure 4: Decay rates ${\rm Im}(\omega)$ corresponding to the frequencies
shown in figure 3. The blue curves correspond to magneto-inertia waves, the
green curve to a mixed magneto-Rossby-gravity wave solution and the black
curve describes the magneto and hydrodynamic Rossby waves with $n=1$.

When compared to the undamped solutions, no difference can be discerned for
all inertia, gravity and Rossby wave solutions. We found that the effect of
$\lambda$ on the eigenfrequencies is entirely negligible except for large
friction coefficients in the order of $\lambda\sim|\omega|$. The corresponding
decay rates, however, show a more complex behavior than the non-equatorial
counterparts. Figure 4 shows the imaginary part of the frequency Im$(\omega)$
of the same solutions. It can be recognized that the resulting damping rates
are smaller than for non-equatorial waves under the same conditions. The decay
behavior of magneto-inertia waves shows an interesting length-scale
dependence. Large-scale Poincaré waves ($k_{x}\rightarrow 0$) diminish with
Im$(\omega)\approx 0.75\lambda$, whereas the small-scale representatives
approach the smaller decay rate of Im$(\omega)=0.5\lambda$ for
$k_{x}\rightarrow\pm\infty$. Furthermore, it can be observed that higher wave
modes tend to dissipate more rapidly. Rossby waves always decay with
Im$(\omega)=0.5\lambda$ independently of the wavelength, which is half as fast
as non-equatorial Rossby waves, see Equation (42). This can be explained by
the fact that Rossby waves are bordered by Alfvén waves and have very similar
eigenfrequencies (Zaqarashvili, 2018). The frequency of Alfvén waves, however,
is governed by the square root term of equation (48)
$\omega^{2}+i\lambda\omega-k_{x}^{2}v_{A}^{2}=0$, yielding a decay rate of
Im$(\omega)=-\lambda/2$ in agreement to the magneto-Rossby waves.

Figure 5: Meridional distributions of $n=1$ Rossby waves in the equatorial
waveguide at $\omega=-1\cdot 10^{-7}\,{\rm s}\equiv k_{x}R=2$ for different
Alfvén velocities $v_{A}$. The latitude coordinate $y$ is normalized by the
inviscid and non-magnetic Rossby scale $\sqrt{C_{0}/\beta}$. The dashed lines
correspond to inviscid waves and the solid lines to damped waves with a
friction coefficient of $\lambda=|\omega|$.

So far we have analyzed damped waves which decrease in time. In many cases, as
for some forced wave systems, it is also interesting to study free wave
motions under the effect of damping. Free waves always have real frequencies,
such that the dispersion relation (48) then necessarily introduces complex
wave numbers $k_{x}$ for $\lambda>0$, letting the waves decay spatially in the
zonal direction. Such waves can evolve under spatially limited but constant
forcing, whereby the wave amplitudes are maintained in the region of forcing
and spatially fade away outside. In geophysics it is well known that friction
can modify the meridional scales and introduce a phase shift in the meridional
structure of planetary waves (Mofjeld, 1981; Yamagata & Philander, 1985).
Therefore, the question arises to what extent meridional scales of magneto-
Rossby waves can be affected. For real frequencies, the wave number of slow
and large-scale magneto-Rossby waves can be approximated as follows:

$\displaystyle
k_{x}^{2}=\frac{(2n+1)^{2}(\omega^{2}+i\lambda\omega)}{C_{0}^{2}+v_{A}^{2}(2n+1)^{2}}.$
(49)

We use this dispersion relation for the following calculations. Figures 5 and
6 show normalized meridional velocity profiles $v(y)$ of $n=1$ and $n=2$
magneto-Rossby waves for different Alfvén velocities $v_{A}$ without (dashed
lines) and with strong damping $\lambda=\omega$ (solid lines) in accordance to
Mofjeld (1981). Two opposing tendencies can be observed. As with classic
Rossby waves, friction always broadens the wave profiles. However, this effect
is increasingly counterbalanced by the presence of toroidal magnetic fields
that tight the magneto-Rossby wave profiles. This was to be expected since the
oscillatory scale (29) more and more decreases as magneto-Rossby frequencies
approach Alfvén wave frequencies for $v_{A}\rightarrow\infty$. Altogether, it
can be concluded that the meridional scale variation due to friction is
negligible relative to the magnetic contribution. Finally, we can also see
that friction causes larger absolute latitudinal velocities (but smaller zonal
velocities, not shown here). This is in conformity with the observations by
Mofjeld (1981), there is no significant difference between classic and
magnetic Rossby waves.

Figure 6: Meridional distributions of $n=2$ Rossby waves in the equatorial
waveguide at $\omega=-1\cdot 10^{-7}\,{\rm s}\equiv k_{x}R=2$ for different
Alfvén velocities $v_{A}$. The latitude coordinate $y$ is normalized by the
inviscid and non-magnetic Rossby scale $\sqrt{C_{0}/\beta}$. The dashed lines
correspond to inviscid waves and the solid lines to damped waves with a
friction coefficient of $\lambda=|\omega|$.

###  3.3 Forced wave dynamics

Building on these preliminary investigations, we can now turn on to the study
of forced wave solutions. For the sake of keeping the analysis as simple, we
focus on wave responses to a single tide generating planet of mass $M_{t}$
moving at a fixed distance $r$ with the angular frequency $\Omega_{t}$ around
a rotating star within the equatorial plane (zero declination). We show in
appendix B that the resulting leading-order tidal potential acting in the
$\beta$-plane approximates to

$\displaystyle
V=K\left(\frac{1}{2}+\frac{y}{R_{0}}\right)\left[1+\cos\left(\frac{2x}{R_{0}}-2(\Omega_{t}-\Omega_{0})t\right)\right]$
(50)

at mid-latitudes $\phi_{0}=45^{\circ}$ and to

$\displaystyle
V=K\left[1+\cos\left(\frac{2x}{R_{0}}-2(\Omega_{t}-\Omega_{0})t\right)\right]$
(51)

in the vicinity of the equator $\phi_{0}=0^{\circ}$. The forcing amplitude is
given as

$\displaystyle
K=\frac{3}{4}\frac{GM_{t}}{R_{0}}\left(\frac{R_{0}}{r}\right)^{3}$ (52)

with $G$ referring to the gravitational constant. Two characteristic forcing
frequencies appear in the potentials; the rotation frequency of the star
$\Omega_{0}$ and an external frequency $\Omega_{t}$ dictated by the tide
generating planet. If we consider the Sun exposed to the tidally dominant
planet Jupiter, we find very different periods of $T_{0}\approx$ 25 days and
of $T_{t}\approx 11$ years, so that, in principle, a high bandwidth of
planetary waves can be excited ranging from high frequency magneto-inertia to
very slow magneto-Rossby waves. In the following sections, we present explicit
solution of wave responses under resonant and non-resonant forcing conditions
for both unbounded and equatorially trapped waves. For didactic considerations
we start this time with equatorial waves.

####  3.3.1 Equatorial waves

Inserting the equatorial potential (51) into the wave equation (12) yields

$\displaystyle\square^{2}_{v_{A}}v-C_{0}^{2}\square_{v_{A}}\Delta
v+\beta^{2}y^{2}\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\beta\frac{\partial}{\partial x}\frac{\partial v}{\partial t}$
$\displaystyle+2\lambda\frac{\partial}{\partial t}\square_{v_{A}}v-\lambda
C_{0}^{2}\Delta\frac{\partial v}{\partial
t}+\lambda^{2}\frac{\partial^{2}v}{\partial t^{2}}=\beta
y\frac{\partial}{\partial x}\frac{\partial^{2}V}{\partial t^{2}}$
$\displaystyle=\beta
y\frac{8K\Omega^{2}}{R_{0}}\sin\left(\frac{2x}{R_{0}}-2\Omega t\right).$ (53)

Figure 7: Normalized velocity amplitudes due to solution (57) for the
hydrodynamic $v_{A}=0\,{\rm cm}\,{\rm s}^{-1}$ (top) and magnetic case
$v_{A}=12.6\cdot 10^{3}\,{\rm cm}\,{\rm s}^{-1}$ (bottom) as a function of the
normalized forcing frequency $\Omega$. The calculated frequency range includes
the first three retrograde propagating magneto-Rossby and magneto-inertia wave
modes. We used $\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$, $R_{0}=5\cdot
10^{10}\,{\rm cm}$, $C_{0}=13\cdot 10^{3}\,{\rm cm}\,{\rm s}^{-1}$,
$\beta=1.04\cdot 10^{-16}\,{\rm cm}^{-1}\,{\rm s}^{-1}$.

The forcing frequencies were unified to $\Omega_{t}-\Omega_{0}=\Omega$ for the
sake of clarity. Since the equatorial potential is independent of the latitude
$y$ only the Coriolis forcing term $\beta y\partial_{x}\partial^{2}_{t}V$
remains such that the leading order tidal action is affected neither by the
Alfvén speed $v_{A}$ nor by damping $\lambda$, and therefore manifests itself
in the same manner as with classic Rossby waves. The Coriolis parameter
$f\approx\beta y$ reintroduces the $y$ coordinate, which requires us to expand
$y$ as an orthogonal series of parabolic functions:

$\displaystyle
y=\sum_{n=1}^{\infty}\frac{2^{\frac{5}{2}-2n}}{\sqrt{\mu}(n-1)!}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y).$ (54)

We have used the practical parabolic orthogonality condition

$\displaystyle\int_{-\infty}^{\infty}\exp(-\mu
y^{2})H_{n}(\sqrt{\mu}y)H_{m}(\sqrt{\mu}y){\rm
d}y=2^{n}n!\sqrt{\frac{\pi}{\mu}}\delta_{nm}$ (55)

for the expansion. As a remarkable intermediate result, Equation (54) implies
that only uneven wave modes $n\rightarrow 2n-1$ can respond to tidal forcing
in the leading order, which are always antisymmetric around the equator. We
show in appendix C that the forced wave problem can be solved explicitly by
introducing the ansatz

$v=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\alpha_{m,n}(t)\exp\left(i\frac{mx}{R_{0}}\right)\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{n}(\sqrt{\mu}y),$ (56)

where $\alpha_{m,n}(t)$ are modal coefficients to be determined. The resulting
solution can be expressed as an infinite series of parabolic functions:

$\displaystyle
v(x,y,t)=\sum_{n=1}^{\infty}\Lambda_{n}\exp\left(i\frac{2x}{R_{0}}-i2\Omega
t\right)\exp\left(-\frac{\mu y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y),$ (57)
$\displaystyle{\rm with}\ \ \
\mu=\frac{\beta|2\Omega|}{C_{0}\sqrt{4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}}}\
\ \ {\rm and}$
$\displaystyle\Lambda_{n}=\frac{8iK\beta\Omega^{2}2^{\frac{5}{2}-2n}(R\sqrt{\mu}(n-1)!)^{-1}}{\left(4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}\right)\left(4\Omega^{2}+2i\lambda\Omega-\frac{4}{R_{0}^{2}}(C_{0}^{2}+v_{A}^{2})\right)-C_{0}^{2}\frac{4\beta}{R_{0}}\Omega-
C_{0}\beta|2\Omega|\sqrt{4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}}(4n-1)}.$

Solution (57) consists of two parts. A real part $\sim\text{Re}[\Lambda_{n}]$
describing non-resonant forced waves between the different wave modes and an
imaginary part $\sim\text{Im}[\Lambda_{n}]$ capturing resonant waves in the
vicinity of the eigenfrequencies described by the dispersion relation (38) for
$\omega=2\Omega$. The velocity amplitude scales with $\sim K\sqrt{\beta
C_{0}}$, showing that pure Alfvén waves ($\beta=C_{0}=0$) do not respond to
tidal forcing at the equator. Besides, we would like to draw attention to the
fact that (57) is formally only a specific solution of Equation (53), the
general solution is composed of the specific and homogeneous solution. The
homogeneous part, however, only describes transient phases, e.g., in classical
sloshing experiments the initial wave motions taking place directly after
switching on the shaking table, and decays exponentially after some settling
time. After the transient phase the saturated, quasi-steady wave responses are
fully governed by the specific solution presented here.

For deeper insights into the dynamics of forced planetary waves, we have
calculated the wave responses as a function of the forcing frequency $\Omega$
within a range comprising the first three (magneto)-Rossby and
(magneto)-inertia wave modes $n=1,2,3\rightarrow H_{1},H_{3},H_{5}$. The upper
plot in Figure 7 shows normalized wave amplitudes of classical planetary waves
$v_{A}=0$ for different dynamic damping coefficients $\lambda$ related to the
forcing frequency $\Omega$, so that $\lambda$ always references to the
characteristic time in which equivalent unforced waves would decay. For
example, the $\lambda=0.2\Omega$ corresponds to a damping rate of $0.1\Omega$,
see the Rossby wave solution in Figure 4, which would let the wave diminish
substantially within around ten periods following an artificial elimination of
the excitation force. Figure 7 reveals that the first Rossby mode $n=1$ is
excited most strongly, about five times stronger than the first inertia-
gravity wave and up to 25 times stronger (for $\lambda=0.1\Omega$) than non-
resonant waves occurring in the frequency band between the low-frequency
Rossby and high-frequency inertia wave solutions. Further it can be recognized
that the peak amplitudes decrease rapidly along with increasing wave numbers
showing that the large-scale modes are always most significant for tidal
interactions. As another important general result, it is evident from the
resonance curves that different planetary waves will respond to different
classes of planet-hosting stars. Is the star’s rotation far higher than the
planet’s orbit frequency $|\Omega_{0}|\gg|\Omega_{t}|$, as it is the case for
our Sun forced by Jupiter, only fast inertia-gravity waves will be exited. Are
both angular frequency very close to each other (but not equal), as often the
case for stars hosting tidally-locked planets or also for binary stars, slow
Rossby modes are expected to be stimulated instead. This result is only true
for stars hosting one single (significant) tide-generating planet. If there
are several planets involved, Rossby waves can still be excited by low-
frequency alignment periodicities, e.g., spring tides, visible in the envelope
of the combined tidal potentials, even if we have
$|\Omega_{0}|\gg|\Omega_{t}|$ for all participating planets, see section 4.3.

For comparison, the lower plot in Figure 7 shows the same frequency spectrum
but for magnetic planetary waves with $v_{A}=12.6\cdot 10^{3}\,{\rm cm}\,{\rm
s}^{-1}$. In accordance with Zaqarashvili (2018), it can be seen that Rossby-
waves below the Alfvén frequency (vertical black line) are largely suppressed.
Magneto-Rossby waves accumulate closely above the Alfvén frequency, where they
respond with much smaller amplitudes. The inertia waves, in contrast, remain
largely unaffected. As a novel and rather striking result, we can also see
that damping allows the waves to overcome the Alfvén frequency barrier.
Inviscid planetary waves below $\omega=\pm 2v_{A}/R_{0}$ do not exists,
whereas damped waves are still able to respond at frequencies far below
$2|\Omega|/\Omega_{0}=10^{-1}$. This behavior is intriguing, but it must be
noted that constant toroidal magnetic fields generally inhibit Rossby waves
and antagonize the tidal response, making these waves unlikely to be of any
practical importance. At this point, it might be of interest for future
studies to analyze the tidal response of Rossby waves subject to more
realistic nonuniform and also oscillatory magnetic fields, which are not
exposed to a cut-off frequency (Zaqarashvili, 2018).

####  3.3.2 Non-equatorial waves

We proceed with solving the forced wave problem at mid-latitudes, which is a
bit more intricate insofar that here the tidal potential (50) depends, in
contrast to the equatorial potential (51), on the local latitude $y$ so that
all three potential terms in Equation (12) must be taken into account. We
obtain the following forced wave equation:

$\displaystyle\square^{2}_{v_{A}}v-C_{0}^{2}\square_{v_{A}}\Delta
v+f_{0}^{2}\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\beta\frac{\partial}{\partial x}\frac{\partial v}{\partial
t}+2\lambda\frac{\partial}{\partial t}\square_{v_{A}}v-\lambda
C_{0}^{2}\Delta\frac{\partial v}{\partial
t}+\lambda^{2}\frac{\partial^{2}v}{\partial
t^{2}}=f_{0}\frac{\partial}{\partial x}\frac{\partial^{2}V}{\partial
t^{2}}-\lambda\frac{\partial}{\partial y}\frac{\partial^{2}V}{\partial
t^{2}}-\frac{\partial}{\partial t}\frac{\partial}{\partial y}\square_{v_{A}}V$
$\displaystyle=\left[f_{0}\Omega+2\Omega^{2}-\frac{2v_{A}^{2}}{R_{0}^{2}}+\frac{2f_{0}\Omega}{R_{0}}y\right]\frac{4K\Omega}{R_{0}}\sin\left(\frac{2x}{R_{0}}-2\Omega
t\right)+\frac{4K\lambda\Omega^{2}}{R_{0}}\cos\left(\frac{2x}{R_{0}}-2\Omega
t\right)$ (58)

Both constant and $y$-proportional terms remain on the right-hand side of
(58), which moreover involves different phases $\sim\sin(2x/R_{0}-2\Omega t)$
and $\sim\cos(2x/R_{0}-2\Omega t)$. Hence, we are required to expand two
different Fourier series for constant terms and $y$, which, however, gives
rise to the difficulty that mid-latitude waves are meridionally unbounded so
that the wave numbers $k_{y}$ are in principle arbitrary. Hence, we need to
constrain the meridional scale in a purposeful way. On the one hand, we are
interested in large-scale responses, but the $\beta$-plane approximation
becomes more and more inaccurate with increasing meridional dimensions on the
other hand. As the best compromise, we confine the wave problem into a
waveguide defined by the interval $-R_{0}/2\leq y\leq R_{0}/2$, with $y=0$
defining the latitude $\phi_{0}=45^{\circ}$. In spherical coordinates, this
interval corresponds approximately to the meridional band
$15^{\circ}\lesssim\phi_{0}\lesssim 75^{\circ}$, the respective latitudinal
range of $60^{\circ}$ is widely regarded as the non-equatorial $\beta$-plane
limit. Within the interval $-R_{0}/2\leq y\leq R_{0}/2$, we can describe the
forcing terms as Fourier series by expanding

$1=\sum_{n=1}^{\infty}\frac{4}{(2n-1)\pi}(-1)^{n-1}\cos\left(\frac{(2n-1)\pi
y}{R_{0}}\right)$ (59)

and

$\displaystyle
y=\sum_{n=1}^{\infty}\frac{4R_{0}}{(2n-1)^{2}\pi^{2}}(-1)^{n-1}\sin\left(\frac{(2n-1)\pi
y}{R_{0}}\right).$ (60)

Interestingly, both series only involve uneven latitudinal wave modes $2n-1$,
which are the only modes that can respond in the leading order. Hence, the
tidal potential, although manifested very differently nearby to and far away
from the equator, imposes the same symmetry on mid-latitude waves as it does
on equatorial waves. The forced wave problem can now be solved most easily by
inserting the ansatz

$\displaystyle v$
$\displaystyle=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\left[\alpha_{m,n}(t)\sin\left(\frac{m}{R_{0}}x\right)+\beta_{m,n}(t)\cos\left(\frac{m}{R_{0}}x\right)\right]$
$\displaystyle\times\exp\left(\frac{in\pi}{R_{0}}y\right)$ (61)

into Equation (58) as we show in Appendix D. Since the right-hand side forcing
terms in Equation (58) contains terms both symmetric and antisymmetric in the
($x,t$) space, we need two independent modal coefficients $\alpha_{m,n}(t)$
and $\beta_{m,n}(t)$ to solve the problem. We find the following explicit
analytic solution:

$\displaystyle v(x,y,t)$
$\displaystyle=\sum_{n=1}^{\infty}\frac{\bar{K}_{2}+i\bar{K}_{1}}{2\left(16\Omega^{4}-4\Omega^{2}\bar{\Lambda}_{1}+\bar{\Lambda}_{3}-\frac{4\beta
C_{0}^{2}\Omega}{R_{0}}\right)-4i\Omega(8\lambda\Omega^{2}-\bar{\Lambda}_{2})}\exp\left(2i\Omega
t-2i\frac{x}{R_{0}}+\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+\sum_{n=1}^{\infty}\frac{\bar{K}_{2}-i\bar{K}_{1}}{2\left(16\Omega^{4}-4\Omega^{2}\bar{\Lambda}_{1}+\bar{\Lambda}_{3}-\frac{4\beta
C_{0}^{2}\Omega}{R_{0}}\right)+4i\Omega(8\lambda\Omega^{2}-\bar{\Lambda}_{2})}\exp\left(-2i\Omega
t+2i\frac{x}{R_{0}}+\frac{i(2n-1)\pi}{R_{0}}y\right).$ (62)

Figure 8: Normalized velocity amplitudes due to solution (62) for
$v_{A}=0\,{\rm cm}\,{\rm s}^{-1}$ (top) and $v_{A}=12.6\cdot 10^{3}\,{\rm
cm}\,{\rm s}^{-1}$ (bottom) as a function of the normalized forcing frequency
$\Omega$. The calculated frequency range includes the first three retrograde
propagating magneto-Rossby and magneto-inertia wave modes. We used
$\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$, $R_{0}=5\cdot 10^{10}\,{\rm cm}$,
$C_{0}=13\cdot 10^{3}\,{\rm cm}\,{\rm s}^{-1}$, $\beta=1.04\cdot
10^{-16}\,{\rm cm}^{-1}\,{\rm s}^{-1}$.

Figure 9: Normalized velocity amplitudes due to solution (62) for
$v_{A}=0\,{\rm cm}\,{\rm s}^{-1}$ (top) and $v_{A}=12.6\cdot 10^{3}\,{\rm
cm}\,{\rm s}^{-1}$ (bottom) as a function of the normalized forcing frequency
$\Omega$. The calculated frequency range includes the first three prograde
propagating magneto-Rossby and magneto-inertia wave modes. We used
$\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$, $R_{0}=5\cdot 10^{10}\,{\rm cm}$,
$C_{0}=13\cdot 10^{3}\,{\rm cm}\,{\rm s}^{-1}$, $\beta=1.04\cdot
10^{-16}\,{\rm cm}^{-1}\,{\rm s}^{-1}$.

The introduced coefficients
$\bar{\Lambda}_{1},\bar{\Lambda}_{2},\bar{\Lambda}_{3},\bar{K}_{1},\bar{K}_{2}$
are specified in Appendix D11. In the solution, the symmetric part
$\sim\cos(...)$ describes non-resonant waves between the resonances, and the
antisymmetric part $\sim\sin(...)$ captures resonant waves in the vicinity of
the eigenfrequencies, respectively. As a remarkable difference to the
equatorial solution (57), the responses (62) do not fade out in the limit of
zero Coriolis forces $f_{0}=\beta=0$ and zero gravity $C_{0}=0$, showing that
Alfvén waves can be directly excited by tidal forces. In the limit of pure
Alfvén waves $\Lambda_{n}$ simplifies to,

$\displaystyle\Lambda_{n}=\frac{8i(-1)^{n-1}K\Omega}{R_{0}(2n-1)\pi\left[4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}\right]},$

which yields, by reintroducing the series (59), the rather simple,
unidirectional wave solution

$\displaystyle
v=\frac{2K\Omega}{R_{0}\left[4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}\right]}\exp\left(i\frac{2x}{R_{0}}-i2\Omega
t\right).$ (63)

This solution is independent of the tachocline layer thickness $H_{0}$ and
does further not rely on the presence of any reduced gravity $g$ such that
Equation (63) is not confined to our shallow water model and can be used in a
more universal way to estimate tidally excited Alfvén waves in stars,
regardless of the exact strata properties. At resonances Alfvén waves have
maximum amplitudes of $v=K/\lambda$ so that further modeling efforts to
estimate the magnetic and viscous dissipation of Alfvén waves
$\lambda_{\text{Alfv\'{e}n}}$ (that is physically different to planetary
waves) seems necessary for future studies.

Figure 10: Resonant velocity amplitudes of the first $n=1$ retrograde wave
mode for Rossby (left) and inertia-gravity (right) waves as a function of
$C_{0}$ and $\lambda$. The parameters $v_{A}=0\,{\rm cm}\,{\rm s}^{-1}$,
$\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$, $R_{0}=5\cdot 10^{10}\,{\rm cm}$,
$\beta=1.04\cdot 10^{-16}\,{\rm cm}^{-1}\,{\rm s}^{-1}$ and $K=504\,{\rm
cm}^{2}\,{\rm s}^{-2}$ were used for the calculation.

Similarly as in section 3.3.1, we study the wave responses graphically through
calculating normalized amplitudes as a function of the tidal frequency
$\Omega$ for the first three wave modes $n=1,2,3$ and for different damping
ratios. Since, in contrast to equatorial waves, the retrograde solution branch
comprising classic Rossby waves and the prograde branch describing magneto-
Rossby waves differ significantly, we study retrograde and prograde waves
separately. Retrograde wave responses are shown in Figure 8 without (top) and
with (bottom) the presence of a toroidal magnetic field. As the most striking
difference to equatorial waves, it appears that the first $n=1$ Rossby mode
resonates with an orders of magnitude higher amplitude than any other wave
mode, which has made it necessary to present the velocity amplitudes on a
logarithmic scale. Higher mode responses and non-resonant waves seem to be
insignificant at mid latitudes. When a magnetic field is applied, the Rossby
eigenfrequencies are shifted towards higher frequencies and, rather
interestingly, we find an inverted dispersion, meaning that the
eigenfrequencies increase with growing wave numbers instead of being reduced.
In Figure 9 we also show the prograde wave responses in the hydrodynamic case
(top) and under the effect of the toroidal magnetic field (bottom). Prograde
Rossby waves do not exist in the hydrodynamic limit, which is why only high-
frequency responses of inertia-gravity waves remain visible. In the magnetic
case, we find almost the same response pattern as for retrograde waves, with
the small difference that the eigenfrequencies are slightly smaller. We can
conclude that both progradely and retrogradely excited Rossby waves may
resonate to a similar extent—an intriguing difference to analogous geophysical
wave responses, where the tidal sense always plays a pivotal role.

##  4 Estimation of wave responses in our sun

In the previous sections we have analyzed characteristic tidal wave
excitations for the general class of solar-like stars. The amplitudes were
normalized to keep the analysis as general as possible, but we had to use
fixed values for the parameters $C_{0}$ and $v_{A}$ to calculate the response
patterns, which both affect the eigenfrequencies and thereby shift the
resonance peaks. Apart from this effect, the presented wave responses are
generally valid and we can easily deduce, e.g., that the first latitudinal
wave mode $n=1$ is invariably showing the most significant response. That was
certainly to be expected due the large-scale nature of the tidal potential,
but the difference to the second mode $n=2$ is substantial. In order to detect
planetary waves in solar-like stars which are mainly excited by some tidally
dominant planet, e.g., (hot) Jupiters, our analysis suggests to search for
$m,n=2,1$ responses, for which the tidal energy input is by far the highest in
the leading order. Starting from this insight, we want to estimate the maximum
attainable velocity amplitudes in our Sun, which can result from tidal
forcing. The tidal forces experienced by the Sun are mainly dictated by a
complex interplay of Jupiter, Venus, Earth (and Mercury), which involves many
different excitation frequencies, among them the solar rotation, the
individual orbit frequencies and, perhaps most intriguing, also tidal
variations with periods of 11 years closely related to the solar cycle
(Okhlopkov, 2016). At time-scales close to the solar rotation $\Omega_{0}$,
waves are always excited retrogradely, whereas the tide-generating planets
force planetary waves progradely at time scales in the order of their orbit
periods. Here, we simplify these complex dynamics by considering Jupiter as
the sole tide-generating body (exactly the case for which we derived the tidal
potentials (50) and (51)) and taking the forcing frequency $\Omega$ as
arbitrary in the first instance. The latter idealization is justified in view
of the argument that the buoyancy frequency (governed by the effective
gravity) increases from zero to several rotations per minute as we descend in
the tachocline from the convection zone down to the radiative interior. This
means there is always a place in the tachocline where the buoyancy period can
match the tidal period, or, to put it another way, the tachocline can in
principal resonate to any given excitation frequency. In the following we
calculate resonant amplitudes as a function of the unspecified quantities
$C_{0}$ and $\lambda$, for which we first calculate the eigenfrequencies
$\omega_{n=1}(C_{0},\lambda)$ using Equations (48) and (38) and inserting them
into the solutions (57) and (62).

Figure 11: Resonant velocity amplitudes of the first $n=1$ retrograde wave
mode for Rossby (left) and inertia-gravity (right) waves as a function of
$C_{0}$ and $\lambda$. The parameters $v_{A}=0\,{\rm cm}\,{\rm s}^{-1}$,
$\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$, $R_{0}=5\cdot 10^{10}\,{\rm cm}$,
$\beta=1.04\cdot 10^{-16}\,{\rm cm}^{-1}\,{\rm s}^{-1}$ and $K=504\,{\rm
cm}^{2}\,{\rm s}^{-2}$ were used for the calculation.

###  4.1 Equatorial waves

We start to analyze resonant equatorial wave excitations. At first it needs to
be noted that non-resonant wave responses, as they occur in the large
frequency range between Rossby and gravity inertia wave, see Figure 7, are
about $v\approx 0.02\,{\rm cm}\,{\rm s}^{-1}$. Such low values are entirely
negligible for the solar dynamics, so that we will indeed have to focus on
resonant excitations. Resonant velocity amplitudes for the full range of
gravity velocities considered in the literature and damping ratios between
$0.01$ and $1$ are shown in Figure 10 for both hydrodynamic Rossby (left) and
gravity-inertia (right) waves. We find a large range of possible velocity
amplitudes $10^{-3}\,{\rm cm}\,{\rm s}^{-1}\lesssim v\lesssim 10\,{\rm
cm}\,{\rm s}^{-1}$, where in particular low-frequency Rossby waves in the
regime $10^{3}\,{\rm cm}\,{\rm s}^{-1}\lesssim C_{0}\lesssim 10^{4}\,{\rm
cm}\,{\rm s}^{-1}$, which correspond to periods being in the order of the
Schwabe cycle, can reach amplitude velocities above $v=1\,{\rm cm}\,{\rm
s}^{-1}$. Such velocities may already be dynamo effective, however, all in
all, tidal excitations of equatorial waves seem to play only a minor role
compared to other mechanism capable of stimulating planetary waves. The
available helioseismic data allows the extraction of $m=2$ Rossby waves only
for RMS velocities larger $v>50\,{\rm cm}\,{\rm s}^{-1}$ (Liang et al., 2019),
such that tidally forced waves at the equator are unverifiable. Up to now,
only sectional modes $m=n$ in the range $3\leq m\leq 15$ have been detected
(Löptien et al., 2018).

Prograde magnetic Rossby waves are largely suppressed by constant toroidal
fields resulting in even weaker wave response not deserving any further
discussion. At this point it seems very promising for future studies to
incorporate nonuniform latitudinal magnetic field profiles into the model as
done by Zaqarashvili (2018), allowing to study very slow prograde magneto-
Rossby waves, which can reach Schwabe-cycle periods for many different
combinations of gravity and Alfvén velocities.

Figure 12: Resonant velocity amplitudes of the first $n=1$ prograde wave mode
for magneto Rossby (left) and magneto inertia-gravity (right) waves as a
function of $C_{0}$ and $\lambda$. The parameters $v_{A}=12.6\cdot
10^{3}\,{\rm cm}\,{\rm s}^{-1}$, $\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$,
$R_{0}=5\cdot 10^{10}\,{\rm cm}$, $\beta=1.04\cdot 10^{-16}\,{\rm
cm}^{-1}\,{\rm s}^{-1}$ and $K=504\,{\rm cm}^{2}\,{\rm s}^{-2}$ were used for
the calculation.

###  4.2 Non-equatorial waves

At mid latitudes, we find still smaller non-resonant velocities about
$v\approx 0.004\,{\rm cm}\,{\rm s}^{-1}$, which about five times smaller than
the corresponding non-resonant equatorial velocities. However, this finding
does not allow us to draw any conclusions about general response tendencies.
As already discussed in section 3.3.2, $n=1$ Rossby waves can reach
extraordinary high amplitudes at resonance despite the low non-resonant
amplitude level. Indeed, Figure 11 confirms that Rossby waves can
(theoretically) resonate with amplitudes of more than $v=1\,{\rm m}\,{\rm
s}^{-1}$, whereas gravity-inertia waves are even less excited than their
equatorial counterparts. In this order of magnitude tidally forced Rossby
waves would indeed be capable of providing sufficient energy to considerably
affect the solar dynamics and then to take part in the synchronization process
of the dynamo, especially since the tidal energy input of the other essential
planets Venus and Earth have been disregarded here.

Figure 13: Resonant velocity amplitudes of the first $n=1$ prograde wave mode
for magneto Rossby (left) and magneto inertia-gravity (right) waves as a
function of $C_{0}$ and $\lambda$. The parameters $v_{A}=1\cdot 10^{3}\,{\rm
cm}\,{\rm s}^{-1}$, $\Omega_{0}=26\cdot 10^{-7}\,{\rm s}^{-1}$, $R_{0}=5\cdot
10^{10}\,{\rm cm}$, $\beta=1.04\cdot 10^{-16}\,{\rm cm}^{-1}\,{\rm s}^{-1}$
and $K=504\,{\rm cm}^{2}\,{\rm s}^{-2}$ were used for the calculation.

At mid latitudes, also prograde magnetic waves can respond significantly to
tidal excitations. For the strong magnetic field scenarios $B_{0}=20\,{\rm
kG}\equiv v_{A}=12.6\cdot 10^{3}\,{\rm cm}\,{\rm s}^{-1}$ discussed before, we
find nevertheless considerably lower velocity amplitudes up to $v=10\,{\rm
cm}\,{\rm s}^{-1}$, see Figure 12, which is comparable to the weak equatorial
wave responses. However, very slow magneto Rossby waves with periods in the
order of the solar cycle and longer are impacted by sufficiently lower
toroidal magnetic fields of $B_{0}\sim 1{\rm kG}$, as they can be generated by
steady non-reversing dynamos or be manifested as an offshoot of the primordial
field in the radiative interior, which may penetrate into the tachocline
(Zaqarashvili et al., 2015). Therefore, in Figure 13 we also show the resonant
responses of prograde waves subject to a slow Alfvén velocity of $v_{A}=1\cdot
10^{3}\,{\rm cm}\,{\rm s}^{-1}$. In this scenario, the amplitudes are
comparable to their retrograde counterparts and may be dynamo-effective in a
larger range of the $\lambda/\Omega$–$C_{0}$ space.

###  4.3 Plausibility of the responses

We have discussed possible wave responses in the solar tachocline under
different scenarios for a wide range of the unspecified quantities $\lambda$
and $C_{0}$. The anticipated amplitudes vary in each case over several orders
of magnitude, so that up to this point it is not possible to make a valid
statement about the specific tidal energy input; only excitation potentials
have been demonstrated. The presented velocity maps also do not allow drawing
any probability information as to which damping regimes can be reached. We
shall therefore briefly discuss the plausibility of the presented scenarios.
First, it is important to be aware that the velocity amplitudes are always
smaller for small $C_{0}$ primarily due to the fact that small $C_{0}$ are
associated with small eigenfrequencies, meaning that the absolute damping
rates $\lambda$ also become smaller since we have considered constant damping
ratios $\lambda/\Omega$. This poses the question of whether this assumption is
realistic, whether low-frequency and high-frequency waves actually have
comparable lifetimes. The answer depends on the nature of damping. Although
solar fluid dynamics is of strongly turbulent nature, leading-order damping of
large-scale Rossby waves is likely governed by laminar viscous boundary
layers. In turbulent convection, the kinematic viscosity should be understood
as an eddy viscosity (Rüdiger, 1989), whose exact value is controversial in
the literature. Gizon et al. (2020) have applied the value $\nu\approx
250\,{\rm km}^{2}\,{\rm s}^{-1}$ associated with supergranular scales to
account for the shearing boundary layers arising from differential rotation.
For turbulent viscosities of such high orders, boundary layers are expected to
be essentially laminar and we can think of oscillatory Stokes-like boundary
layers forming on top of the rigid radiative interior, in which the horizontal
velocities fall off exponentially to fulfill the no-slip boundary condition.
When damping is dominated by such oscillatory layers, the damping rate scales
as $\lambda\sim\sqrt{\Omega}$, such that low frequency waves (corresponding to
small $C_{0}$) are indeed subject to weaker damping. However, the quantities
$\lambda/\Omega$ and $C_{0}$ are only independent if $\lambda\sim\Omega$,
which leads us to conclude that $\lambda/\Omega$ will be distinctly higher for
small $C_{0}$ than for large gravity velocities. Indeed, Liang et al. (2019)
found, in average, longer lifetimes for high-frequency waves at lower Rossby
modes than for low-frequency waves at higher modes. Nevertheless, the overall
damping behavior remains largely unsettled, all the more so since thermal and
magnetic dissipation may also make significant contributions. Therefore, it
seems expedient for future studies to explicitly calculate the different
Stokes and Ekman boundary layers, similar to the authors Bildsten & Ushomirsky
(2000) for the case of neutron stars, in order to provide better estimates of
viscous damping rates, which, apart from the calculation of resonant wave
responses, could also allow for better modeling of horizontal eigenfunctions
observed at the solar surface (Proxauf et al., 2020).

Finally, we need to take a closer look at the direction of excitation. We
presented response scenarios for both prograde and retrograde planetary waves,
however, when focusing on the tidally dominant planet Jupiter, we only have
two frequencies $\Omega_{0}$ and $\Omega_{t}$ involved. Obviously, the sun’s
angular frequency is orders of magnitude higher than Jupiter’s orbit frequency
$\Omega_{0}\gg\Omega_{t}$, letting us approximate
$\Omega=\Omega_{t}-\Omega_{0}\approx-\Omega_{0}$. Apparently only the
retrograde branches, comprising classic Rossby and gravity-inertial waves,
will be stimulated directly. The most relevant and promising scenario to reach
significantly high forced velocity amplitudes is therefore represented in
Figure 11 in the form of hydrodynamic waves excited at mid latitudes. However,
when considering the tidal forces of Jupiter, Venus and Earth in conjunction,
far lower forcing frequencies appear in envelope curves of the potentials,
most relevantly the 11 year period closely matching to the solar activity
periodicity (Okhlopkov, 2016). Envelope variations of the combined tidal
potentials can also excite planetary waves progradely in the most susceptible
low $C_{0}$ regimes—best requirements for the excitation of slow magneto-
Rossby waves. At this point it seems worthwhile to incorporate the tidal
forcing of all three planets and to calculate multifrequential wave responses
numerically. We plan to accomplish this in a future study. As a second
interesting scenario of low-frequency forcing, we can also think off solar-
like stars forced by tidally-synchronized planets, or likewise binary stars,
where the dominant frequency $\Omega_{0}$ is absent and slow frequencies, both
prograde and retrograde, may remain as a result of orbital changes. For all
these systems it is conceivable that low-frequency tidal forcing may dictate
the activity cycles, be it by synchronizing the dynamo or by directly
implanting large-scale atmospheric motions, provided that the forcing
amplitude $K$ is large enough (resp. the damping rate $\lambda$ small enough)
to induce sufficiently high wave amplitudes.

##  5 Concluding remarks

By equipping the magnetohydrodynamic “shallow water” equations with a tidal
potential term and linear friction, we have constructed a first theoretical
set-up for the study of damped and forced planetary waves in the tachocline
layer of solar-like stars. The governing equations were projected onto two
different Cartesian planes in the vicinity of the equator and at mid
latitudes, allowing us to describe both equatorially trapped and locally
unbounded waves in the most approachable way. As a key result, we have shown
that the linearized system of governing equations can be combined into one
decoupled wave equation (12) for the local latitudinal velocity component,
which markedly simplifies the Fourier analysis for extracting the
characteristic dispersion relations in different wave limits of interest.

We solve this wave equation analytically within different regimes, starting
with the known free wave solutions via damped wave dynamics up to the complete
forced wave problem. The analysis revealed that the damping behavior of
magnetohydrodynamic planetary waves is more intricate than the damping of
classical geophysical waves since the introduced damping parameter can
translate into very different decay rates for the different wave types. Most
interesting here is our finding that the damping rates of retrograde Rossby
waves at mid-latitudes correspond exactly to the decay rates, whereas
progradely propagating magnetic Rossby waves are predicted to decay
considerably slower by a factor of the squared natural eigenfrequency divided
by the squared Alfvén frequency. Further, the damped wave solutions derived
for equatorial waves facilitate to calculate the meridional scales as a
function of the damping rate and the toroidal magnetic field, with the result
that damping always widens and the presence of magnetic fields always, in an
even overcompensating way, narrows the equatorial waveguide.

The forced wave problem is solved analytically for the idealized case of a
single tide-generating body prescribing a perfect circular orbit around the
central star. The solutions can describe both non-resonant and resonant wave
responses, the latter, however, are largely determined by the a priori unknown
damping coefficient. We found that for fixed damping ratios equatorial waves
always respond with higher velocity amplitudes than mid-latitudes waves under
non-resonant conditions, whereas mid-latitude waves have higher peak
velocities in proximity to resonances. Among all types of planetary waves, the
first large-scale Rossby mode, be it the classic retrograde or magnetic
prograde Rossby wave, is found to always resonate with the highest amplitudes
when considering fixed lifetimes. Rossby waves are therefore confirmed to be
indeed a most promising candidate to potentially act as a resonance ground for
low-frequency tidal excitations.

Finally, we applied the solutions to the specific scenario of our Sun tidally
forced by Jupiter for estimating possible velocity responses. We obtained non-
resonant amplitudes of $v\approx 0.02\,{\rm cm}\,{\rm s}^{-1}$ at the equator
and of $v\approx 0.004\,{\rm cm}\,{\rm s}^{-1}$ at mid-latitudes, which are,
for the solar dynamo, completely negligible. Resonant amplitudes strongly
depend on both damping and the effective gravity (the considered region in the
tachocline layer) so that resonant velocities deviate by several orders of
magnitude with respect to these parameters. Consequently, a reliable
prediction of the anticipated responses cannot yet be made without further
ado. However, we found that for low-frequency excitations in the order of the
11 year solar cycle, with particular view on the 11.07-years alignment
periodicity of the tidally dominant planets Venus, Earth, and Jupiter, the
tidal energy input of Jupiter alone can evoke high Rossby wave amplitudes of
$v\gtrsim 100\,{\rm cm}\,{\rm s}^{-1}$ for fairly small damping ratios in the
range $0.01\lesssim\lambda/\Omega\lesssim 0.1$. We can conclude that, despite
the fact that tidal accelerations are very small, significant velocities can
potentially be induced through resonant amplification if dissipation is
sufficiently small.

To draw more definite conclusions, our model must be extended in different
directions. Our pilot analysis has stressed the potential that tidal forcing
may induce significant tachoclinic wave motions, but to irrevocably verify or
disprove this possibility further modeling of the different dissipation
sources, including viscous, turbulent, thermal and magnetic damping, is
required. In addition, it is expedient to solve the problem in spherical
coordinates and to include the full potentials of all significant planets,
which would allow to determine the exact response to alignment periodicities
and spring tides. Regardless of the exact attainable amplitude levels, there
also remains the question of how and in what way Rossby waves may affect the
solar dynamo. The linear Rossby wave solutions do not produce net kinetic
helicity, which is the main ingredient for the $\alpha$ effect. Hence, it
appears instructive to look at the nonlinear evolution of tidally excited
Rossby waves, e.g., tachocline nonlinear oscillations (TNOs) arising from the
energy exchange between Rossby waves and differential rotation (Dikpati et
al., 2018b). Moreover, it might be promising to estimate how wavelike
displacements of the tachocline layer affect the entropy stratification to
check if tide stimulated flux tube instabilities can possibly encroach into
the dynamo process (Ferriz-Mas, 1996; Charbonneau, 2022). Finally, coupling
the planetary wave equations with dynamo models may provide valuable insights
into possible synchronization mechanisms imparted by Rossby waves.

This work has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No 787544). It was also supported in frame of the Helmholtz - RSF
Joint Research Group “Magnetohydrodynamic instabilities: Crucial relevance for
large scale liquid metal batteries and the sun-climate connection”, contract
No HRSF-0044. TVZ was supported by the Austrian Fonds zur Förderung der
Wissenschaftlichen Forschung (FWF) project P30695-N27 and by Shota Rustaveli
National Science Foundation of Georgia (project FR-21-467).

##  Appendix A Derivation of the wave equation

Differentiation of equations (8) and (9) with respect to time and using
equations (10) and (11) yields the coupled differential equations for the
velocities

$\displaystyle\frac{\partial^{2}u}{\partial t^{2}}-f\frac{\partial v}{\partial
t}=v_{A}^{2}\frac{\partial^{2}u}{\partial
x^{2}}+C_{0}^{2}\left(\frac{\partial^{2}u}{\partial
x^{2}}+\frac{\partial^{2}v}{\partial x\partial
y}\right)-\frac{\partial^{2}V}{\partial x\partial t}-\lambda\frac{\partial
u}{\partial t},$ (A1) $\displaystyle\frac{\partial^{2}v}{\partial
t^{2}}+f\frac{\partial u}{\partial t}=v_{A}^{2}\frac{\partial^{2}v}{\partial
x^{2}}+C_{0}^{2}\left(\frac{\partial^{2}u}{\partial x\partial
y}+\frac{\partial^{2}v}{\partial^{2}y}\right)-\frac{\partial^{2}V}{\partial
y\partial t}-\lambda\frac{\partial v}{\partial t},$ (A2)

to be decoupled in the following. For this we can replace the term
$f\partial_{t}u$ in (A2) by equation (8), giving

$\displaystyle\square_{v_{A}}v+f^{2}v-fg\frac{\partial\eta}{\partial
x}+f\frac{B_{0}}{4\pi\rho}\frac{\partial b_{x}}{\partial x}-f\frac{\partial
V}{\partial x}-f\lambda u=C_{0}^{2}\frac{\partial}{\partial
y}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial
y}\right)-\frac{\partial^{2}V}{\partial y\partial t}-\lambda\frac{\partial
v}{\partial t}.$ (A3)

As an intermediate step, we calculate $\partial_{y}(\ref{eq:TD1})$ \-
$\partial_{x}(\ref{eq:TD2})$ to extract the following expression

$\displaystyle\frac{\partial\eta}{\partial
t}=\frac{H_{0}}{f}\frac{\partial\zeta}{\partial
t}+\lambda\frac{H_{0}}{f}\zeta+\frac{H_{0}}{f}\beta
v-\frac{H_{0}}{f}\frac{B_{0}}{4\pi\rho}\frac{\partial}{\partial
x}\left(\frac{\partial b_{y}}{\partial x}-\frac{\partial b_{x}}{\partial
y}\right),$ (A4)

where $\zeta=\partial_{x}v-\partial_{y}u$ denotes the horizontal vorticity. We
can now differentiate (A3) with respect to time and insert (A4) to eliminate
$\eta$

$\displaystyle\frac{\partial}{\partial t}\square_{v_{A}}v+f^{2}\frac{\partial
v}{\partial t}-C_{0}^{2}\frac{\partial^{2}\zeta}{\partial t\partial
x}-C_{0}^{2}\lambda\frac{\partial\zeta}{\partial
x}-C_{0}^{2}\beta\frac{\partial V}{\partial
x}+C_{0}^{2}\frac{B_{0}}{4\pi\rho}\frac{\partial^{2}}{\partial
x^{2}}\left(\frac{\partial b_{y}}{\partial x}-\frac{\partial b_{x}}{\partial
y}\right)+fv_{A}^{2}\frac{\partial^{2}u}{\partial
x^{2}}-f\frac{\partial^{2}V}{\partial t\partial x}-f\lambda\frac{\partial
u}{\partial t}$ $\displaystyle=C_{0}^{2}\frac{\partial}{\partial
t}\frac{\partial}{\partial y}\left(\frac{\partial u}{\partial
x}+\frac{\partial v}{\partial y}\right)-\frac{\partial^{3}V}{\partial
y\partial t^{2}}-\lambda\frac{\partial^{2}v}{\partial t^{2}}.$ (A5)

By noting that

$\displaystyle C_{0}^{2}\frac{\partial^{2}\zeta}{\partial t\partial
x}+C_{0}^{2}\frac{\partial}{\partial t}\frac{\partial}{\partial
y}\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial
y}\right)=C_{0}^{2}\frac{\partial}{\partial t}\Delta v$ (A6)

and taking again the time derivative, we find

$\displaystyle\frac{\partial^{2}}{\partial
t^{2}}\square_{v_{A}}v+f^{2}\frac{\partial^{2}v}{\partial
t^{2}}-C_{0}^{2}\frac{\partial^{2}}{\partial t^{2}}\Delta
v-C_{0}^{2}\lambda\frac{\partial^{2}\zeta}{\partial t\partial
x}-C_{0}^{2}\beta\frac{\partial^{2}V}{\partial t\partial
x}+C_{0}^{2}v_{A}^{2}\left(\frac{\partial^{4}v}{\partial
x^{4}}-\frac{\partial^{4}u}{\partial x^{3}\partial
y}\right)+fv_{A}^{2}\frac{\partial^{3}u}{\partial t\partial x^{2}}$
$\displaystyle-f\frac{\partial^{3}V}{\partial t^{2}\partial
x}-f\lambda\frac{\partial^{2}u}{\partial t^{2}}+\frac{\partial^{4}V}{\partial
y\partial t^{3}}+\lambda\frac{\partial^{3}v}{\partial t^{3}}=0.$ (A7)

Let us collect all remaining $u$ velocity terms in (A7), reading

$\displaystyle C_{0}^{2}\lambda\frac{\partial^{3}u}{\partial x\partial
y\partial t}-C_{0}^{2}v_{A}^{2}\frac{\partial^{4}u}{\partial x^{3}\partial
y}+fv_{A}^{2}\frac{\partial^{3}u}{\partial t\partial
x^{2}}-f\lambda\frac{\partial^{2}u}{\partial
t^{2}}=\left(v_{A}^{2}\frac{\partial^{2}}{\partial
x^{2}}-\lambda\frac{\partial}{\partial t}\right)\left(f\frac{\partial
u}{\partial t}-C_{0}^{2}\frac{\partial^{2}u}{\partial x\partial y}\right).$
(A8)

Thankfully, the two terms in the right bracket can be eliminated simply by
rearranging equation (A2) in the following form

$\displaystyle f\frac{\partial u}{\partial
t}-C_{0}^{2}\frac{\partial^{2}u}{\partial x\partial
y}=v_{A}^{2}\frac{\partial^{2}v}{\partial x^{2}}-\frac{\partial^{2}v}{\partial
t^{2}}+C_{0}^{2}\frac{\partial^{2}v}{\partial^{2}y}-\frac{\partial^{2}V}{\partial
y\partial t}-\lambda\frac{\partial v}{\partial t}.$ (A9)

This way, we have finally eliminated all $u$ velocity terms. Inserting
Equation (A9) into Equation (A7) yields the wave equation (12) presented in
section 2.

##  Appendix B Derivation of the tidal potential

We consider a central body with radius $R_{0}$ with its center of mass being
the origin labeled $O$ and a tide generating body with mass $M_{\rm{t}}$ at a
distance $r=r(t)$ from the $O$, see Figure 14a. The tidal force originating
from $M_{\rm{t}}$ results from the combined action of the pseudoforce of
inertia due to the free fall motion of the central body around the center of
mass of the system and the gravitational force of the perturbing body. At
location $A$ on the surface of the central object the tidal force can be
computed from the vector difference of the gravitational pull the perturbing
body exerts on a test object at position $A$ and the gravitational pull the
perturbing body would exert on this object at the center of the central body
(which corresponds to the pseudo force of inertia of the free fall motion).

Figure 14:  (a) Sketch of a two body system with definition of quantities that
are involved in the calculation of the tidal potential caused by $M_{\rm{t}}$
on the surface of a central body. (b) Definition of angles in a spherical
system of reference.

The tidal force is a conservative force and can be expressed in terms of a
scalar potential, $\mbox{\boldmath$F$}_{\rm{t}}=-\nabla V_{\rm{t}}$, which can
be developped in terms of an infinite sum of Legendre polynomials labeled by
their degree $l$. Since the terms with $l=0$ and $l=1$ do not contribute to
the force and the relation of the radius $R_{0}$ to the distance $r$ is small,
only one term with $l=2$ has to be considered (see, e.g., Agnew (2015)). Then
the tidal potential is

$V_{\rm{t}}(A)=\frac{GM_{\rm{t}}R_{0}^{2}}{r^{3}(t)}P_{2}(\cos\alpha),$ (B1)

where $G$ is the gravitational constant and $P_{2}$ is the Legendre polynom
with degree $l=2$. The angle $\alpha$ represents the solid angle between the
line connecting the origin $O$ and the observation position $A$ and the line
connecting $O$ and the center of the tide generating object (Figure 14b). It
is convenient to transfer the problem into a spherical coordinate system,
where we have two sets of variables: $\varphi$ and $\theta$, which denote the
longitude and the latitude coordinate of the observation point $A$, and
$\varTheta$ and $\varPhi$, which describe the location of the tide generating
object $M_{\rm{t}}$. The Legendre polynom $P_{l}(\cos\alpha)$ can then be
expressed as a sum of associated Legendre polynomials in
$\vartheta,\varphi,\varTheta$ and $\varPhi$ via the expansion in terms of
spherical harmonics:

$P_{l}(\cos\alpha)=\frac{1}{2l+1}\sum\limits_{m=-l}^{l}\overline{Y_{l}^{m}}(\vartheta,\varphi)Y_{l}^{m}(\varTheta,\varPhi)\mbox{
with
}Y_{l}^{m}(\vartheta,\varphi)=(-1)^{m}\left[\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}\right]^{{1}/{2}}P_{l}^{m}(\cos\vartheta)e^{im\varphi},$
(B2)

where $P_{l}^{m}$ are the associated Legendre functions taken from Munk &
Cartwright (1966). Since we only perform a decomposition for $l=2$, we only
need the following five associated Legendre polynomials

$P_{2}^{-2}=\frac{1}{8}\sin^{2}\vartheta,\qquad
P_{2}^{-1}=\frac{1}{2}\sin\vartheta\cos\vartheta,\qquad
P_{2}^{0}=\frac{1}{2}(3\cos^{2}\vartheta-1),\qquad
P_{2}^{1}=-3\sin\vartheta\cos\vartheta,\qquad P_{2}^{2}=3\sin^{2}\vartheta,$
(B3)

so that after some calculations using standard relations for the trigonometric
functions we obtain

$P_{2}(\cos\alpha)=\frac{3}{4}\left[\frac{1}{3}(3\cos^{2}\varTheta-1)(3\cos^{2}\vartheta-1)+\sin
2\varTheta\sin
2\vartheta\cos(\varphi-\varPhi)+\sin^{2}\varTheta\sin^{2}\vartheta\cos(2\varphi-2\varPhi)\right].$
(B4)

We now consider a tide generating object that moves exactly in the plane
defined by the equatorial plane of the central body so that
$\varTheta=90^{\circ}$. Thus the second term in the brackets on the right side
of Eq. (B4) ($\propto\sin 2\varTheta$) vanishes and the first term simplyfies
to $\sin^{2}\vartheta$ (after dropping the constant terms that do not
contribute to the tidal force. We further assume that the distance between
central body and perturbing body is constant, i.e. $r(t)=r$. In reality the
two bodies propagate around a common center of mass on a path that is the
solution of the corresponding Kepler problem. The resulting variation of the
distance between both bodies introduces a slow time scale for the amplitude of
the tidal force, which can be neglected with regard to the current problem.
Hence, the tidal potential finally reads:

$V_{\rm{t}}(R,\vartheta,\varphi,t)=\underbrace{\frac{3}{4}\frac{GM_{\rm{t}}}{R}\left(\frac{R}{r}\right)^{3}}_{\displaystyle
K}\left[\sin^{2}\vartheta\left[1+\cos\left(2\varphi-2\varPhi\right)\right]\right],$
(B5)

where we have introduced the constant term $K$. The only time-dependent
quantity in Equation (B5) is the angle $\varPhi$, which describes the
longitude of the tide generating body with respect to the origin. The time-
dependence of $\varPhi$ is simply given by
$\varPhi(t)=\left(\varOmega_{\rm{t}}-\varOmega_{\rm{0}}\right)t$ where
$\varOmega_{\rm{0}}$ denotes the angular frequency of the rotation of the
central body (i.e. the sun in our particular case) and $\Omega_{\rm{t}}$ the
angular frequency of the motion of the perturbing body around the central
object. In the following we translate the angular coordinates $\varphi$ and
$\vartheta$ into the Cartesian coordinates $x$ and $y$ as they appear in the
$\beta$-plane model. We consider a small regime around a fixed latitudinal
coordinate $\vartheta_{0}$ so that we can write
$\vartheta=\vartheta_{0}+\vartheta^{\prime}$ with
$\vartheta^{\prime}=y/R_{0}\ll 1$. The corresponding expansions of the
expression that describes the latitudinal dependence is then

$\sin^{2}(\vartheta)=\sin^{2}(\vartheta_{0}+\vartheta^{\prime})\approx\sin^{2}(\vartheta_{0})+\sin(2\vartheta_{0})\frac{y}{R_{0}}+\cos(2\vartheta_{0})\frac{y^{2}}{R_{0}^{2}}+...$
(B6)

In order to keep the same approximation used for the $\beta$-plane approach we
ignore the term quadratic in $y$. We examine two particular cases. At the
equator we have $\vartheta_{0}=\pi/2$ so that $\sin^{2}(\vartheta_{0})=1$ and
$\sin(2\vartheta_{0})=0$ and we end up with a forcing independent of $y$.

$V_{\rm{t}}^{\rm{eq}}(x,y,t)=K\left[1+\cos\left[2x-2\left(\varOmega_{\rm{t}}-\varOmega_{\rm{0}}\right)t\right]\right].$
(B7)

For non-equatorial solutions it is reasonable to choose $\vartheta_{0}=\pi/4$
so that $\sin^{2}(\vartheta_{0})=1/2$ and $\sin(2\vartheta_{0})=1$ and the
corresponding tidal potential reads

$V_{\rm{t}}^{45^{\circ}}(x,y,t)=K\left(\frac{1}{2}+\frac{y}{R_{0}}\right)\left[1+\cos\left[2x-2\left(\varOmega_{\rm{t}}-\varOmega_{\rm{0}}\right)t\right]\right].$
(B8)

##  Appendix C Forced wave solutions at the equator

The forced wave problem is solved by inserting the ansatz

$v=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\alpha_{m,n}(t)\exp\left(\frac{im}{R_{0}}x\right)\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{n}(\sqrt{\mu}y),$ (C1)

and the coordinate expansion

$y=\sum_{n=1}^{\infty}\frac{2^{\frac{5}{2}-2n}}{\sqrt{\mu}(n-1)!}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y).$ (C2)

into the wave equation (53). The comparison of the modal coefficients
$\alpha_{m,n}(t)$ with the forcing potential $\sim\cos(2x/R_{0}-2\Omega t)$
directly reveals the zonal wave number to be fixed at $m=2/R_{0}$, which
describes planetary waves with two crests and troughs fitting around the
equator. All other modal coefficients are not subject to any forcing and
consequently obey damped wave equations, letting them decay exponentially in
time. Linear stationary solutions can only be found for $m=2$. In the same way
we can argue that latitudinal wave numbers corresponding to stationary
solutions must be uneven $n\rightarrow 2n+1$ since only uneven Hermitian
polynomials appear in the expansion (C2). The remaining modal coefficients
must satisfy the conditional equation

$\displaystyle\sum_{n=1}^{\infty}\ddddot{\alpha}_{2,2n-1}(t)\frac{4\Omega^{2}}{C_{0}^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\dddot{\alpha}_{2,2n-1}(t)\frac{8\lambda\Omega^{2}}{C_{0}^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\ddot{\alpha}_{2,2n-1}(t)\frac{32v_{A}^{2}\Omega^{2}}{C_{0}^{2}R_{0}^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\alpha_{2,2n-1}(t)\frac{64v_{A}^{2}\Omega^{2}}{R_{0}^{4}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\alpha_{2,2n-1}(t)\frac{64v_{A}^{4}\Omega^{2}}{C_{0}^{2}R_{0}^{4}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\ddot{\alpha}_{2,2n-1}(t)\frac{16\Omega^{2}}{R_{0}^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle-$
$\displaystyle\sum_{n=1}^{\infty}\ddot{\alpha}_{2,2n-1}(t)\left[4\Omega^{2}+2\lambda\Omega
i-\frac{4v_{A}^{2}}{R_{0}^{2}}\right]{\color[rgb]{0,0,1}\frac{\partial^{2}}{\partial
y^{2}}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\ddot{\alpha}_{2,2n-1}(t){\color[rgb]{0,0,1}\frac{4\beta^{2}\Omega^{2}}{C_{0}^{2}}y^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\ddot{\alpha}_{2,2n-1}(t)\frac{32\lambda\Omega^{2}v_{A}^{2}}{C_{0}^{2}R_{0}^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\dot{\alpha}_{2,2n-1}(t)\frac{16\lambda\Omega^{2}}{R_{0}^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\ddot{\alpha}_{2,2n-1}(t)\frac{4\lambda^{2}\Omega^{2}}{C_{0}^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\dot{\alpha}_{2,2n-1}(t)\frac{8i\beta\Omega^{2}}{R_{0}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}\right)$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{32K\beta\Omega^{4}}{C_{0}^{2}R}\frac{2^{\frac{5}{2}-2n}}{\sqrt{\mu}(n-1)!}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)\exp\left(i\frac{2x}{R_{0}}-i2\Omega
t\right).$ (C3)

We have used $\alpha_{2,2n-1}=i\dot{\alpha}_{2,2n-1}/2\Omega$ and
$\alpha_{2,2n-1}=-\ddot{\alpha}_{2,2n-1}/4\Omega^{2}$ to evaluate the three
terms involving $y$ derivatives. From equation (C3) we want to extract
ordinary differential equations for the modal coefficients
$\alpha_{2,2n-1}(t)$. All terms are proportional to $\exp(i2x/R_{0})$ so that
the $x$ coordinate drops out by default. This is unfortunately not the case
for the $y$ coordinate since the blue-marked terms do not fit with the
Hermitian base. However, we can bypass this issue by applying the identity

$\displaystyle\mu^{2}y^{2}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{n}(\sqrt{\mu}y)-\frac{\partial^{2}}{\partial
y^{2}}\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{n}(\sqrt{\mu}y)=(2n+1)\mu\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{n}(\sqrt{\mu}y)$ (C4)

that projects the two terms on the Hermitian eigenbasis. Inserting $\mu$ (47)
allows to rearrange (C4) into

$\displaystyle\left[\frac{4\Omega^{2}+2i\lambda\Omega\beta^{2}}{C_{0}^{2}}y^{2}-\left(4\Omega^{2}-\frac{4v_{A}^{2}}{R_{0}^{2}}\right)\frac{\partial^{2}}{\partial
y^{2}}\right]\exp\left(-\frac{\mu y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)$
$\displaystyle=\left[\frac{2\Omega\beta}{C_{0}}\sqrt{4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}}(4n-1)\right]\exp\left(-\frac{\mu
y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)$ (C5)

we insert into (C3) to replace the blue marked terms, yielding

$\displaystyle\sum_{n=1}^{\infty}\Bigg{[}\ddddot{\alpha}_{2,2n-1}(t)+2\lambda\dddot{\alpha}_{2,2n-1}(t)+\left(\frac{C_{0}\beta}{2\Omega}\sqrt{4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}}(4n-1)+\frac{4}{R_{0}^{2}}(C_{0}^{2}+2v_{A}^{2})+\frac{8\lambda
v_{A}^{2}}{R_{0}^{2}}+\lambda^{2}\right)\ddot{\alpha}_{2,2n-1}(t)$
$\displaystyle+\left(\frac{4\lambda
C_{0}^{2}}{R_{0}^{2}}-\frac{2iC_{0}^{2}\beta}{R_{0}}\right)\dot{\alpha}_{2,2n-1}(t)+\left(\frac{16v_{A}^{2}C_{0}^{2}}{R_{0}^{4}}+\frac{16v_{A}^{4}}{R_{0}^{4}}\right)\alpha_{2,2n-1}(t)$
$\displaystyle-\frac{8K\beta\Omega^{2}}{R}\frac{2^{\frac{5}{2}-2n}}{\sqrt{\mu}(n-1)!}\exp(-i2\Omega
t)\Bigg{]}\exp\left(-\frac{\mu y^{2}}{2}\right)H_{2n-1}(\sqrt{\mu}y)=0$ (C6)

The infinite sums can only yield zero if all individual summands disappear.
Therefore, each coefficient $\alpha_{2,2n-1}(t)$ must satisfy the following
set of ordinary differential equations

$\displaystyle\ddddot{\alpha}_{2,2n-1}(t)+2\lambda\dddot{\alpha}_{2,2n-1}(t)+\left(\frac{C_{0}\beta}{2\Omega}\sqrt{4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}}(4n-1)+\frac{4}{R_{0}^{2}}(C_{0}^{2}+2v_{A}^{2})+\frac{8\lambda
v_{A}^{2}}{R_{0}^{2}}+\lambda^{2}\right)\ddot{\alpha}_{2,2n-1}(t)$
$\displaystyle+\left(\frac{4\lambda
C_{0}^{2}}{R_{0}^{2}}-\frac{2iC_{0}^{2}\beta}{R_{0}}\right)\dot{\alpha}_{2,2n-1}(t)+\left(\frac{16v_{A}^{2}C_{0}^{2}}{R_{0}^{4}}+\frac{16v_{A}^{4}}{R_{0}^{4}}\right)\alpha_{2,2n-1}(t)-\frac{8K\beta\Omega^{2}}{R}\frac{2^{\frac{5}{2}-2n}}{\sqrt{\mu}(n-1)!}\exp(-i2\Omega
t)=0$ (C7)

We have now reduced the partial differential Equation (53) into an infinite
set of decoupled ordinary differential equations to be solved readily.
Equation (C7) has the following stationary solution:

$\displaystyle\alpha_{2,2n-1}(t)=\frac{8iK\beta\Omega^{2}2^{\frac{5}{2}-2n}(R\sqrt{\mu}(n-1)!)^{-1}\exp(-i2\Omega
t)}{\left[\left(4\Omega^{2}+2i\lambda\Omega-\frac{4v_{A}^{2}}{R_{0}^{2}}\right)\left(4\Omega^{2}+2i\lambda\Omega-\frac{4}{R_{0}^{2}}(C_{0}^{2}+v_{A}^{2})\right)-C_{0}^{2}\frac{4\beta}{R_{0}}\Omega-
C_{0}\beta|2\Omega|\sqrt{4\Omega^{2}+2i\lambda\omega-\frac{4v_{A}^{2}}{R_{0}^{2}}}(4n-1)\right]}.$
(C8)

Inserting the modal coefficients back into the ansatz (C1) yields the solution
(57) we present in section 3.3.1.

##  Appendix D Forced wave solutions at mid latitudes

We solve the forced wave problem by inserting the ansatz

$v=\sum_{m=0}^{\infty}\sum_{n=0}^{\infty}\left[\alpha_{m,n}(t)\sin\left(\frac{m}{R_{0}}x\right)+\beta_{m,n}(t)\cos\left(\frac{m}{R_{0}}x\right)\right]\exp\left(\frac{in\pi}{R_{0}}y\right)$
(D1)

to be confined within a mid-latitude wave channel defined by the interval
$-R_{0}/2\leq y\leq R_{0}/2$, into the wave equation (58). The forcing terms
on the right-hand side of (58) involve both constant and $y$-proportional
terms, which both must be expanded as Fourier series (59) and (60) in order to
project them onto the harmonic basis used in the ansatz. The direct comparison
of the modal coefficients $\alpha_{m,n}(t)$ and $\beta_{m,n}(t)$ with the
forcing potential $\sim\cos(2x/R_{0}-2\Omega t)$ directly reveals the zonal
wave number to be fixed at $m=2/R_{0}$, which describes planetary waves with
two crests and troughs fitting around the equator. All other modal
coefficients are not subject to any forcing and consequently obey freely
damped wave equations, letting them decay exponentially in time. Linear
stationary solutions can only be found for $m=2$. In the same way we can argue
that latitudinal wave numbers corresponding to stationary solutions must be
uneven $n\rightarrow 2n+1$ since only uneven terms appear in the Fourier
expansion (59) and (60). The remaining modal coefficients must satisfy the
conditional equation

$\displaystyle\sum_{n=1}^{\infty}\left[\ddddot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\ddddot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\dddot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\dddot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]2\lambda\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\ddot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\ddot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{8v_{A}^{2}}{R_{0}^{2}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\ddot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\ddot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{4C_{0}^{2}}{R_{0}^{2}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\ddot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\ddot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{C_{0}^{2}(2n-1)^{2}\pi^{2}}{R_{0}^{2}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\ddot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\ddot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]f_{0}^{2}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\ddot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\ddot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\lambda^{2}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\dot{\beta}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)-\dot{\alpha}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{2C_{0}^{2}\beta}{R_{0}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\dot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\dot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{8\lambda
v_{A}^{2}}{R_{0}^{2}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\dot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\dot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{4\lambda
C_{0}^{2}}{R_{0}^{2}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\dot{\alpha}_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\dot{\beta}_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{(2n-1)^{2}\pi^{2}\lambda
C_{0}^{2}}{R_{0}^{2}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\alpha_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\beta_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{16C_{0}^{2}v_{A}^{2}}{R_{0}^{4}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\alpha_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\beta_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{4(2n-1)^{2}\pi^{2}C_{0}^{2}v_{A}^{2}}{R_{0}^{4}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\alpha_{2,2n-1}(t)\sin\left(\frac{2x}{R_{0}}\right)+\beta_{2,2n-1}(t)\cos\left(\frac{2x}{R_{0}}\right)\right]\frac{16v_{A}^{4}}{R_{0}^{4}}\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\left[\sin\left(\frac{2x}{R_{0}}\right)\cos\left(2\Omega
t\right)-\cos\left(\frac{2x}{R_{0}}\right)\sin\left(2\Omega
t\right)\right]\left(f_{0}\Omega-\frac{2if_{0}\Omega}{(2n-1)\pi}+2\Omega^{2}-\frac{2v_{A}^{2}}{R_{0}^{2}}\right)\frac{16K\Omega}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi}\exp\left(\frac{i(2n-1)\pi
y}{R_{0}}\right)$ $\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\left[\sin\left(\frac{2x}{R_{0}}\right)\sin\left(2\Omega
t\right)+\cos\left(\frac{2x}{R_{0}}\right)\cos\left(2\Omega
t\right)\right]\frac{16K\lambda\Omega^{2}}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi}\exp\left(\frac{i(2n-1)\pi
y}{R_{0}}\right)$ (D2)

All summands are either proportional to $\sim\sin(2x/R_{0})\exp(i(2n-1)\pi
y/R_{0})$ or to $\sim\cos(2x/R_{0})\exp(i(2n-1)\pi y/R_{0})$ so that we can
rearrange the terms as follows:

$\displaystyle\sum_{n=1}^{\infty}\Bigg{[}\ddddot{\alpha}_{2,2n-1}(t)+2\lambda\dddot{\alpha}_{2,2n-1}(t)+\left(\frac{4C_{0}^{2}+8v_{A}^{2}+(2n-1)^{2}\pi^{2}C_{0}^{2}}{R_{0}^{2}}+f_{0}^{2}+\lambda^{2}\right)\ddot{\alpha}_{2,2n-1}(t)$
$\displaystyle+$ $\displaystyle\left(\frac{4\lambda C_{0}^{2}+8\lambda
v_{A}^{2}+(2n-1)^{2}\pi^{2}\lambda
C_{0}^{2}}{R_{0}^{2}}\right)\dot{\alpha}_{2,2n-1}(t)+\frac{2\beta
C_{0}^{2}}{R_{0}}\dot{\beta}_{2,2n-1}(t)+\left(\frac{16v_{A}^{2}C_{0}^{2}+16v_{A}^{4}+4(2n-1)^{2}\pi^{2}v_{A}^{2}C_{0}^{2}}{R_{0}^{4}}\right)\alpha_{2,2n-1}(t)$
$\displaystyle-$
$\displaystyle\left(f_{0}\Omega-\frac{2if_{0}\Omega}{(2n-1)\pi}+2\Omega^{2}-\frac{2v_{A}^{2}}{R_{0}^{2}}\right)\frac{16K\Omega}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi}\cos\left(2\Omega
t\right)-\frac{16K\lambda\Omega^{2}}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi}\sin\left(2\Omega
t\right)\Bigg{]}$ $\displaystyle\times$
$\displaystyle\sin\left(\frac{2}{R_{0}}x\right)\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)$
$\displaystyle+$
$\displaystyle\sum_{n=1}^{\infty}\Bigg{[}\ddddot{\beta}_{2,2n-1}(t)+2\lambda\dddot{\beta}_{2,2n-1}(t)+\left(\frac{4C_{0}^{2}+8v_{A}^{2}+(2n-1)^{2}\pi^{2}C_{0}^{2}}{R_{0}^{2}}+f_{0}^{2}+\lambda^{2}\right)\ddot{\beta}_{2,2n-1}(t)$
$\displaystyle+$ $\displaystyle\left(\frac{4\lambda C_{0}^{2}+8\lambda
v_{A}^{2}+(2n-1)^{2}\pi^{2}\lambda
C_{0}^{2}}{R_{0}^{2}}\right)\dot{\beta}_{2,2n-1}(t)-\frac{2\beta
C_{0}^{2}}{R_{0}}\dot{\alpha}_{2,2n-1}(t)+\left(\frac{16v_{A}^{2}C_{0}^{2}+16v_{A}^{4}+4(2n-1)^{2}\pi^{2}v_{A}^{2}C_{0}^{2}}{R_{0}^{4}}\right)\alpha_{2,2n-1}(t)$
$\displaystyle+$
$\displaystyle\left(f_{0}\Omega-\frac{2if_{0}\Omega}{(2n-1)\pi}+2\Omega^{2}-\frac{2v_{A}^{2}}{R_{0}^{2}}\right)\frac{16K\Omega}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi}\sin\left(2\Omega
t\right)-\frac{16K\lambda\Omega^{2}}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi}\cos\left(2\Omega
t\right)\Bigg{]}$ $\displaystyle\times$
$\displaystyle\cos\left(\frac{2}{R_{0}}x\right)\exp\left(\frac{i(2n-1)\pi}{R_{0}}y\right)=0$
(D3)

The infinite sums can only yield zero if all individual summands disappear
independently. Therefore, each pair of coefficients $\alpha_{2,2n-1}(t)$ and
$\beta_{2,2n-1}(t)$ must satisfy the following set of coupled ordinary
differential equations

$\displaystyle\ddddot{\alpha}_{2,2n-1}(t)+2\lambda\dddot{\alpha}_{2,2n-1}(t)+\tilde{\Lambda}_{1}\ddot{\alpha}_{2,2n-1}(t)+\tilde{\Lambda}_{2}\dot{\alpha}_{2,2n-1}(t)+\frac{2\beta
C_{0}^{2}}{R_{0}}\dot{\beta}_{2,2n-1}(t)+\tilde{\Lambda}_{3}\alpha_{2,2n-1}(t)$
$\displaystyle-\tilde{K}_{1}\cos\left(2\Omega
t\right)-\tilde{K}_{2}\sin\left(2\Omega t\right)=0,$ (D4)
$\displaystyle\ddddot{\beta}_{2,2n-1}(t)+2\lambda\dddot{\beta}_{2,2n-1}(t)+\tilde{\Lambda}_{1}\ddot{\beta}_{2,2n-1}(t)+\tilde{\Lambda}_{2}\dot{\beta}_{2,2n-1}(t)-\frac{2\beta
C_{0}^{2}}{R_{0}}\dot{\alpha}_{2,2n-1}(t)+\tilde{\Lambda}_{3}\beta_{2,2n-1}(t)$
$\displaystyle+\tilde{K}_{1}\sin\left(2\Omega
t\right)-\tilde{K}_{2}\cos\left(2\Omega t\right)=0,$ (D5)

where we have introduced the coefficients $\tilde{\Lambda}_{1}$,
$\tilde{\Lambda}_{2}$, $\tilde{\Lambda}_{3}$, $\tilde{K}_{1}$ and
$\tilde{K}_{2}$ for better readability, which are given as

$\displaystyle\tilde{\Lambda}_{1}=\frac{4C_{0}^{2}+8v_{A}^{2}+(2n-1)^{2}\pi^{2}C_{0}^{2}}{R_{0}^{2}}+f_{0}^{2}+\lambda^{2},$
$\displaystyle\tilde{\Lambda}_{2}=\frac{4\lambda C_{0}^{2}+8\lambda
v_{A}^{2}+(2n-1)^{2}\pi^{2}\lambda C_{0}^{2}}{R_{0}^{2}},$
$\displaystyle\tilde{\Lambda}_{3}=\frac{16v_{A}^{2}C_{0}^{2}+16v_{A}^{4}+4(2n-1)^{2}\pi^{2}v_{A}^{2}C_{0}^{2}}{R_{0}^{4}},$
$\displaystyle\tilde{K}_{1}=\left(f_{0}\Omega-\frac{2if_{0}\Omega}{(2n-1)\pi}+2\Omega^{2}-\frac{2v_{A}^{2}}{R_{0}^{2}}\right)\frac{16K\Omega}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi},\
\ \
\tilde{K}_{2}=\frac{16K\lambda\Omega^{2}}{R_{0}}\frac{(-1)^{n-1}}{(2n-1)\pi}.$

We have now reduced the partial differential Equation (58) into an infinite
set of ordinary differential equations to be solved readily. It is convenient
to introduce auxiliary variables

$a_{1}=\alpha_{2,2n-1}+i\beta_{2,2n-1},~{}~{}~{}a_{2}=\alpha_{2,2n-1}-i\beta_{2,2n-1},$
(D6)

allowing to decouple the Equations (D4) and (D5) into the following form:

$\displaystyle\ddddot{a_{1}}+2\lambda\dddot{a_{1}}+\bar{\Lambda}_{1}\ddot{a}_{1}+\left(\bar{\Lambda}_{2}-\frac{2i\beta
C_{0}^{2}}{R_{0}}\right)\dot{a}_{1}+\bar{\Lambda}_{3}a_{1}=(\bar{K}_{1}+i\bar{K}_{2})\exp(-2i\Omega
t),$ (D7)
$\displaystyle\ddddot{a_{2}}+2\lambda\dddot{a_{2}}+\bar{\Lambda}_{1}\ddot{a}_{2}+\left(\bar{\Lambda}_{2}+\frac{2i\beta
C_{0}^{2}}{R_{0}}\right)\dot{a}_{2}+\bar{\Lambda}_{3}a_{2}=(\bar{K}_{1}-i\bar{K}_{2})\exp(2i\Omega
t).$ (D8)

The particular, stationary solutions of these equations can be readily
obtained

$\displaystyle
a_{1}=\frac{\bar{K}_{1}+i\bar{K}_{2}}{16\Omega^{4}-4\Omega^{2}\bar{\Lambda}_{1}+\bar{\Lambda}_{3}-\frac{4\beta
C_{0}^{2}\Omega}{R_{0}}+2i\Omega(8\lambda\Omega^{2}-\bar{\Lambda}_{2})}\exp(-2i\Omega
t)$ (D9) $\displaystyle
a_{2}=\frac{\bar{K}_{1}-i\bar{K}_{2}}{16\Omega^{4}-4\Omega^{2}\bar{\Lambda}_{1}+\bar{\Lambda}_{3}-\frac{4\beta
C_{0}^{2}\Omega}{R_{0}}-2i\Omega(8\lambda\Omega^{2}-\bar{\Lambda}_{2})}\exp(2i\Omega
t).$ (D10)

As the last step, it remains to express the original coefficients
$\alpha_{2,2n-1}$ and $\beta_{2,2n-1}$ as functions of $a_{1}$ and $a_{2}$ as

$\alpha_{2,2n-1}=\frac{a_{1}+a_{2}}{2},~{}~{}~{}\beta_{2,2n-1}=\frac{a_{1}-a_{2}}{2i}.$
(D11)

Inserting these modal coefficients back into the ansatz (D1) finally yields
the solution (62) we present in section 3.3.2.

## References

  * Abreu et al. (2012) Abreu, J. A., Beer, J., Ferriz-Mas, A., McCracken, K. G., & Steinhilber, F. 2012, A&A, 548, 9, doi: 10.1051/0004-6361/201219997
  * Agnew (2015) Agnew, D. C. 2015, Treatise on Geophysics, 2nd edn. (Elsevier), 151–178
  * Bildsten & Ushomirsky (2000) Bildsten, L., & Ushomirsky, G. 2000, ApJ, 529, L33, doi: 10.1086/312454
  * Chang & Lim (1982) Chang, C.-P., & Lim, H. 1982, J. Atmos. Sci., 39, 1726, doi: 10.1175/1520-0469(1982)039<1726:OTEOVD>2.0.CO;2
  * Charbonneau (2022) Charbonneau, P. 2022, Front. Astron. Space Sci., 9:853676, doi: 10.3389/fspas.2022.853676
  * Dikpati et al. (2018a) Dikpati, M., Belucz, B., Gilman, P. A., & McIntosh, S. W. 2018a, ApJ, 862, 159, doi: 10.3847/1538-4357/aacefa
  * Dikpati et al. (2017) Dikpati, M., Cally, P. S., McIntosh, S. W., & Heifetz, E. 2017, Scientific Reports, 7, 14750, doi: 10.1038/s41598-017-14957-x
  * Dikpati & Gilman (2001) Dikpati, M., & Gilman, P. A. 2001, ApJ, 551, 536, doi: 10.1086/320080
  * Dikpati et al. (2020) Dikpati, M., Gilman, P. A., Chatterjee, S., McIntosh, S. W., & Zaqarashvili, T. V. 2020, ApJ, 896, 141, doi: 10.3847/1538-4357/ab8b63
  * Dikpati & McIntosh (2020) Dikpati, M., & McIntosh, S. W. 2020, Space Weather, 18, e2018SW002109, doi: 10.1029/2018SW002109
  * Dikpati et al. (2018b) Dikpati, M., McIntosh, S. W., Bothun, G., et al. 2018b, The Astrophysical Journal, 853, 144, doi: 10.3847/1538-4357/aaa70d
  * Fedotova et al. (2021) Fedotova, M., Klimachkov, D., & Petrosyan, A. 2021, Universe, 7, 87, doi: 10.3390/universe7040087
  * Ferriz-Mas (1996) Ferriz-Mas, A. 1996, ApJ, 458, 802, doi: 10.1086/176859
  * Gill (1980) Gill, A. E. 1980, Quarterly Journal of the Royal Meteorological Society, 106, 447, doi: 10.1002/qj.49710644905
  * Gilman (1968) Gilman, P. A. 1968, Science, 160, 760, doi: 10.1126/science.160.3829.760.b
  * Gilman (2000) —. 2000, ApJL, 544, L79, doi: 10.1086/317291
  * Gilman & Dikpati (2002) Gilman, P. A., & Dikpati, M. 2002, ApJ, 576, 1031, doi: 10.1086/341799
  * Gizon et al. (2020) Gizon, L., Fournier, D., & Albekioni, M. 2020, A&A, 642, 10, doi: 10.1051/0004-6361/202038525
  * Horstmann et al. (2021) Horstmann, G. M., Anders, S., Kelley, D. H., & Weier, T. 2021, J. Fluid Mech., 925, doi: 10.1017/jfm.2021.686
  * Horstmann et al. (2020) Horstmann, G. M., Herreman, W., & Weier, T. 2020, J. Fluid Mech., 891, A22, doi: 10.1017/jfm.2020.163
  * Hunter (2015) Hunter, S. 2015, PhD thesis, University of Leeds
  * Klimachkov & Petrosyan (2016) Klimachkov, D. A., & Petrosyan, A. S. 2016, J. Exp. Theor. Phys., 122, 832, doi: 10.1134/S1063776116040178
  * Klimachkov & Petrosyan (2017) —. 2017, Soviet Journal of Experimental and Theoretical Physics, 125, 597, doi: 10.1134/S1063776117090059
  * Li et al. (2006) Li, L., Ingersoll, A. P., Vasavada, A. R., et al. 2006, Icarus, 185, 416, doi: 10.1016/j.icarus.2006.08.005
  * Liang et al. (2019) Liang, Z.-C., Gizon, L., Birch, A. C., & Duvall Jr., T. L. 2019, A&A, 626, 11, doi: 10.1051/0004-6361/201834849
  * Longuet-Higgins (1964) Longuet-Higgins, H. C. 1964, Proc. R. Soc. London, Ser. A, 279, 446, doi: 10.1098/rspa.1964.0116
  * Löptien et al. (2018) Löptien, B., Gizon, L., Birch, A. C., et al. 2018, Nat Astron, 2, 568, doi: 10.1038/s41550-018-0460-x
  * Lyra & Umurhan (2019) Lyra, W., & Umurhan, O. M. 2019, Publ. Astron. Soc. Pac., 131, 072001, doi: 10.1088/1538-3873/aaf5ff
  * Mandal & Hanasoge (2020) Mandal, K., & Hanasoge, S. 2020, ApJ, 891, 125, doi: 10.3847/1538-4357/ab7227
  * Matsuno (1966) Matsuno, T. 1966, Journal of the Meteorological Society of Japan. Ser. II, 44, 25, doi: 10.2151/jmsj1965.44.1_25
  * McIntosh et al. (2017) McIntosh, S. W., Cramer, W. J., Pichardo Marcano, M., & Leamon, R. J. 2017, Nature Astronomy, 1, 1, doi: 10.1038/s41550-017-0086
  * Mofjeld (1981) Mofjeld, H. O. 1981, J. Phys. Oceanogr., 11, 1585, doi: 10.1175/1520-0485(1981)011<1585:AATOHF>2.0.CO;2
  * Munk & Cartwright (1966) Munk, W. H., & Cartwright, D. E. 1966, Philos. Trans. R. Soc. Lond. A, 259, 533, doi: 10.1098/rsta.1966.0024
  * Okhlopkov (2016) Okhlopkov, V. P. 2016, Moscow Univ. Phys., 71, 440, doi: 10.3103/S0027134916040159
  * Pedlosky (1987) Pedlosky, J. 1987, Geophysical Fluid Dynamics, 2nd edn. (New York: Springer-Verlag), doi: 10.1007/978-1-4612-4650-3
  * Petrosyan et al. (2020) Petrosyan, A., Klimachkov, D., Fedotova, M., & Zinyakov, T. 2020, Atmosphere, 11, 314, doi: 10.3390/atmos11040314
  * Proxauf et al. (2020) Proxauf, B., Gizon, L., Löptien, B., et al. 2020, A&A, 634, 16, doi: 10.1051/0004-6361/201937007
  * Raphaldini & Raupp (2015) Raphaldini, B., & Raupp, C. F. M. 2015, The Astrophysical Journal, 799, 78, doi: 10.1088/0004-637X/799/1/78
  * Raphaldini et al. (2019) Raphaldini, B., Teruya, A. S., Raupp, C. F. M., & Bustamante, M. D. 2019, ApJ, 887, 1, doi: 10.3847/1538-4357/ab5067
  * Read et al. (2006) Read, P., Dowling, T., & Schubert, G. 2006, Nature, 460, 608, doi: 10.1038/nature08194
  * Rossby (1939) Rossby, C.-G. 1939, J. Mar. Res., 2, 38
  * Rüdiger (1989) Rüdiger, G. 1989, Differential Rotation and Stellar Convection (Berlin: Akademie Verlag), doi: 10.1515/9783112532126
  * Schecter et al. (2001) Schecter, D. A., Boyd, J. F., & Gilman, P. A. 2001, ApJ, 551, L185, doi: 10.1086/320027
  * Stefani et al. (2016) Stefani, F., Giesecke, A., Weber, N., & Weier, T. 2016, Sol Phys, 291, 2197, doi: 10.1007/s11207-016-0968-0
  * Stefani et al. (2018) —. 2018, Sol Phys, 293, 12, doi: 10.1007/s11207-017-1232-y
  * Stefani et al. (2019) Stefani, F., Giesecke, A., & Weier, T. 2019, Sol Phys, 294, 60, doi: 10.1007/s11207-019-1447-1
  * Stefani et al. (2021) Stefani, F., Stepanov, R., & Weier, T. 2021, Sol Phys, 296, 88, doi: 10.1007/s11207-021-01822-4
  * Tsai et al. (2014) Tsai, S.-M., Dobbs-Dixon, I., & Gu, P.-G. 2014, ApJ, 793, 141, doi: 10.1088/0004-637X/793/2/141
  * Weber et al. (2015) Weber, N., Galindo, V., Stefani, F., & Weier, T. 2015, New J. Phys., 17, 113013, doi: 10.1088/1367-2630/17/11/113013
  * Wu et al. (2001) Wu, Z., Sarachik, E. S., & Battisti, D. S. 2001, J. Atmos. Sci., 58, 724, doi: 10.1175/1520-0469(2001)058<0724:TDTCUR>2.0.CO;2
  * Yamagata & Philander (1985) Yamagata, T., & Philander, S. G. H. 1985, J. Oceanogr. Soc. Jpn., 41, 345, doi: 10.1007/BF02109241
  * Zaqarashvili (2018) Zaqarashvili, T. V. 2018, ApJ, 856, 32, doi: 10.3847/1538-4357/aab26f
  * Zaqarashvili et al. (2010a) Zaqarashvili, T. V., Carbonell, M., Oliver, R., & Ballester, J. L. 2010a, ApJ, 709, 749, doi: 10.1088/0004-637X/709/2/749
  * Zaqarashvili et al. (2010b) —. 2010b, ApJL, 724, L95, doi: 10.1088/2041-8205/724/1/L95
  * Zaqarashvili et al. (2009) Zaqarashvili, T. V., Oliver, R., & Ballester, J. L. 2009, ApJ, 691, L41, doi: 10.1088/0004-637X/691/1/L41
  * Zaqarashvili et al. (2007) Zaqarashvili, T. V., Oliver, R., Ballester, J. L., & Shergelashvili, B. M. 2007, A&A, 470, 815, doi: 10.1051/0004-6361:20077382
  * Zaqarashvili et al. (2015) Zaqarashvili, T. V., Oliver, R., Hanslmeier, A., et al. 2015, ApJL, 805, L14, doi: 10.1088/2041-8205/805/2/L14
  * Zaqarashvili et al. (2021) Zaqarashvili, T. V., Albekioni, M., Ballester, J. L., et al. 2021, Space Sci Rev, 217, 15, doi: 10.1007/s11214-021-00790-2