<figure><img src="image_1.png" /></figure>

HAL
open science

# Layered semi-convection and tides in giant planet interiors. I. Propagation of internal waves

Q. André, A. Barker, S. Mathis

<figure><img src="image_2.png" /></figure>

► To cite this version:

Q. André, A. Barker, S. Mathis. Layered semi-convection and tides in giant planet interiors. I. Propagation of internal waves. Astronomy & Astrophysics - A&A, 2017, 605, pp.A117. 10.1051/0004-6361/201730765. obspm-02193087

HAL Id: obspm-02193087

https://hal-obspm.ccsd.cnrs.fr/obspm-02193087v1

Submitted on 21 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Layered semi-convection and tides in giant planet interiors

## I. Propagation of internal waves

Q. André<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, A. J. Barker<fnref n="2" />,<fnref n="4" />, and S. Mathis<fnref n="1" />,<fnref n="5" />

Received 10 March 2017 / Accepted 27 April 2017

## ABSTRACT

*Context.* Layered semi-convection is a possible candidate to explain Saturn's luminosity excess and the abnormally large radius of some hot Jupiters. In giant planet interiors, it could lead to the creation of density staircases, which are convective layers separated by thin stably stratified interfaces. These are also observed on Earth in some lakes and in the Arctic Ocean.

*Aims.* We aim to study the propagation of internal waves in a region of layered semi-convection, with the aim to predict energy transport by internal waves incident upon a density staircase. The goal is then to understand the resulting tidal dissipation when these waves are excited by other bodies such as moons in giant planets systems.

*Methods.* We used a local Cartesian analytical model, taking into account the complete Coriolis acceleration at any latitude, thus generalising previous works. We used a model in which stably stratified interfaces are infinitesimally thin, before relaxing this assumption with a second model that assumes a piecewise linear stratification.

*Results.* We find transmission of incident internal waves to be strongly affected by the presence of a density staircase, even if these waves are initially pure inertial waves (which are restored by the Coriolis acceleration). In particular, low-frequency waves of all wavelengths are perfectly transmitted near the critical latitude, defined by $\theta_c = \sin^{-1}(\omega/2\Omega)$, where $\omega$ is the wave's frequency and $\Omega$ is the rotation rate of the planet. Otherwise, short-wavelength waves are only efficiently transmitted if they are resonant with a free mode (interfacial gravity wave or short-wavelength inertial mode) of the staircase. In all other cases, waves are primarily reflected unless their wavelengths are longer than the vertical extent of the entire staircase (not just a single step).

*Conclusions.* We expect incident internal waves to be strongly affected by the presence of a density staircase in a frequency-, latitude- and wavelength-dependent manner. First, this could lead to new criteria to probe the interior of giant planets by seismology; and second, this may have important consequences for tidal dissipation and our understanding of the evolution of giant planet systems.

**Key words.** methods: analytical – planets and satellites: dynamical evolution and stability – planets and satellites: interiors – hydrodynamics – waves – planet-star interactions

## 1. Introduction

Since the first discovery of a planet orbiting a star outside our solar system (Mayor & Queloz 1995), astronomy has experienced an epoch of remarkable expansion: more than 3000 extrasolar planets are now confirmed, including more than 500 planets in multi-planetary systems<fnref n="1" />. Planet formation is thus a universal physical process. Planetary systems subsequently evolve dynamically by gravitational and magnetic interactions over astronomical timescales (for example Laskar et al. 2012; Bolmont & Mathis 2016; Strugarek 2016). The induced evolution of the orbital, rotational, thermal and compositional properties of the planets due to these interactions depends strongly on the internal structure of the planets involved (for example Ogilvie & Lin 2004; Efroimsky & Lainey 2007; Auclair-Desrotour et al. 2014). For example, the convective

instability, which is expected to operate in giant planet gaseous envelopes (which are the focus of attention here), is efficient at transporting heat, homogenises mean density profiles and mixes chemical elements. This may strongly impact tidal friction in these planets (Zahn 1966, 1989; Ogilvie & Lin 2004; Ogilvie & Lesur 2012; Mathis et al. 2016). However, whether the gaseous envelope is fully convective remains an open question.

Planetary interiors are poorly constrained. The JUNO spacecraft, orbiting Jupiter since July, 4, 2016, should provide high precision measurements of Jupiter's gravitational potential, aiming to constrain its interior (Militzer et al. 2016). Giant planet seismology, on the other hand, is very difficult because radial velocities associated with modes that can potentially be observed in Jupiter or Saturn are of very small amplitude (Gaulme et al. 2011). Based on Saturn's ring seismology, Fuller (2014) inferred that there could be a region of stable stratification in the deep

<footnote n="1">¹ Laboratoire AIM Paris-Saclay, CEA/DRF, CNRS, Univ. Paris-Diderot, IRFU/SAp Centre de Saclay, 91191 Gif-sur-Yvette, France
e-mail: quentin.andre@cea.fr</footnote>

<footnote n="2">² Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>

<footnote n="3">³ Département de Physique, ENS Paris-Saclay, Université Paris-Saclay, 61 Avenue du Président Wilson, 94230 Cachan, France</footnote>

<footnote n="4">⁴ Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds, LS2 9JT, UK
e-mail: a.j.barker@leeds.ac.uk</footnote>

<footnote n="5">⁵ LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris-Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France
e-mail: stephane.mathis@cea.fr</footnote>

<footnote n="1">¹ See http://exoplanet.eu</footnote>
interior of the giant planet, departing from the standard model of planetary interiors, which considers a large H/He convective gaseous envelope sitting on a rocky and/or icy core (that could be either fluid or solid, see for example Mazevet et al. 2015) expected from planet formation by core accretion (see for example Pollack et al. 1996).

Moreover, it has been shown that a stabilising compositional gradient could exist in certain regions of giant planet interiors as a natural outcome of planet formation and thermal evolution, thus competing with the destabilising entropy gradient that drives the convective instability. Namely, this is expected to occur in two different regions.

Firstly, just outside the core, the erosion of part of the core of giant planets was shown to be energetically plausible in Guillot et al. (2004). This is because in the conditions of temperature and pressure that reign in their central regions, some of the heavy elements composing the core (for example silicates) are thermodynamically unstable (for example Wilson & Militzer 2012b,a; Wahl et al. 2013; González-Cataldo et al. 2014; Mazevet et al. 2015). Thus, erosion and redistribution of core materials in the envelope must be taken into account, and could provide a stable compositional gradient. We also note that Stevenson (1985) suggested that impacts of planetesimals with a giant planet (this could happen before the proto-planetary disc clears) could lead to the formation of a stabilising gradient of heavy elements just outside the core.

Secondly, in the transition region between molecular and metallic H/He ices: the helium rain region (see Salpeter 1973; Stevenson 1975), the phase separation between H and He could provide a stabilising compositional gradient (Stevenson & Salpeter 1977) which could in turn trigger double diffusive convection (Nettelmann et al. 2015).

Furthermore, laboratory experiments have shown that the presence of a compositional gradient (i.e. a gradient of the mean molecular weight) can change the mean density profile that develops in a stratified fluid. This is due to the fact that, like temperature, the mean molecular weight influences the buoyancy of the fluid (Ledoux 1947). But because of diffusive processes, even a density stratification stable with respect to the convective instability can be unstable. This so-called double-diffusive instability, first theorised by Stern (1960), can arise if the diffusivity of one of the quantities (in general heat) is significantly greater than the other (heavy elements). In the case of oscillatory double-diffusive convection (also referred to as semi-convection, see Garaud 2013, for a review), the entropy gradient is destabilising while the compositional gradient is stabilising (the gradient of heavy elements is directed towards the planet centre). We note that without this stable chemical gradient, the envelope would be convectively unstable in the usual sense.

For particular parameter values, as the instability grows, the system can quickly develop a layered structure of well mixed, convective layers separated by thin stably stratified interfaces where both temperature and density undergo a sudden jump (Radko 2003). This so-called layered semi-convection is thus associated with a density staircase-like profile. This is confirmed by local three-dimensional non-linear numerical simulations (for example Rosenblum et al. 2011; Stellmach et al. 2011; Mirouh et al. 2012; Wood et al. 2013) and by direct observation on Earth, for example in the Canada basin in the Arctic Ocean (Ghaemsaidi et al. 2016). We note that in local numerical experiments, these layers are observed to merge (see for example Garaud 2013), but the long-term evolution of such a configuration is not currently understood.

Could this kind of layered structure exist in giant planet interiors and affect the dissipative processes at play? Several decades ago, Stevenson (1985) pointed out that such layered interior profiles could be relevant for solar system giant planets. Following that suggestion, Leconte & Chabrier (2012) proposed a giant planet interior model involving layered semi-convection, with the aim to verify whether such a model would be consistent with observational constraints. The new picture obtained depends from the standard picture of giant planets that assumes a three-layer structure composed of a rocky/icy core surrounded by a metallic H/He layer with a molecular H₂/He envelope on top. Indeed, because of its ability to hamper large-scale convection, we expect the presence of layered semi-convection in giant planet interiors to deeply modify the long-term interior evolution of planets. Saturn's infrared luminosity shows an excess compared to what is expected from the inherent gravitational contraction and cooling of the body, which cannot be explained invoking standard models of giant planet interiors (Pollack et al. 1977; Fortney et al. 2011). Leconte & Chabrier (2013) proposed that this could be explained invoking layered semi-convection. Before that, Chabrier & Baraffe (2007) showed that layered semi-convection could also play a role in explaining the abnormally large radius of some hot Jupiters (first noticed by Bodenheimer et al. 2001; and Guillot & Showman 2002; see also Baraffe et al. 2005), though its efficiency in practice has been questioned by Kurokawa & Inutsuka (2015). In any case, those findings seem to suggest that layered semi-convection could be a crucial ingredient in realistic models of giant planet internal structures.

If layered semi-convection and its associated density staircases are present in giant planet interiors, it is very important to determine ultimately how this complex structure affects the rates of tidal dissipation, which has not been studied before. In particular, it is crucial to determine whether layered semi-convection in giant planet interiors could account for the higher tidal dissipation than previously thought in Jupiter and Saturn found by Lainey et al. (2009, 2012, 2017) based on astrometric measurements spanning more than a century, including some from the Cassini spacecraft. In a series of papers, we will thus be driven by the following question: how would the presence of density staircases in giant planet interiors affect the propagation of internal waves and modify the rates of tidal dissipation?

We also note that semi-convection is thought to be able to produce a layered state in massive stars of mass $M_* \ge 15 M_\odot$, ($M_\odot$ denoting the mass of the Sun) outside their convective core, which contracts over time, potentially leaving a stabilising He gradient in the hydrogen envelope of the star (see Schwarzschild & Härm 1958; Sakashita & Hayashi 1959). As a consequence, the results of these papers, which are focused on giant planets, may also be relevant for those massive stars.

Before evaluating the rates of tidal dissipation in a layered profile (which will be the focus of attention of a second paper), we need to understand how density staircases affect the propagation and transmission of gravito-inertial waves that are potentially excited by tidal forcing. This paper aims to determine how density staircases modify the linear propagation of internal waves in a rotating planet. This is done by extending and generalising two previous studies: Belyaev et al. (2015; BQF15 hereafter), who derived the dispersion relation for the free modes of a staircase and considered the effects of rotation at the pole and equator; and Sutherland (2016; S16 hereafter), who studied the transmission of an incident internal wave upon a density staircase embedded in a stably stratified medium under the traditional approximation. The traditional approximation consists in
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Overview of our model. In the helium rain region or in the region just outside the core (where heavy elements, symbolised here by Z, could be released into the gaseous envelope), layered semi-convection could operate. The resulting density profile is staircase-like.</figcaption></figure>

neglecting the horizontal component of the rotation vector in the Coriolis acceleration, which is mostly valid in strongly stratified fluids (see for example Friedlander 1987). Based on a similar model of a plane-parallel density staircase, we study in detail the effects of rotation (by including the complete Coriolis force) at any latitude, to determine its effects on the free modes and transmission of incident waves. In Fig. 1, we give an overview of our reference physical model (and its regions of applicability discussed above), Convective layers of size *d*, in which density is uniform, are separated by infinitesimally thin stably stratified interfaces across which the density undergoes a discontinuous jump by a value $\Delta\rho$.

The outline of this paper is as follows. In Sect. 2, we present the important mathematical and physical aspects of the linear propagation of internal waves, using a formalism first introduced by Gerkema & Shrira (2005). Some energetical aspects are discussed in Sect. 2.4, to yield a general expression of the transmission coefficient. Section 3 then presents our study of the layered case. In particular, we generalise the dispersion relation obtained by BQF15 in Sect. 3.3, and derive a series of analytical expression for the transmission coefficient of an internal wave incident upon a density staircase in Sect. 3.4, with the aim to predict which waves will be able to penetrate into deeper regions of giant planets. A link is made between the bands of perfect transmission that arise and the free modes of the staircase given by the dispersion relation. In addition, we extend in Sect. 3.5 the physical model to a more realistic one for which stably stratified interfaces have a finite size. Finally, we summarise our main results and discuss their astrophysical implications, particularly for giant planets seismology, in Sect. 4.

## 2. Internal waves in giant planet interiors

### 2.1. Main assumptions

We wish to study the propagation of short-wavelength internal waves. Therefore we adopt a local Cartesian model (Gerkema & Shrira 2005; Mathis et al. 2014) that represents a small-patch of a giant planet (see the appendix of Ogilvie & Lin 2004; Auclair Desrotour et al. 2015, for the cases of pure inertial waves and gravito-inertial waves, respectively), simplifying the global spherical geometry. We centre our box on a point M of the gaseous envelope (see Fig. 2).

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Left: global view of a giant planet: the gaseous envelope (in yellow, the shading denoting density), lies on top of the core (in red). Right: magnified picture of the local Cartesian box, centred on a point M of a giant planet envelope, corresponding to a colatitude $\Theta$. The local box is tilted with respect to the spin axis, and its vertical axis $z$, corresponding to the local radial direction, is thus anti-aligned with gravity. The $x$ and $y$ axes correspond to the local azimuthal and latitudinal directions, respectively, while the $\chi$ axis makes an angle $\alpha$ with respect to the $x$-axis.</figcaption></figure>

The local system of coordinates $(x, y, z)$ corresponds to the local azimuthal, latitudinal and radial directions, respectively. The rotation vector $\Omega$ makes an angle $\Theta$ with respect to the gravity vector $g$, aligned with the vertical direction. Thus, in our local system of coordinates, the latitudinal and vertical components of the rotation vector are, respectively,

$$ \tilde{f} = 2\Omega \sin \Theta, \qquad (1) $$

$$ f = 2\Omega \cos \Theta, \qquad (2) $$

so that $2\Omega = (0, \tilde{f}, f)$.

Both the rotation rate of the planet $\Omega$ and the local gravity $g$ are assumed to be uniform and constant. We assume the rotation rate to be far below the breakup angular velocity ($\Omega_K = \sqrt{GM_p/R_p^3}$, where $G$ is the gravitational constant, $M_p$ and $R_p$ are the mass and radius of the planet), and accordingly we can ignore the centrifugal acceleration. We also introduce a reduced horizontal coordinate, $\chi$, that makes an angle $\alpha$ with respect to the $x$-axis: $\chi = x \cos \alpha + y \sin \alpha$.

One must keep in mind that the local approach is valid only for a box size of negligible extent compared to the characteristic length scale of the planet,

$$ L_x, L_y, L_z \ll R, \qquad (3) $$

otherwise curvature effects due to the spherical geometry should be taken into account. Here, $L_x$, $L_y$ and $L_z$ are the lengths of the box in the $x$, $y$ and $z$ directions, respectively, and $R$ is the radius of the planet. In addition, considering a constant gravity vector restricts us to consider dynamical phenomena with length scales $\lambda$ far below the pressure scale height $H_p$,

$$ \lambda \ll H_p. \qquad (4) $$

In giant planet deep interiors, we have $H_p \sim R$ (Leconte & Chabrier 2012), so that if condition (3) is fulfilled, condition (4) is as well, since $\lambda < L_x, L_y, L_z$. Finally, such
an approach is suitable because it is expected that tidally excited waves have a small-scale structure (Ogilvie & Lin 2004).

Our other main assumptions are the followings:

- We adopt the Boussinesq approximation. The fluid is assumed to be quasi-incompressible with a reference density value $\rho_0$, and accordingly we restrict our study to low Mach numbers, i.e. $|u| \ll c_s$, where $|u|$ is the velocity and $c_s$ is the sound speed. In addition, the vertical extent occupied by the fluid is far below the pressure scale height: $L_z \ll H_p$, a condition that is fulfilled because we use a local approach.

- Dissipative processes (viscosity and thermal diffusion) are not taken into account, but they will be in our second paper, which will be focused on tidal dissipation.

- The spatial dependence of the background quantities is assumed to be fixed, as resulting from double-diffusive instabilities as described by for example Leconte & Chabrier (2012). Thus, the back-reaction of internal waves on the layered structure is not taken into account, nor is the possible excitation of internal waves by double-diffusive convection (see for example Moll et al. 2016).

- Non-linear effects are entirely neglected.

## 2.2. Equations of motion

Before studying the propagation (and transmission) of internal waves in density staircases associated with layered semi-convection, we need to introduce the formalism that allows us to treat gravito-inertial waves (GIWs) with the complete Coriolis acceleration, as well as compute the corresponding energetic quantities. In Sect. 3, which is the heart of this paper, we will use these results.

We study the linear propagation of GIWs in the local Cartesian model. First, let us introduce the velocity field,

$$u(\mathbf{r}, t) = \begin{pmatrix} u(\mathbf{r}, t) \\ v(\mathbf{r}, t) \\ w(\mathbf{r}, t) \end{pmatrix}, \qquad (5)$$

where $u$, $v$ and $w$ are the components of the velocity perturbation in the local azimuthal, latitudinal and radial directions, respectively. Next, we define the fluid buoyancy,

$$b(\mathbf{r}, t) = -g \frac{\rho(\mathbf{r}, t)}{\rho_0}, \qquad (6)$$

where $\rho(\mathbf{r}, t)$ is the density fluctuation field. Then, we define the buoyancy frequency in the Boussinesq approximation,

$$N^2(z) = - \frac{g}{\rho_0} \frac{d\bar{\rho}}{dz}, \qquad (7)$$

where $\bar{\rho}(z)$ is the resulting background density profile. We stress that $\rho_0$ is the reference Boussinesq density value, while $\bar{\rho}$ is the background density distribution associated with the layered density profile. We assume $\max_{|z|<L_z/2} \bar{\rho}(z) - \min_{|z|<L_z/2} \bar{\rho}(z) \ll \rho_0$. In Sect. 3,

we will model layered semi-convection by a succession of convective layers, in which we assume $N^2 = 0$, separated by infinitesimally thin stably stratified interfaces, in which $N^2 > 0$.

The three linearised components of the momentum equation are given by

$$\frac{\partial u}{\partial t} - f v + \tilde{f} w = -\frac{1}{\rho_0} \frac{\partial p}{\partial x}, \qquad (8)$$

$$\frac{\partial v}{\partial t} + f u = -\frac{1}{\rho_0} \frac{\partial p}{\partial y}, \qquad (9)$$

$$\frac{\partial w}{\partial t} - \tilde{f} u = -\frac{1}{\rho_0} \frac{\partial p}{\partial z} + b, \qquad (10)$$

where $p(\mathbf{r}, t)$ is the pressure fluctuation, and $\tilde{f}$ and $f$ are expressed by Eqs. (1) and (2), respectively. Then, we write the continuity equation,

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0. \qquad (11)$$

Finally, we write the thermal energy equation in the adiabatic limit,

$$\frac{\partial b}{\partial t} + N^2 w = 0. \qquad (12)$$

## 2.3. Propagation of gravito-inertial waves

### 2.3.1. Dispersion relation in a uniformly stratified medium

We now introduce some key aspects to analyse the linear propagation of GIWs. We consider monochromatic plane wave solutions of the form

$$x(\mathbf{r}, t) = \Re \{X \exp [i(\mathbf{k} \cdot \mathbf{r} - \omega t)]\}, \qquad (13)$$

where $\omega$ is the wave frequency, and $\mathbf{k}$ the wave vector. Here, $x$ stands for either $\rho$, $p$, $b$, $u$, $v$ or $w$.

Substituting the above solution into Eqs. (8)–(12), we obtain the dispersion relation for GIWs in a uniformly stratified ($N(z) = N_0$) and rotating medium:

$$\omega^2 = N_0^2 \frac{k_z^2}{k^2} + \frac{(2\Omega \cdot \mathbf{k})^2}{k^2}. \qquad (14)$$

One can then derive the frequency domain for which GIWs propagate, corresponding to a real frequency $\omega$. This will be done in Sect. 2.3.4.

In a convective region, in which we assume isentropy ($N_0 = 0$), we recover the dispersion relation for pure inertial waves (IWs) given by

$$\omega^2 = \frac{(2\Omega \cdot \mathbf{k})^2}{k^2}. \qquad (15)$$

Similarly, when rotation is absent, we recover the dispersion relation for pure internal gravity waves given by

$$\omega^2 = N_0^2 \frac{k_z^2}{k^2}. \qquad (16)$$

### 2.3.2. Group velocities

From the dispersion relation given by Eq. (14), we can obtain that the group velocity of GIWs is

$$v_g^{(\text{GIW})} = \frac{1}{\omega} \left[ N_0^2 \left( \frac{k_z}{k} \right) \frac{\mathbf{k} \times (-\hat{\mathbf{e}}_z \times \mathbf{k})}{k^3} + \left( \frac{2\Omega \cdot \mathbf{k}}{k} \right) \frac{\mathbf{k} \times (2\Omega \times \mathbf{k})}{k^3} \right]. \qquad (17)$$

Thus, the energy carried by GIWs can propagate along two directions, corresponding to the two signs of $\omega$.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Propagation of energy carried by pure inertial waves at the frequency $\omega = f$. The two allowed group velocity vectors, $v_{g,\pm}$, form an angle $\lambda_{\pm} = \pm(\pi/2 - \Theta)$ with respect to the rotation axis, represented in red. Thus, the energy propagates along two characteristics, represented in blue. At the frequency $\omega = f$, one of the characteristics is perpendicular to the vertical axis, so that the energy propagates in the local horizontal direction.</figcaption></figure>

The case of pure inertial waves. Let us focus more precisely on the case of pure inertial waves, which propagate in convective regions. Setting $N_0 = 0$ in Eq. (17) and using Eq. (15), we get that the group velocity of pure inertial waves is given by

$$v_g = \operatorname{sgn}(\omega) \frac{k \times (2\Omega \times k)}{k^3}. \qquad (18)$$

It is worth considering along which direction the energy of a pure inertial wave propagates, because in a region of layered semi-convection, we expect the volume to be mostly convective (the interfaces are very thin). This can be done with a little algebra, which yields that the group velocity of pure inertial waves makes an angle

$$\lambda = \pm \sin^{-1} \left( \frac{\omega}{2\Omega} \right) \qquad (19)$$

with respect to the rotation vector.

At the frequency $\omega = f = 2\Omega \cos \Theta$, we have

$$\lambda_{\pm} = \pm \left(\frac{\pi}{2} - \Theta\right). \qquad (20)$$

The solutions are displayed in Fig. 3. At the frequency $\omega = f$, one of the energy propagation directions is perpendicular to the vertical axis, so that the energy propagates along the local horizontal and thus does not propagate towards deeper regions of the giant planet.

Recalling the simplified physical model to be adopted (shown in Fig. 1), this means that at this particular frequency, one of the directions of energy propagation is parallel to the stably stratified interfaces (which lie in the local horizontal plane).

### 2.3.3. Poincaré equation

Now, we introduce a key partial differential equation (PDE) to study the dynamics of GIWs. By reducing the system (8)–(12), we can derive a PDE solely for the vertical velocity $w$ (see Appendix A),

$$\frac{\partial^2}{\partial t^2} \nabla^2 w + (f \cdot \nabla)^2 w + N^2 \nabla_{\perp}^2 w = 0, \qquad (21)$$

where $\nabla_{\perp}^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$ is the horizontal Laplacian, and we recall that $f = (0, \tilde{f}, f)$.

We want to study the propagation of a given monochromatic GIW with a frequency $\omega$, that propagates in the direction $(\cos \alpha, \sin \alpha)$ in the $(Mxy)$ plane (see Fig. 2). Substituting $w = W(\chi, z) \exp[i\omega t]$ where $\chi = x \cos \alpha + y \sin \alpha$ is the reduced horizontal coordinate, we obtain the Poincaré equation for GIWs:

$$[N^2(z) - \omega^2 + \tilde{f}_s^2] \frac{\partial^2 W}{\partial \chi^2} + 2f \tilde{f}_s \frac{\partial^2 W}{\partial \chi \partial z} + [f^2 - \omega^2] \frac{\partial^2 W}{\partial z^2} = 0, \quad (22)$$

where

$$\tilde{f}_s = \tilde{f} \sin \alpha. \qquad (23)$$

The case $\omega = f$. Let us focus first on the case $\omega = f$. Substituting

$$W(\chi, z) = \tilde{W}(z) \exp[ik_{\perp}\chi], \qquad (24)$$

where $k_{\perp}$ the wave number in the $\chi$ direction, we obtain

$$2f \tilde{f}_s \frac{d\tilde{W}}{dz} + ik_{\perp} [N^2 - f^2 + \tilde{f}_s^2] \tilde{W} = 0. \qquad (25)$$

This equation being of first order in $z$, it has only one solution, which we can easily calculate in the case of a uniformly stratified medium ($N = N_0$):

$$\tilde{W}(z) \propto \exp \left[ i \left( \frac{N_0^2 - f^2 + \tilde{f}_s^2}{2f \tilde{f}_s} \right) k_{\perp} z \right]. \qquad (26)$$

The case $\omega \neq f$. For $\omega \neq f$, following Gerkema & Shrira (2005) we introduce the transformation

$$w = \hat{W}(z) \exp \left[ i(k_{\perp}(\chi + \delta z) - \omega t) \right], \qquad (27)$$

where

$$\delta = \frac{f \tilde{f}_s}{\omega^2 - f^2}. \qquad (28)$$

Substituting Eq. (27) into (21) leads to

$$\frac{d^2 \hat{W}}{dz^2} + k_z^2 \hat{W} = 0, \qquad (29)$$

where

$$k_z^2 = k_{\perp}^2 \left[ \frac{N^2 - \omega^2}{\omega^2 - f^2} + \left( \frac{\omega \tilde{f}_s}{\omega^2 - f^2} \right)^2 \right]. \qquad (30)$$

In a uniformly stratified or in a convective medium (which both have $N(z) = \text{const}$), solutions of Eq. (29) have the form $\hat{W}(z) \propto \exp(\pm ik_z z)$. We have $k_z^2 > 0$ in the propagative regime and $k_z^2 < 0$ in the evanescent regime.

We stress that the transformation given by Eq. (27) has the effect of splitting the vertical wave number into two parts: $\pm k_z$, contained in the $z$-dependence of the $\hat{W}$ function, and $k_{\perp} \delta$, contained in the exponential factor of Eq. (27). Thus, we define the total vertical wave number as

$$\tilde{k}_z = \pm k_z + k_{\perp} \delta, \qquad (31)$$

so that $w \propto \exp[i(k_{\perp}\chi + \tilde{k}_z z - \omega t)]$. The non-traditional component $k_{\perp} \delta$ corresponds to the intrinsic 2D behaviour of GIWs when taking the complete Coriolis acceleration into account. It vanishes at the pole or for a null rotation, for which the problem is separable.
### 2.3.4. Frequency spectrum

We now describe GIWs propagation as a function of their frequency $\omega$, along the lines of Mathis et al. (2014). At a given $z$, GIWs are propagative if $k_z^2 > 0$, where $k_z^2$ is defined by Eq. (30), which occurs when

$$ \omega_{-} < \omega < \omega_{+}, \qquad (32) $$

with

$$ \omega_{\pm} = \frac{1}{\sqrt{2}} \sqrt{[N^2 + 4\tilde{\Omega}^2] \pm \sqrt{[N^2 + 4\tilde{\Omega}^2]^2 - (2fN)^2}}, \quad (33) $$

where we have defined a modified rotation rate of the planet

$$ \tilde{\Omega} = \frac{1}{2} \sqrt{f^2 + \tilde{f}_s^2} = \Omega \sqrt{1 - \sin^2 \Theta \cos^2 \alpha}. \quad (34) $$

In convective regions, the local vertical wave number given by Eq. (30) is obtained by setting $N = 0$ and becomes

$$ k_{z,c}^2 = k_{\perp}^2 \frac{\omega^2(4\tilde{\Omega}^2 - \omega^2)}{(\omega^2 - f^2)^2}. \quad (35) $$

In Fig. 4, we illustrate the spectrum of internal waves, depending on their frequency $\omega$ and on the type of layer they propagate in: convective or stably stratified. In a stably stratified layer, in which both rotation and stratification are present, gravito-inertial waves propagate for $\omega_{-} < \omega < \omega_{+}$, as explained above. In a convective layer, we have pure inertial waves propagative for $\omega < 2\tilde{\Omega}$, and evanescent for $\omega > 2\tilde{\Omega}$.

### 2.3.5. Properties of the reflected wave upon an interface

Goodman & Lackner (2009) have considered the reflection of a monochromatic plane pure inertial wave upon a rigid wall, with normal $\hat{e}_z$. In this section, we extend their work to the more general case of a gravito-inertial wave, which obeys the dispersion relation given by Eq. (14). The wave vector is expanded as $\mathbf{k} = k_{\perp} \hat{e}_{\chi} + \tilde{k}'_{z} \hat{e}_{z}$.

The goal of the following calculation is to understand how the vertical wave vector of the reflected wave is expected to change as a result of the non-specular reflection from an interface. In our case, the interfaces are in the horizontal plane, with $\hat{e}_z$ being the vector normal to the interfaces. We thus consider incident and reflected waves of the form

$$ \hat{W}_{\text{in}} \propto \exp [i(\mathbf{k} \cdot \mathbf{r} - \omega t)], \quad (36) $$

$$ \hat{W}_{\text{re}} \propto \exp [i(\mathbf{k}' \cdot \mathbf{r} - \omega't)], \quad (37) $$

where the scattered (outgoing) wave vector $\mathbf{k}' = \mathbf{k}'_{\perp} \hat{e}_{\chi} + \tilde{\mathbf{k}}'_{z} \hat{e}_{z}$ is determined by two conditions. First, in order that the incident and reflected waves have the same relative phase at all points along the interface, as required by the impermeability condition, it is necessary that $\mathbf{k}$ and $\mathbf{k}'$ have the same components parallel to the interface. That means that $\hat{e}_z \times (\mathbf{k}' - \mathbf{k}) = \mathbf{0}$, or equivalently $\mathbf{k}'_{\perp} = \mathbf{k}_{\perp}$. Thus,

$$ \mathbf{k}' = k_{\perp} \hat{e}_{\chi} + \tilde{k}'_{z} \hat{e}_{z}, \quad (38) $$

where $\tilde{k}'_z$ (the total vertical wave number) remains to be determined. Similarly, in order that the relative phase stays constant in time, the two waves must have the same frequency, i.e $\omega' = \omega$, so that (see Eq. (14))

$$ \omega^2 = \frac{k_{\perp}^2}{k'^2} N^2 + \frac{(2\Omega \cdot \mathbf{k}')^2}{k'^2}. \quad (39) $$

Substituting (38) into (39) leads to

$$ (\omega^2 - f^2) \tilde{k}_z^2 - 2f \tilde{f}_s k_{\perp} \tilde{k}'_z + [\omega^2 - (\tilde{f}_s^2 + N^2)] k_{\perp}^2 = 0, \quad (40) $$

where we have used that $(2\Omega \cdot \mathbf{k}) = \tilde{f}_s k_{\perp} + f \tilde{k}_z$. From Eq. (40), it is possible to get the two roots, which are

$$ \tilde{k}'_z = \left( \frac{f \tilde{f}_s}{\omega^2 - f^2} \right) k_{\perp} \pm \frac{\sqrt{(f \tilde{f}_s)^2 - (\omega^2 - f^2)[\omega^2 - (\tilde{f}_s^2 + N^2)]}}{\omega^2 - f^2} k_{\perp} \\ = \tilde{\delta} k_{\perp} \pm k_z, \quad (41) $$

where $\tilde{\delta}$ and $k_z$ are given by Eqs. (28) and (30), respectively. We recognise that the expression of $\tilde{k}'_z$ given above matches the definition of the total vertical wave number defined by Eq. (31), the plus sign corresponding to the incident (ingoing) wave, and the minus sign corresponding to the reflected (outgoing) wave. Therefore, this calculation provides an independent check that the formalism introduced with Eq. (27), which splits the vertical wave number into two parts, is correct.

### 2.3.6. Polarization relations

The Poincaré equation given by Eq. (21) is a PDE solely for the vertical velocity $w$. In order to study the dynamics of the other fields: the horizontal components of the velocity, $u$ and $v$, the pressure $p$, and the buoyancy $b$, we can derive analytic formulae to get $u, v, p$ and $b$ in term of $w$ and its derivatives. To do that, we first express each field $x(\chi, z, t)$ describing the perturbed flow using the transformation introduced in Eq. (27) for the vertical velocity,

$$ x = \Re \{X(\chi, z)e^{-i\omega t}\}, \quad (42) $$

where

$$ X = \hat{X}(z) \exp [i k_{\perp} (\chi + \tilde{\delta} z)]. \quad (43) $$

We stress again that $\hat{X}$ carries only part of the vertical dependence, the other part being included in $\exp(ik_{\perp}\tilde{\delta}z)$. The system of equations in Sect. 2.2 thus becomes

$$ -i\omega U - fV + \tilde{f}W = -\frac{ik_{\perp} \cos \alpha}{\rho_0} P, \quad (44) $$

$$ -i\omega V + fU = -\frac{ik_{\perp} \sin \alpha}{\rho_0} P, \quad (45) $$

$$ -i\omega W - \tilde{f}U = -\frac{1}{\rho_0} \frac{\partial P}{\partial z} + B, \quad (46) $$

for the three components of the equation of momentum, and

$$ ik_{\perp} \cos \alpha U + ik_{\perp} \sin \alpha V + \frac{\partial W}{\partial z} = 0, \quad (47) $$

$$ -i\omega B + N^2 W = 0, \quad (48) $$

for the equations of conservation of mass and energy, respectively. Noting that

$$ \frac{\partial W}{\partial z} = (\hat{W}' + ik_{\perp} \tilde{\delta} \hat{W}) \exp [ik_{\perp}(\chi + \tilde{\delta} z)], $$
<figure><img src="image_7.png" /><figcaption>Fig. 4. Low-frequency spectrum for internal waves in a rotating giant planet. Waves in the convective and stably stratified regions are indicated at the top and bottom, respectively. The purple box corresponds to sub-inertial gravito-inertial waves that are propagative in both convective and stably stratified regions. Adapted from Mathis et al. (2014).</figcaption></figure>

all the fields can be expressed in term of $\hat{W}$ as follows:

$$ \hat{U} = \frac{\tilde{f}_s(f \cos \alpha - i\omega \sin \alpha)}{f^2 - \omega^2} \hat{W} + \frac{f \sin \alpha + i\omega \cos \alpha}{\omega k_{\perp}} \hat{W}', \quad (49) $$

$$ \hat{V} = \frac{\tilde{f}_s(f \sin \alpha + i\omega \cos \alpha)}{f^2 - \omega^2} \hat{W} - \frac{f \cos \alpha - i\omega \sin \alpha}{\omega k_{\perp}} \hat{W}', \quad (50) $$

$$ \hat{P} = i \frac{\rho_0 \tilde{f}_c}{k_{\perp}} \hat{W} + i \frac{\rho_0 (f^2 - \omega^2)}{\omega k_{\perp}^2} \hat{W}', \quad (51) $$

$$ \hat{B} = i \frac{N^2}{\omega} \hat{W}, \quad (52) $$

where the prime denotes differentiation with respect to z, and $\tilde{f}_c = f \cos \alpha$. Then, each field has to be multiplied by $\exp[i(k_{\perp}(\chi + \delta z) - \omega t)]$ to get the complete solution. The physical solution is then the real part of the complete complex solution.

## 2.4. Energetical aspects

One of our motivations is to predict energy transport by internal waves in regions of layered semi-convection. In this section, we thus focus on the energetics of GIWs propagation.

### 2.4.1. Expression of the kinetic and potential energies in term of $\hat{W}$

Kinetic and potential energy densities, averaged over a wave period $2\pi/\omega$, can be written as

$$ \begin{aligned} E_k &= \frac{1}{2} \rho_0 \left\langle \mathfrak{R}(u)^2 + \mathfrak{R}(v)^2 + \mathfrak{R}(w)^2 \right\rangle \\ &= \frac{1}{4} \rho_0 (UU^* + VV^* + WW^*), \end{aligned} \quad (53) $$

and

$$ \begin{aligned} E_p &= \frac{1}{2} \rho_0 \frac{\langle \mathfrak{R}(b)^2 \rangle}{N^2} \\ &= \frac{1}{4} \rho_0 \frac{BB^*}{N^2}, \end{aligned} \quad (54) $$

where the brackets $\langle \cdot \rangle$ denotes time averaging, and the asterisks, the complex conjugate. To get Eqs. (53) and (54), we have used the identity $\mathfrak{R}(A)\mathfrak{R}(B) = \frac{1}{2}\mathfrak{R}(AB + AB^*)$, so that provided the transformation given by Eq. (42), $\langle \mathfrak{R}(x)^2 \rangle = \frac{1}{2}\langle \mathfrak{R}(X^2e^{-2i\omega t} + XX^*) \rangle = \frac{1}{2}XX^*$. Then, using the transformation given by Eq. (43), we have $XX^* = \hat{X}\hat{X}^*$, so that the kinetic and potential energies can be written as

$$ E_k = \frac{1}{4} \rho_0 (\hat{U}\hat{U}^* + \hat{V}\hat{V}^* + \hat{W}\hat{W}^*) \quad (55) $$

and

$$ E_p = \frac{1}{4} \rho_0 \frac{\hat{B}\hat{B}^*}{N^2}. \quad (56) $$

Then, using the polarization relations given by Eqs. (49), (50) and (52), we can express the kinetic and potential energies in term of $\hat{W}$ only as

$$ E_k = \frac{1}{4} \rho_0 \left\{ \left[ \frac{\tilde{f}_s^2 (f^2 + \omega^2)}{(\omega^2 - f^2)^2} + 1 \right] \hat{W} \hat{W}^* + \frac{1}{k_{\perp}^2} \left[ \left( \frac{f}{\omega} \right)^2 + 1 \right] \hat{W}' (\hat{W}')^* + \frac{4 f \tilde{f}_s}{k_{\perp} (f^2 - \omega^2)} \mathfrak{F} [\hat{W} (\hat{W}')^*] \right\} \quad (57) $$

and

$$ E_p = \frac{1}{4} \rho_0 \left(\frac{N}{\omega}\right)^2 \hat{W} \hat{W}^*, \quad (58) $$

where we recall that $\hat{W}$ is a function of z only. Assuming a plane wave form in the vertical direction also (as suggested by Eq. (29)),

$$ \hat{W} = \mathcal{A} \exp(ik_z z), \quad (59) $$

we find that the total energy is

$$ E_k + E_p = \frac{1}{4} \rho_0 \left\{ \left[ \frac{\tilde{f}_s^2 (f^2 + \omega^2)}{(\omega^2 - f^2)^2} + \left(\frac{N}{\omega}\right)^2 + 1 \right] + \left[ \left(\frac{f}{\omega}\right)^2 + 1 \right] \left(\frac{k_z}{k_{\perp}}\right)^2 + \frac{4 f \tilde{f}_s}{(\omega^2 - f^2) k_{\perp}} \frac{k_z}{k_{\perp}} \right\} |\mathcal{A}|^2. \quad (60) $$

### 2.4.2. Expression of the vertical energy flux density in term of $\hat{W}$

In our setup, the energy flux density is $\langle \rho u \rangle$ in real variables, so that the vertical energy flux density is given by

$$ F_z = \langle \mathfrak{R}(w)\mathfrak{R}(p) \rangle = \frac{1}{2} \mathfrak{R}(\hat{W}\hat{P}^*). \quad (61) $$

Using Eq. (51) for $\hat{P}$, we get

$$ F_z = \frac{1}{2} \rho_0 \left( \frac{f^2 - \omega^2}{\omega k_{\perp}^2} \right) \mathfrak{F} [\hat{W} (\hat{W}')^*]. \quad (62) $$

Finally, using Eq. (59) for $\hat{W}$, we obtain that the vertical energy flux density is given by

$$ F_z = \frac{1}{2} \rho_0 \left( \frac{\omega^2 - f^2}{\omega k_{\perp}^2} \right) k_z |\mathcal{A}|^2, \quad (63) $$
<figure><img src="image_8.png" /><figcaption>Fig. 5. Illustration of the system to be considered in Sect. 3, consisting of a density staircase (S) of vertical extent D. An incident wave of amplitude A<sub>in</sub>, thus carrying an energy density ∝|A<sub>in</sub>|<sup>2</sup>, comes from above (S). The reflected wave has an amplitude A<sub>re</sub> and thus carries an energy density ∝|A<sub>re</sub>|<sup>2</sup> = R|A<sub>in</sub>|<sup>2</sup>, while the transmitted wave has an amplitude A<sub>tr</sub> and thus carries an energy density ∝|A<sub>tr</sub>|<sup>2</sup> = T|A<sub>in</sub>|<sup>2</sup>. R and T are the reflection and transmission coefficients, respectively. The arrows indicate the schematic vertical direction of propagation of energy. The horizontal components are not shown.</figcaption></figure>

where we recall that A is the amplitude of the vertical compo-
nent of the velocity, and k_z the vertical wave number.

It is interesting to note that the vertical energy flux density
can be expressed in an alternative manner, noting that the energy
is transported vertically with the vertical component of the group
velocity, v_g · ˆe_z. That suggests that the vertical energy flux density
can be expressed as

$$
F_z = (v_g \cdot \hat{e}_z) (E_k + E_p), \quad (64)
$$

where the group velocity is defined by Eq. (15) and the sum of
the mean kinetic and potential energies has been derived in the
case of a plane wave in Eq. (60). We stress that this expression
only makes sense when E<sub>k</sub> and E<sub>p</sub> are mean quantities, aver-
aged over a wavelength and period, as it has been done to obtain
Eq. (60). It has been checked that both expressions are equal,
providing that k<sub>z</sub> takes the expression given by Eq. (30).

2.4.3. Expression of the transmission coefficient

In Sect. 3, we will study in detail the transmission of inter-
nal waves through a portion of a density staircase produced
by semi-convection: an incident (downward propagating) wave
carrying an energy density ∝|A<sub>in</sub>|<sup>2</sup> enters the staircase (S), of
vertical extent D. At the top of the staircase, a reflected (up-
ward propagating) wave is created, carrying an energy density
∝ |A<sub>re</sub>|<sup>2</sup> = R|A<sub>in</sub>|<sup>2</sup>. At the bottom of the staircase, the transmit-
ted (downward propagating) wave comes out, carrying an en-
ergy density ∝ |A<sub>tr</sub>|<sup>2</sup> = T|A<sub>in</sub>|<sup>2</sup>. Here, R and T are the reflection
and transmission coefficients, respectively; this is illustrated in
Fig. 5.

The transmission coefficient itself is defined to be the ratio of
transmitted (labelled by (tr)) to incident (labelled by (in)) energy
flux densities, for which we have found an analytic expression in
the previous section:

$$
T = \frac{F_z^{(\text{tr})}}{F_z^{(\text{in})}}. \qquad (65)
$$

In our case, using Eq. (63) for the vertical energy flux density
gives

$$
T = \frac{k_z^{(\text{tr})}}{k_z^{(\text{in})}} \left| \frac{\mathcal{A}_{\text{tr}}}{\mathcal{A}_{\text{in}}} \right|^2 . \qquad (66)
$$

This is the expression we will use to calculate transmission coef-
ficients. The regions above and below (S), respectively defined
by z > 0 and z < −D (see Fig. 5), are a priori different: in par-
ticular, the stratification (stable or unstable) in those regions is
not necessarily the same, so that the vertical wave numbers k_z
will not be either. We note that here, k_z is not the total vertical
wavenumber, that we denote by ˜k_z (see Eq. (31)), but only the
part defined through Eqs. (29)–(30).

In the absence of energy sources such as background shear,
we expect that the energy is conserved. Therefore, the reflection
coefficient is given by

$$
R = 1 - T. \tag{67}
$$

3. Propagation of internal waves in layered semi-convection

Armed with the results of the previous section, we can now study
the propagation of internal waves in a idealised model describing
layered semi-convection.

3.1. Physical set up

We continue to work in the local Cartesian box described in
Sect. 2.1. The background state is assumed to be in hydrostatic
equilibrium with constant gravity pointing in the −ê<sub>z</sub> direction,
so that dP/dz = −ρ<sub>0</sub>g.

We consider a region of a giant planet envelope in which double-diffusive convection has produced a layered density profile, as described in the introduction. Thus, the density profile is close to a density staircase in which convective layers of size *d* are separated by infinitesimally thin stably stratified interfaces. This idealised reference model is similar to that considered by BQF15 and S16. The distance between adjacent interfaces, *d*, is assumed constant. We will later discuss the effects of relaxing some of these assumptions in Sects. 3.4.5 and 3.5. At each interface, the density undergoes a discontinuous jump by a value Δρ > 0. This is illustrated in Fig. 6, on which are displayed the physical quantities introduced so far.

The Boussinesq approximation restricts us to have *d* ≪ *H*, where *H* is the characteristic length over which the background quantities vary. Under this condition, the magnitude of the density jump, Δρ, between adjacent steps (which we obtain by integrating Eq. (7) over one step including the interface) is given by

$$
\bar{N}^2 = \frac{g\Delta\rho}{\rho_0 d} \qquad (68)
$$

(BQF15). Here, N̄ corresponds to a mean stratification: typi-
cally, an internal wave with a large wavelength compared to the
size of the steps will see the staircase as a continuously strati-
fied medium characterised by the buoyancy frequency N̄. Since
N = 0 within a convective layer, the only contribution to N̄
comes from the density jump across the discrete interfaces. Also,
because d ≪ H by assumption, we generally have Δρ ≪ ρ₀, so
that the background density profile ρ̄(z) does not depart signifi-
cantly from the constant reference density value ρ₀.
<figure><img src="image_9.png" /><figcaption>Fig. 6. Model of the internal structure of a giant planet hosting layered semi-convection. The red area represents the core, while the white one represents the gaseous envelope. Within the latter, double-diffusive convection acts to create a staircase-like profile for the density (in orange). Convective layers of vertical extent *d* are separated by infinitesimally thin stably stratified interfaces. The density undergoes a density jump by a value Δρ at each interface. The rotation axis $\Omega$ (in red) forms an angle $\Theta$ with respect to the gravity *g* (in blue), directed along the (Mz) direction.</figcaption></figure>

In Fig. 7, three panels are displayed to summarise our reference model. Panel 7a shows the general scheme: *m* convective steps, each of vertical extent *d*, are separated by infinitesimally thin stably stratified interfaces. The incident (ingoing) wave – whose vertical velocity has an amplitude $\mathcal{A}_0$ – enters the staircase from above (region $z > 0$); a reflected wave with amplitude $\mathcal{B}_0$ is created, and a transmitted (outgoing) wave with amplitude $\mathcal{A}_{m+1}$ comes out of the staircase (region $z < -md$). By causality, there is no upward propagating wave below the staircase (S16). The regions above and below the staircase can either be convective ($N = 0$), or stably stratified ($N > 0$). We will denote the stratification above and below the staircase by $N_a$ and $N_b$, respectively. We note that Fig. 7 corresponds to the particular case where $N_a = N_b = 0$, i.e. a staircase embedded in a convective medium, but we also consider the general case.

In each convective step, labelled by the integer *n* ranging from 0 to *m*, the solution is the sum of an upward and a downward propagating (or evanescent) wave, whose vertical velocities have amplitudes $\mathcal{A}_n$ and $\mathcal{B}_n$, respectively.

We recall that the propagative or evanescent behaviour of a pure inertial wave depends on the value of $|\omega/2\Omega| (<1$ for propagative and $>1$ for evanescent). Panel 7b shows the corresponding density profile: the density undergoes a discontinuous jump by a value Δρ across each interface, and in between is uniform. This creates a mean density gradient

$$ \left| \frac{d\bar{\rho}}{dz} \right| = \frac{\Delta\rho}{d}. \qquad (69) $$

Finally, the buoyancy frequency profile $N^2(z)$ is displayed on panel 7c. It consists of a sum of Dirac distributions centred on the interfaces, each with an integrated value of $\bar{N}^2$, such as to create a profile

$$ N^2(z) = \begin{cases} N_a^2 & \text{for } z > 0, \\ \bar{N}^2 \sum_{n=0}^{m} \delta(z + nd) & \text{for } 0 < z > -md, \\ N_b^2 & \text{for } z < -md, \end{cases} \qquad (70) $$

where $\delta$ denotes the Dirac distribution. We stress that the vertical wave number, $k_z$, is a function of the buoyancy frequency $N$. This means that the incident gravito-inertial wave (or pure inertial wave if $N_a = 0$) has a vertical wave number $k_{z,a}$ which is in general different from the one of the transmitted wave, $k_{z,b}$, which are both a priori different from the wave number of purely inertial waves inside the convective steps, $k_{z,c}$ ($N_c = 0$). Remembering Eq. (30), we have defined

$$ k_{z,\alpha}^2 = -k_{\perp}^2 \left[ \frac{N_{\alpha}^2 - \omega^2}{\omega^2 - f^2} + \left( \frac{\omega \tilde{f}_s}{\omega^2 - f^2} \right)^2 \right]. \qquad (71) $$

Here, $\alpha$ stands for either a, b or c. We note that this definition is minus the one given by Eq. (30), a choice we will explain in the following section. Once again, we stress that the total vertical wave number is given by $\tilde{k}_z = \pm k_z + \delta k_{\perp}$.

In the following sections, we will mostly consider a finite staircase constituted of *m* steps, as described above, but we will also consider as a first approach in Sect. 3.3.1 an infinite staircase, for which $m \gg 1$.

## 3.2. Mathematical statement of the problem

The quantity we will manipulate is the GW's vertical velocity $w$, whose dynamics is governed by the Poincaré equation (see Sect. 2.3). Above and below the staircase (where $N = N_a, N_b$ respectively) and within each step (where $N = 0$), the differential equation given by Eq. (25) for $\omega = f$ and by Eq. (29) otherwise (modified by the definition of $k_z$ given by Eq. (71)), can be solved explicitly.

The case $\omega = f$. First, let us focus on the case $\omega = f$. In that case, the solution of the equation governing the vertical dependence of the vertical velocity is given by Eq. (25). The solution for prescribed $k_{\perp}$ and $\omega$ is

$$ \tilde{W}(z) = \begin{cases} \mathcal{E}_0 e^{i\gamma_a z} & \text{for } z > 0, \\ \mathcal{E}_n e^{i\gamma(z+nd)} & \text{for } -nd < z < -(n+1)d, n = \{1, \dots, m\} \\ \mathcal{E}_{m+1} e^{i\gamma_b(z+md)} & \text{for } z < -md, \end{cases} \qquad (72) $$

where we have defined

$$ \gamma_{\alpha} = \left( \frac{N_{\alpha}^{2} - f^{2} + \tilde{f}_{s}^{2}}{2f \tilde{f}_{s}} \right) k_{\perp} \qquad (73) $$

(we have omitted the subscript “c” in the convective steps to make notations lighter). Equation (25) being of first order in $z$, the only boundary condition we need to get the full solution comes from the requirement that the two fluids above and below the interface stay in contact: $w_n = w_{n+1}$, for $n = \{0, \dots, m+1\}$. It is then straightforward that the amplitudes of the incident and transmitted waves, $\mathcal{E}_0$ and $\mathcal{E}_{m+1}$ respectively, are related by

$$ \mathcal{E}_0 = e^{i(m-1)\gamma_d} \mathcal{E}_{m+1}. \qquad (74) $$

The coefficient $\gamma$ being real, we conclude that

$$ |\mathcal{E}_0| = |\mathcal{E}_{m+1}|. \qquad (75) $$
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Summary of our physical model: *m* convective steps of constant size *d* and indexed by the integer *n* = {1, ..., *m*} are separated by discrete interfaces. *a*) General scheme: the incident, reflected and transmitted waves have amplitudes *A*₀, *B*₀ and *A*<sub>*m*+1</sub>, respectively. In the *n*th convective step, the ingoing wave has an amplitude *A*<sub>n</sub>, and the outgoing has an amplitude *B*<sub>n</sub>. *b*) The corresponding density profile: at each interface, the density undergoes a discontinuous jump by a value Δρ > 0. In the convective steps, it follows an adiabatic gradient. This creates a mean density gradient |dρ/dz| = Δρ/d. *c*) The corresponding buoyancy frequency profile: *N* = 0 everywhere, except at the location of the interfaces, where it is a Dirac distribution, which creates a mean stratification N̄² = gΔρ/ρ₀d.</figcaption></figure>

The case $\omega \neq f$. The solution for prescribed $k_{\perp}$ and $\omega$ is then

$$ \hat{W}(z) = \begin{cases} \mathcal{A}_0 e^{k_{\perp} z} + \mathcal{B}_0 e^{-k_{\perp} z} & \text{for } z > 0, \\ \mathcal{A}_n e^{k_{\perp} [z+(n-1/2)d]} + \mathcal{B}_n e^{-k_{\perp} [z+(n-1/2)d]} & \text{for } -nd < z < -(n+1)d, n = \{1, \dots, m\} \\ \mathcal{A}_{m+1} e^{k_{\perp} [z+md]} + \mathcal{B}_{m+1} e^{-k_{\perp} [z+md]} & \text{for } z < -md, \end{cases} \quad (76) $$

where

$$ k_z \equiv \sqrt{k_{z,c}^2} = i \sqrt{-k_{z,c}^2}, \quad (77) $$

and the addition of $(n - 1/2)d$ in the exponents have been included to take advantage of symmetry in finding analytic solutions (S16). Thus, defining $k_z^2$ as in Eq. (71) allows us to treat the propagative versus evanescent behaviour of the ingoing and outgoing waves self consistently.

We note that under the traditional approximation, $\tilde{f} = 0$ so that $\tilde{\delta} = 0$. In this sense, the term $\tilde{\delta}$ accounts for non-traditional effects, together with the other non-traditional terms contained in Eq. (71). In the general case, these effects vanish when $f = 0$ or $\tilde{f} = 0$, i.e. at the equator ($\Theta = \pi/2$) and at the pole ($\Theta = 0$), respectively.

Equation (29) being of second order, we need a second boundary condition to get the complete solution. It arises from the requirement that the momentum flux is continuous across the interface. It can be obtained by integrating Eq. (21) across the interface, situated at $z_n = -nd$ for the example case that follows:

$$ \underbrace{\int_{-nd-\epsilon d}^{-nd+\epsilon d} (f^2 - \omega^2) \frac{\partial^2 w}{\partial z^2} dz}_{\alpha_1} + \underbrace{\int_{-nd-\epsilon d}^{-nd+\epsilon d} 2i f \tilde{f}_s k_{\perp} \frac{\partial w}{\partial z} dz}_{\alpha_2} - \underbrace{\int_{-nd-\epsilon d}^{-nd+\epsilon d} k_{\perp}^2 (N^2 - \omega^2 + \tilde{f}_s^2) w dz}_{\alpha_3} = 0, \quad (78) $$

where we have introduced a dimensionless parameter $\epsilon < 1$. We have

$$ \alpha_1 = (f^2 - \omega^2) (w'_n - w'_{n+1}), \quad (79) $$

where the primes denote differentiation with respect to *z*, and

$$ \alpha_2 = 2i f \tilde{f}_s k_{\perp} (w_n - w_{n+1}) = 0, \quad (80) $$

where we have used the first boundary condition $w_n = w_{n+1}$. Finally, we have

$$ \alpha_3 = \int_{-nd-\epsilon d}^{-nd+\epsilon d} k_{\perp}^2 (N^2 - \omega^2 + \tilde{f}_s^2) w dz \\ = k_{\perp}^2 w_n d \left[ \underbrace{\int_{-\epsilon}^{+\epsilon} N^2 d\zeta + (\tilde{f}_s^2 - \omega^2) \times (2\epsilon)}_{\bar{N}^2} \right], \quad (81) $$

where we used the change of variable $\zeta = (z + nd)/d$, and took $w$ out of the integral as it does not vary across the interface. Replacing $\bar{N}^2$ by its expression given by Eq. (68), and taking the limit $\epsilon \to 0$ lead to

$$ \alpha_3 = k_{\perp}^2 g \frac{\Delta \rho}{\rho_0} w_n. \quad (82) $$

Thus, we finally obtain a set of two boundary conditions involving the vertical velocity $w_n$ and its z-derivative $w'_n$,

$$ w_{n+1} = w_n, \quad (83) $$

$$ w'_{n+1} = w'_n + \frac{g \bar{k}^2}{\omega^2} \frac{\Delta \rho}{\rho_0} w_n, \quad (84) $$

where we have defined

$$ \bar{k} = k_{\perp} \frac{\omega}{\sqrt{\omega^2 - f^2}}. \quad (85) $$
These results agree with those of BQF15 for the two cases they consider of the pole and the equator.

Let us consider the interior of the staircase, i.e. consider $n = \{1, \dots, m\}$. Substituting the expression of $w_n$ given by Eq. (76) into the two boundary conditions given by Eqs. (83) and (84), we obtain two recurrence relations between the coefficients $(\mathcal{A}_n, \mathcal{B}_n)$,

$$ \mathcal{A}_n = e^{-i\tilde{\varphi}} [\Delta(1-\Gamma)\mathcal{A}_{n+1} - \Gamma\mathcal{B}_{n+1}], \quad (86) $$

$$ \mathcal{B}_n = e^{-i\tilde{\varphi}} [\Gamma\mathcal{A}_{n+1} + \Delta^{-1}(1+\Gamma)\mathcal{B}_{n+1}], \quad (87) $$

where the following dimensionless quantities have been defined:

$$ \tilde{\varphi} \equiv \delta k_{\perp} d, \quad (88) $$

$$ \Delta \equiv e^{k_{\perp} d}, \quad (89) $$

$$ \Gamma \equiv \frac{1}{2} \frac{g \tilde{k}^2}{k_z \omega^2} \frac{\Delta \rho}{\rho_0} = \frac{1}{2} \left( \frac{\tilde{N}}{\omega} \right)^2 \frac{(\tilde{k}d)^2}{k_z d}. \quad (90) $$

The coefficient $\Gamma$ can further be expressed as a function of more appropriate variables for our problem, as

$$ \Gamma = \frac{1}{2} (k_{\perp} d) \left(\frac{\tilde{N}}{\omega}\right)^2 \left[1 - \left(\frac{2\tilde{\Omega}}{\omega}\right)^2\right]^{-1/2}. \quad (91) $$

Inside the staircase, the coefficients of adjacent convective steps are then related by

$$ \begin{bmatrix} \mathcal{A}_n \\ \mathcal{B}_n \end{bmatrix} = \tilde{\mathbf{T}} \begin{bmatrix} \mathcal{A}_{n+1} \\ \mathcal{B}_{n+1} \end{bmatrix}, \quad (92) $$

where the transfer matrix $\tilde{\mathbf{T}}$ is defined by

$$ \tilde{\mathbf{T}} = e^{-i\tilde{\varphi}} \begin{bmatrix} \Delta(1-\Gamma) & -\Gamma \\ \Gamma & \Delta^{-1}(1+\Gamma) \end{bmatrix}. \quad (93) $$

Above and below the staircase, we have to perform a separate calculation because the stratification is not the same a priori. This will give us the boundary conditions of the entire staircase.

### 3.3. Dispersion relation for gravito-inertial waves in a staircase density profile

The dispersion relation for pure gravity modes in a staircase modelled as described in Sect. 3.1 has been derived by BQF15. This was done in detail in the case without rotation, and they discussed the effects of rotation for the two particular cases where the spin axis is parallel or perpendicular to $\hat{e}_z$ (respectively $\Theta = 0$ at the pole and $\Theta = \pi/2$ at the equator). In this section, we extend their calculation to the general case where the spin axis makes an arbitrary angle $\Theta$ to the local radial direction, in order to obtain the free modes of a density staircase at any latitude.

#### 3.3.1. Infinite staircase

Following BQF15, we first consider an infinite staircase. The local model strictly loses validity in this case, but this idealised calculation will provide us useful insights. We first consider strictly periodic boundary conditions. That is, we assume that there is an integer $m$ so that $\mathcal{A}_n = \mathcal{A}_{n+m}$ and $\mathcal{B}_n = \mathcal{B}_{n+m}$, so that $m$ is the

periodicity of the infinite staircase. Recalling Eq. (92), we obtain

$$ \begin{bmatrix} \mathcal{A}_n \\ \mathcal{B}_n \end{bmatrix} = \tilde{\mathbf{T}}^m \begin{bmatrix} \mathcal{A}_n \\ \mathcal{B}_n \end{bmatrix}. \quad (94) $$

Non trivial solutions of this equation exist only if

$$ \det(\tilde{\mathbf{T}}^m - \delta) = 0, \quad (95) $$

where $\delta$ is the $2 \times 2$ identity matrix. With some algebra, following BQF15, we obtain the dispersion relation for periodic solutions with rotation,

$$ \omega^2 = \bar{N}^2 \left( \frac{(\bar{k}d)^2/k_z d}{2\coth(k_z d) - 2 \cos \theta \censch(k_z d)} \right), \quad (96) $$

where

$$ \theta = \frac{2\pi n}{m} \pm \delta(k_{\perp} d), \quad (97) $$

$n$ being an integer that ranges from 0 to $m-1$, and $\bar{N}^2$ is the background buoyancy frequency that corresponds to the averaged density gradient, defined by Eq. (68). Equation (96) gives the frequencies for the modes of the staircase.

In the case without rotation, $\delta = 0$ and $k_z d = \bar{k}d = k_{\perp}d$, so that we recover BQF15's dispersion relation (see their Eq. (18)):

$$ \omega^2 = \bar{N}^2 \left( \frac{k_{\perp} d}{2\coth(k_{\perp} d) - 2 \cos \theta \censch(k_{\perp} d)} \right). \quad (98) $$

Equation (96) also agrees with BQF15 for the specific cases with rotation that they study, i.e. at the pole and the equator.

#### 3.3.2. Finite staircase embedded in a convective medium

In this section, we extend the calculation made in BQF15 for the case of a finite staircase embedded in a convective medium, including rotation at any colatitude $\Theta$. This is done in order to get the free modes of oscillation of a finite staircase in that setup. The perturbations are assumed to decay as $z \to \pm\infty$. The boundary conditions are then $\mathcal{A}_0 = 0$ at the top, and $\mathcal{B}_{m+1} = 0$ at the bottom.

Applying those boundary conditions and using Eq. (92), we can write

$$ \begin{bmatrix} \mathcal{A}_{m+1} \\ 0 \end{bmatrix} = \tilde{\mathbf{T}}^m \begin{bmatrix} 0 \\ \mathcal{B}_0 \end{bmatrix}. \quad (99) $$

The equation above is true in general if the lower-right corner of the $2 \times 2$ matrix $\tilde{\mathbf{T}}^m$ is zero. The dispersion relation thus becomes

$$ (\tilde{\mathbf{T}}^m)_{22} = 0. \quad (100) $$

We can diagonalise this matrix, and with some algebra (following BQF15), we can obtain a similar dispersion relation:

$$ \omega^2 = \bar{N}^2 \left( \frac{(\bar{k}d)^2/k_z d}{2\coth(k_z d) - 2 \cos \theta \censch(k_z d)} \right), \quad (101) $$

but now $\cos \theta$ is one of the roots of the polynomial

$$ T_m(\cos \theta) + [\cos \theta \coth(k_z d) - \censch(k_z d)] U_{m-1}(\cos \theta) = 0, \quad (102) $$

where $T_m$ and $U_m$ are Chebyshev polynomials of the first and second kinds, respectively, defined by

$$ \cos(m\theta) = T_m(\cos \theta), \quad (103) $$

$$ \sin(m\theta) = U_{m-1}(\cos \theta) \sin \theta. \quad (104) $$

We refer the reader to BQF15 to see the details of the calculation in the case without rotation. Here, the same lines are reproduced. However, taking into account rotation at any colatitude somewhat complicates the analysis.
### 3.4. Transmission of an incident (gravito-)inertial wave

In this section, we aim to answer the following question: what is the effect of a density staircase on the transmission of an incident internal wave with amplitude $\mathcal{A}_0$ taking into account the complete Coriolis acceleration? This should help us to understand how the presence of layered semi-convection and associated density staircases would affect the fate of tidally excited waves launched in one region as it propagates towards another. It could also be of great interest when studying the seismology of giant planets, to predict what modes can be observed when looking at the oscillations of their surfaces (Gaulme et al. 2011; Fuller 2014).

We thus analyse the properties of the transmission of an internal wave upon a finite-length density staircase. We will use the same formalism as previously, except for the boundary conditions. Indeed, because we consider the transmission of an internal wave upon the staircase, we refer to Fig. 5: an incident (ingoing) wave with amplitude $\mathcal{A}_0$ enters the staircase, creating a reflected (outgoing) wave with amplitude $\mathcal{B}_0$. Only a downward propagating wave is assumed to exist below the staircase. This is the transmitted wave, with amplitude $\mathcal{A}_{m+1}$ ($\mathcal{B}_{m+1} = 0$).

The first work of this type was S16, who studied internal wave transmission through a density staircase in the ocean using the traditional approximation. If this is often appropriate for a thin oceanic layer on the Earth, this is not suitable to model low-frequency internal waves in the deep envelopes of giant planets (Ogilvie & Lin 2004). In this section, we generalise S16 to arbitrary top and bottom layer properties, and include the complete Coriolis acceleration.

#### 3.4.1. Analytic expression of the transmission coefficient

Our aim is now to calculate a transmission coefficient in the general case of arbitrary boundary conditions and number of steps. This will allow us to determine what fraction of the incident wave energy makes its way through the staircase and propagates to deeper regions of a giant planet. Recalling the expression we derived in Sect. 2.4.3, the transmission coefficient is defined by

$$T = \frac{k_{z,b}}{k_{z,a}} \left| \frac{\mathcal{A}_{m+1}}{\mathcal{A}_0} \right|^2, \qquad (105)$$

which simply reduces to $T = |\mathcal{A}_{m+1}/\mathcal{A}_0|^2$ when the stratification is the same above and below the staircase, for example in the case of a staircase embedded in a stably stratified medium (S16).

Substituting Eq. (76) (including the boundary condition discussed above, namely $\mathcal{B}_{m+1}$) into the interface conditions given by Eqs. (83) and (84), we obtain after some algebra the following set of equations:

$$\begin{aligned} \mathcal{A}_0 = & \frac{1}{2} \Delta^{1/2} [1 - K_a^c (1 - \Gamma)] \mathcal{A}_1 \\ & + \frac{1}{2} \Delta^{-1/2} [1 + K_a^c (1 + \Gamma)] \mathcal{B}_1 \end{aligned} \qquad (106)$$

$$\mathcal{A}_n = \Delta(1 - \Gamma) \mathcal{A}_{n+1} - \Gamma \mathcal{B}_{n+1} \qquad (107)$$

$$\mathcal{B}_n = \Gamma \mathcal{A}_{n+1} + \Delta^{-1}(1 + \Gamma) \mathcal{B}_{n+1} \qquad (108)$$

$$\mathcal{A}_m = \frac{1}{2} \Delta^{1/2} [1 - \Gamma + K_c^b] \mathcal{A}_{m+1} \qquad (109)$$

$$\mathcal{B}_m = \frac{1}{2} \Delta^{-1/2} [1 + \Gamma - K_c^b] \mathcal{A}_{m+1}, \qquad (110)$$

where $n = \{1, \dots, m-1\}$, and we have defined

$$K_{\beta}^{\alpha} = \frac{k_{z,\alpha}}{k_{z,\beta}}. \qquad (111)$$

We recall that it has been defined that $k_z \equiv k_{z,c}$ in convective steps, to make notations lighter. Those equations can be combined to express $\mathcal{A}_0$ in term of solely $\mathcal{A}_{m+1}$,

$$\mathcal{A}_0 = [b_a^T \tilde{T}^{m-1} b_b] \mathcal{A}_{m+1}, \qquad (112)$$

where $^T$ denotes transposition, and the left and right vectors are defined by

$$b_a = \frac{1}{2} \left[ \frac{\Delta^{1/2} [1 - K_a^c (1 - \Gamma)]}{\Delta^{-1/2} [1 + K_a^c (1 + \Gamma)]} \right], \qquad (113)$$

$$b_b = \frac{1}{2} \left[ \frac{\Delta^{1/2} [1 - \Gamma + K_c^b]}{\Delta^{-1/2} [1 + \Gamma - K_c^b]} \right]. \qquad (114)$$

Therefore, the transmission coefficient is given by

$$T = \frac{k_{z,b}}{k_{z,a}} |b_a^T \tilde{T}^{m-1} b_b|^2, \qquad (115)$$

with

$$\frac{k_{z,b}}{k_{z,a}} = \left( \frac{\omega^4 + \omega^2[N_b^2 - 4\tilde{\Omega}^2] + (N_b f)^2}{\omega^4 + \omega^2[N_a^2 - 4\tilde{\Omega}^2] + (N_a f)^2} \right)^{1/2}, \qquad (116)$$

where we recall that $2\tilde{\Omega} \equiv \sqrt{f^2 + f_c^2}$ (i.e. $2\Omega$ when $\alpha = \pi/2$).

Equation (115) is a general analytic expression of the transmission coefficient, that can now be used in order to analyse the behaviour of the transmission in the parameter space (wave frequency, wavelength, number of steps, boundary conditions, etc.).

#### 3.4.2. Transmission across one step

Staircase embedded in a stably stratified medium. S16 have considered the transmission of an internal wave through a density staircase embedded in a stably stratified medium ($N_a, N_b > 0$). Because the aim of his work is to predict energy transport by internal waves incident upon observed density staircases in the strongly stratified ocean, the traditional approximation is adopted (even though Gerkema & Shrira 2005, have demonstrated that this is not always appropriate, even for the ocean). We begin by recovering their results in the case of the transmission across one step (see his Sect. II.B) under the traditional approximation, which corresponds to a colatitude $\Theta = 0$, where the rotation axis and gravity are aligned. This will provide us a first check that our mathematical formalism, in which the traditional approximation is not assumed, is correct. The result is shown in Fig. 8.

The transmission coefficient $T$ is displayed as a function of the wave frequency normalised by the mean buoyancy frequency, $\omega/\tilde{N}$, and of the horizontal wave number rendered dimensionless by multiplication by the size of the convective layers, $k_{\perp}d$. In this and subsequent figures, the range of frequency over which the transmission coefficient is calculated is chosen so that both the incident and reflected waves (carrying energy densities $\propto|\mathcal{A}_0|^2$ and $\propto|\mathcal{A}_{m+1}|^2$, respectively) are propagative in both regions $z > 0$ and $z < -D$, respectively. This gives a range $\omega_- < \omega < \omega_+$, where

$$\omega_- = \max(\omega_-(N_a), \omega_-(N_b)), \qquad (117)$$

$$\omega_+ = \min(\omega_+(N_a), \omega_+(N_b)), \qquad (118)$$

where $\omega_{\pm}(N_{\alpha})$ is given by Eq. (33). Figure 8 can be directly compared to Fig. 2a of S16 for parameters $\Omega = 0.05\tilde{N}$, $\Theta = 0$,
<figure><img src="image_11.png" /><figcaption>Fig. 8. Transmission coefficient as a function of normalised frequency, $\omega/\bar{N}$ and horizontal wave number, $k_{\perp}d$, with $\Omega = 0.05\bar{N}$ and $\Theta = 0$ ($f = 0.1\bar{N}$), for comparison with S16. The transmission coefficient, $T$, lies between 0 and 1 and is calculated for the range of frequency over which the incident wave is propagative, which we denote by $\omega_{-} < \omega < \omega_{+}$.</figcaption></figure>

giving $f = 0.1\bar{N}$. No transmission occurs for $\omega < f$ because the incident wave would be evanescent in that case, as explained above (this occurs in the traditional approximation).

We can see that the qualitative behaviour of the transmission coefficient is the same: transmission is enhanced for large wave-lengths (low wave numbers) compared to the size of the staircase ($D = d$ in this section because there is only one step), and there is a branch of enhanced transmission near $\omega_{+}$, where transmission can be large even with $k_{\perp}d$ of order unity.

However, a quantitative discrepancy is found between the two figures: perfect transmission for incident waves with $\omega \approx \omega_{+}$ occur for $k_{\perp}d \approx 2.4$ in the case of S16, and $k_{\perp}d \approx 1.5$ in our case. This difference is due to the fact that the density jumps at the first and the last interfaces is taken to be $\Delta\rho/2$ in S16, and $\Delta\rho$ everywhere in our work in order to treat all cases self-consistently, regardless of the boundary conditions applied to the staircase. However, it has been checked that the model of S16 with our formalism yields his results. The origin of this band of enhanced transmission is a resonance between the ingoing wave with interfacial gravity waves on either side of the convective step. The dispersion relation for these waves on an interface with density jump $\Delta\rho$ (assuming infinitely deep layer) is

$$ \omega^2 = \frac{1}{2} \frac{k_{\perp} g \Delta \rho}{\rho_0} = \frac{1}{2} (k_{\perp} d) \bar{N}^2. \qquad (119) $$

So when $\omega \approx \bar{N}$, we expect enhanced transmission when $k_{\perp}d \approx 2$, which is what we (and S16) observe, at least approximately. However, given that our layer is not infinitely deep, we expect quantitative discrepancies from this simple estimate.

We note that waves with low wave numbers, and thus large wavelengths

$$ \lambda \gg d \qquad (120) $$

are expected to be unaffected by the staircase. This is confirmed by the fact that we get perfect transmission at low wave numbers, for any frequency. Those waves will see the staircase as a continuously stratified medium with $N_0 = \bar{N}$. However, shorter

wavelength waves, with $\lambda \sim d$, are strongly affected by the presence of a staircase, only being efficiently transmitted in certain regions of parameter space.

Now, to quantify the importance of the non-traditional effects, we move the box away from the pole and choose a colatitude $\Theta = \pi/4$. We calculate the transmission coefficient embedded in a stably stratified medium with $N_a = N_b = \bar{N}$ for various rotation rates to quantify the influence of rotation. The results are displayed in Fig. 9 for $\Omega = 0, 0.4\bar{N}$ and $0.6\bar{N}$. Thus, in an astrophysical context we mostly consider fast rotators, so that the effects of the complete Coriolis acceleration are clearer, especially near $\omega = f$. In Fig. 9, we see that new features arise compared to the traditional case (see Fig. 8), in the form of bands of perfect transmission departing from $\omega = f$ ($\approx 0.57\bar{N}$ for $\Omega = 0.4\bar{N}$, and $\approx 0.85\bar{N}$ for $\Omega = 0.6\bar{N}$), these getting thicker with increasing rotation (but not more numerous, as we will see later). We note that the discontinuous (scattered, even) appearance of some of these branches is a plotting issue, because the branches are very narrow. However, the underlying branches are physical and not due to numerical errors, as we will discuss below.

Figure 10 shows the behaviour of the transmission coefficient as a function of the vertical wave number $k_z$ (instead of $k_{\perp}$), rendered dimensionless by multiplication by $d$. It will be seen that as soon as rotation is non-zero (see Figs. 10b, c), we obtain perfect transmission around the Coriolis frequency $\omega = f$, which is indicated by a white arrow. Furthermore, two regions can be identified:

- the region $\omega_{-} < \omega < 2\bar{\Omega}$ (delimited by the vertical dashed red line) – in which transmission is enhanced even for large vertical wave-lengths (meaning low vertical wave numbers) – shows bands of perfect transmission departing from the vertical transmission line at $\omega = f$, to reach $k_zd$ equalling multiple of $\pi$ lines (indicated by the horizontal dashed black lines) near $2\bar{\Omega}$;

- the region $2\bar{\Omega} < \omega < \omega_{+}$, which shows a similar behaviour to Fig. 9 as a function of the horizontal wave number $k_{\perp}$.

We note that the region $\omega_{-} < \omega < 2\bar{\Omega}$ corresponds to the frequency window for which both inertial waves in the convective layers, and GIWs in the stably stratified ones, are propagative. Thus, it makes sense to obtain enhanced transmission for those frequencies. We note that $\bar{\Omega} = \Omega$ for $\alpha = \pi/2$ that is adopted throughout this and subsequent sections.

Thus, the transmission differs mostly for GIWs, which have frequency close to $f$. We find that the qualitative behaviour remains the same for $\omega > 1.5f$. However, because including the complete Coriolis acceleration extends the range of admissible frequencies below $f$, there is a new band of perfect transmission for $k_{\perp}d$ of order unity arising for sub-inertial waves (i.e. for $\omega_{-} < \omega < f$). Similarly, the range of admissible frequencies is extended above $N_{a,b} = \bar{N}$, so that the band of perfect transmission for $\omega \approx \omega_{+}$ reaches higher values of $k_{\perp}d$ with increasing rotation rates.

The features described above only arise when the complete Coriolis force is taken into account, and therefore they are intrinsically caused by non-traditional effects. We will give them a physical interpretation in Sect. 3.4.4.

Staircase embedded in a convective medium. Now, we consider a staircase embedded in a convective medium, which is probably more relevant to a portion of a giant planet interior in which double-diffusive convection operates. In this case, the incident and transmitted waves are pure inertial waves.
<figure><img src="image_12.png" /><figcaption>Fig. 9. Transmission coefficient for one step embedded in a stably stratified medium, in an inclined box with $\Theta = \pi/4$, for different rotation frequencies a) $\Omega = 0$; b) $\Omega = 0.4\bar{N}$ and c) $\Omega = 0.6\bar{N}$, as a function of frequency and perpendicular wave number. The white and dashed red arrows indicate frequencies $\omega = f$ and $\omega = 2\Omega$, respectively. A set of bands depart from $f = 2\Omega \cos \Theta (=\sqrt{2}\Omega$ here) corresponding with a resonance with short-wavelength inertial waves. The transmission coefficients are calculated over a frequency range $\omega_{-} < \omega < \omega_{+}$, calculated such that both the incident and transmitted wave are propagative.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 10. Transmission coefficient for one step embedded in a stably stratified medium, in an inclined box with $\Theta = \pi/4$, for different rotation frequencies a) $\Omega = 0$, b) $\Omega = 0.4\bar{N}$ and c) $\Omega = 0.6\bar{N}$, as a function of frequency and of the vertical wave number. The white arrow indicates the frequency $\omega = f$, while the vertical dashed red line indicates the frequency $\omega = 2\Omega$. The horizontal dashed black lines indicates multiple of $\pi$.</figcaption></figure>

Considering a staircase embedded in a convective medium can also be seen as considering a portion of a vertically (more) extended staircase.

In the special case of a single convective layer embedded in convective medium ($m = 1, N_a = N_b = 0$), the expression of the transmission coefficient, $T$ given by Eq. (115) takes the following analytical expression:

$$T = \frac{1}{1 + 4\Gamma^2 [\cos(k_z d) - \Gamma \sin(k_z d)]^2}, \quad (121)$$

where $k_z$ is the vertical wave number of inertial waves,

$$k_z = k_{\perp} \sqrt{\frac{\omega^2 (4\bar{\Omega}^2 - \omega^2)}{(\omega^2 - f^2)^2}}. \quad (122)$$

From this formula, it is obvious that $T$ lies between 0 and 1, as we expect for a transmission coefficient.

In this case, the range of frequency over which the transmission coefficient is calculated is the range of frequencies over which inertial waves are propagative in a convective medium, i.e.

between 0 and $2\bar{\Omega}$ (see Eq. (34)). In the case where $\alpha = \pi/2$ that we consider here, $\bar{\Omega} = \Omega$. The results are displayed in Fig. 11 for different rotation rates $\Omega = 0.2\bar{N}$, $\Omega = 0.4\bar{N}$, $\Omega = 0.6\bar{N}$, and the colatitude $\Theta = \pi/4$. Again, we get perfect transmission for large enough wavelengths (small wave numbers), and a set of bands depart from $\omega = f (\approx 0.28\bar{N}$ for $\Omega = 0.2\bar{N}$, $\approx 0.57\bar{N}$ for $\Omega = 0.4\bar{N}$, and $\approx 0.85\bar{N}$ for $\Omega = 0.6\bar{N}$). What differs from the previous case is that transmission stays close to unity for higher and higher wave numbers as rotation is increased, which made us choose a range in $k_{\perp}d$ from 1.0 to 2.5 for $\Omega = 0.2\bar{N}$ and $\Omega = 0.6\bar{N}$, in order to observe this effect.

Figure 12 shows the behaviour of the transmission coefficient as a function of the vertical wave number, $k_z$. Again, we find perfect transmission near $\omega = f$ (white arrow), and a set of bands of perfect transmission are departing from this frequency to reach lines of $k_zd = 0, \pi, 2\pi, \dots$ (black dashed lines) for any vertical wavelength.

Staircase embedded in convective medium at the top and stably stratified medium at the bottom. The boundary conditions
<figure><img src="image_14.png" /><figcaption>Fig. 11. Same as Fig. 9, but for one step embedded in a convective medium. The frequency domain extends from 0 to 2Ω.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Same as Fig. 10, but for one step embedded in a convective medium. The frequency domain extends from 0 to 2Ω.</figcaption></figure>

of the staircase are changed once again to consider a staircase embedded in a convective medium at the top and a stably stratified medium at the bottom. Thus, the incident wave is a pure inertial wave, and the transmitted wave is a gravito-inertial wave. This situation is interesting to consider for astrophysical applications because density staircases might develop in the deep interiors of giant planets where there could be a stably stratified region just outside the core (Fuller 2014), on top of which could sit a region of layered semi-convection and associated density staircases.

This is found to be equivalent to the opposite situation where $N_a$ and $N_b$ are interchanged, which could correspond to a stably stratified layer near the surface of a hot Jupiter matching the convective envelope of deeper regions through a region of layered semi-convection. The fact that both cases are identical likely results from the symmetry of the Boussinesq equations under the transformation {$z \to -z, \Theta \to -\Theta, b \to -b, w \to -w$} (i.e. “up \to down”, “hot \to cold”).

The results are displayed in Figs. 13 and 14, as a function of $k_{\perp}d$ and $k_zd$, respectively, for different rotation rates $\Omega = 0.2\bar{N}$, $\Omega = 0.4\bar{N}$, $\Omega = 0.6\bar{N}$. The behaviour of the transmission coefficient is similar to the two previous cases. One difference that arises, however, concerns the band of perfect transmission corresponding to resonance with interfacial gravity waves on either side of the convective step, which do not extend to $\omega_+$ like in the two previous cases.

### 3.4.3. Transmission across *m* steps

Density staircases that might exist in giant planet interiors are expected to have a large number of steps, with each step being much smaller than the planet's radius (Leconte & Chabrier 2012). Therefore, we now study the effects of having more than one step. We choose again to perform our calculations at a colatitude $\Theta = \pi/4$, for illustration and comparison with the results of the previous section. The transmission coefficient is given by Eq. (115), and we checked that our results agreed with S16 assuming the traditional approximation, and for similar parameters. We recall that here, unlike in S16, the complete Coriolis acceleration is taken into account.

Figure 15 shows the results for $m = 2, 5$ and 10 steps embedded in a stably stratified medium with $\Omega = 0.4\bar{N}$ (upper panels) and in a convective medium with $\Omega = 0.6\bar{N}$ (bottom panels). The upper and bottom panels can directly be compared to Figs. 9b and 11c, respectively, on which is displayed the corresponding single step case with same colatitude and rotation rates.

In the regions $\omega > f$ and $\omega < f$, it will be seen that as the number of steps increases, more and more transmission peaks appear. These are particularly visible for $\omega \le \omega_+$ ($\omega \ge \omega_-$, respectively), and propagate back (forward, respectively) to $\omega = f$. Confirming the results of S16 in the case of a staircase embedded in a stably stratified medium, we find that if there are *m* steps, the number of peaks equals *m*. This will be explained in the following section. Also, as the numbers of steps increases,
<figure><img src="image_16.png" /><figcaption>Fig. 13. Same as Fig. 9, but for one step embedded in a convective medium at the top and a stably stratified medium at the bottom.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Same as Fig. 10, but for one step embedded in a convective medium at the top and a stably stratified medium at the bottom.</figcaption></figure>

those bands of perfect transmission become narrower, and transmission becomes inhibited for larger and larger wavelengths ($k \to 0 \Rightarrow \lambda \to \infty$).

As a function of the vertical wave number (see the bottom panel of Fig. 15), those observations hold. Namely, in each region of enhanced transmission previously described, one can see m bands of perfect transmission for a staircase with m steps.

### 3.4.4. Understanding the transmission coefficient

Perfect transmission at $\omega = f$. We have seen that the transmission is perfect at the Coriolis frequency $\omega = f$ for almost any wavelength (see for example Fig. 12). This result can be interpreted by noting that when $\omega = f$ there is only one direction for wave propagation that is allowed by the dispersion relation that is not parallel to the interface (see Fig. 3: the group velocity $v_g$ is directed along the perpendicular to the z-axis), so there is only one wave solution with non-zero vertical velocity. This means that the boundary condition at the interface that matches the vertical velocity between the incoming and transmitted waves (see Eq. (131)) requires the transmitted wave to have the same amplitude, so that $T = 1$. Furthermore, from the separate calculation we have made in Sect. 3.2 for $\omega = f$ (see Eq. (75)), it is straightforward that $T = |\varepsilon_{m+1}/\varepsilon_0|^2 = 1$ (at least when the stratification is the same above and below the staircase).

Resonances with short-wavelength inertial modes. In Sects. 3.4.2 and 3.4.3, we have seen that, regardless of the properties of the regions that surround the staircase, we always obtained a set of bands of perfect transmission departing from the inertial frequency $\omega = f$. We interpret those peaks of transmission as being caused by resonances between the incident internal wave and the short-wavelength inertial waves that exist within a convective step. This happens when a multiple of the vertical semi wavelength ($\lambda_z/2$) of the internal waves that propagate inside the step fits inside one step, i.e. for

$$ \frac{n\lambda_z}{2} = d, \qquad (123) $$

where n is an integer. This is illustrated in Fig. 17 for $n = \{1, 2, 3\}$. The condition above can be rewritten in term of the vertical wave number of the inertial wave inside the steps, to read

$$ k_z d = n\pi. \qquad (124) $$

In Fig. 3.4.5, curves of equation $k_z d = n\pi$ (in dashed red) are overplotted to Fig. 9c for $n = \{1, \dots, 6\}$. As one can see on this figure, the matching between the bands of transmission and the curves of equation $k_z d = n\pi$ is satisfying. The largest discrepancy is obtained for the peak of transmission in the bottom left corner of the figure, which corresponds to very low frequency.
<figure><img src="image_18.png" /><figcaption>Fig. 15. Transmission coefficient for a staircase embedded in a stably stratified medium, in an inclined box with $\Omega = 0.4\bar{N}$ and $\Theta = \pi/4$, for different number of steps a) $m = 2$; b) $m = 5$ and c) $m = 10$, as a function of dimensionless frequency and horizontal wave number (upper panel) or vertical wave number (bottom panel). The white arrow indicates the frequency $\omega = f$. A set of bands extending vertically depart from $f = 2\Omega \cos \Theta$ ($\approx 0.57\bar{N}$). In addition, m bands of perfect transmission appear in the regions $\omega < f$ and $\omega > f$, each becoming narrower with increasing m.</figcaption></figure>

Excitation of free gravity modes of the staircase. Other features we give a physical interpretation for are the bands of perfect transmission, more noticeable for the $m > 1$ steps case with m bands of perfect transmission starting from $\omega = \omega_+$ and $k_{\perp}d$ of order unity, and propagating back to $\omega = f$ and $k_{\perp}d \ll 1$ (see Fig. 15). In Sect. 3.3, we have derived dispersion relations, whose branches give us the free modes of oscillation of the staircase. Therefore, incident waves with a prescribed frequency can resonate with a free mode of oscillation of the staircase if its frequency is such that it matches a root of the dispersion relation given by Eq. (101) (in the case of a finite staircase).

Since it is very challenging to extract the roots of the dispersion relation with rotation (in particular because there are infinitely many solutions near $\omega = f$, which we have just identified), we have focused on the case without rotation, for which the dispersion relation given by Eq. (101) gives us a direct relation between $\omega/\bar{N}$ and $k_{\perp}d$. In Fig. 18b, the transmission coefficient in the case of 10 steps without rotation is displayed, together with the branches of the corresponding dispersion relation. As one can see, the bands of perfect transmission are matched satisfyingly well by the branches of the dispersion relation (in light blue). They correspond to incident waves that excite a free mode of oscillation of the staircase. In principle, the same could be done for the case with rotation, the only reason this has not been done here being that it is hard to extract all of the roots of

the dispersion relation, but we have nevertheless identified the branches near $\omega = f$ as explained above.

### 3.4.5. Non-uniform step sizes

In this section, we relax the assumption of having equally sized steps. Observed density staircases in the ocean have approximately equally-sized steps with a typical length-scale, $d \sim 2.5$ m (see Ghaemsaidi et al. 2016), which determines the vertical extent of the convective steps. We might expect density staircases in giant planet interiors also to have approximately equally-sized steps, though this is not clear from theory or observations. However, we can expect (as is observed in the case of the ocean, see Ghaemsaidi et al. 2016) some inherent variations in the vertical extent of each step around a mean value d. To model them, we adopt the same procedure as in S16: inside the staircase, the steps have a vertical extent

$$d_n = d(1 + \epsilon\sigma_n), \quad (125)$$

where, for $1 \le n < m$, $\sigma_n$ are random numbers between -1 and 1, $\epsilon$ being the amplitude of the fluctuations of the vertical extent of the convective layers about the mean value d. Typically, we set $0 \le \epsilon < 0.5$, the value $\epsilon = 0$ corresponding to equally spaced interfaces, i.e. our reference model. It is interesting to see how
<figure><img src="image_19.png" /><figcaption>Fig. 16. Same as Fig. 15, but for a staircase embedded in a convective medium, with a rotation rate $\Omega = 0.6\bar{N}$ ($f \approx 0.85\bar{N}$).</figcaption></figure>

the value of the parameter $\epsilon$ affects the results that have been described so far, because having uneven steps is a more realistic situation. The transmission coefficient can be calculated numerically, and is given by Eq. (115), where $\bar{T}^{m-1}$ is now given by

$$ \prod_{n=1}^{m-1} \left[ \begin{array}{cc} (1 + \Gamma) \Delta^{-1} & \Gamma \Delta^{-2\epsilon\sigma_n} \\ -\Gamma \Delta^{2\epsilon\sigma_n} & (1 - \Gamma) \Delta \end{array} \right] e^{-i\bar{\phi}_n}, \quad (126) $$

where $\bar{\phi}_n = \delta k_{\perp} d_n$. The matrix given by Eq. (126) reduces to $\bar{T}^{m-1}$ when $\epsilon = 0$.

The results obtained in the case of five steps ($m = 5$) for $\Omega = 0.4\bar{N}$, $\Theta = \pi/4$ and different values of $\epsilon = \{0, 0.1, 0.5\}$ are displayed in Fig. 19. The left panel ($\epsilon = 0$) is used as a test of the ability of the method described above to reproduce the analytical solution given by Eq. (115), and is useful for comparison with cases which have $\epsilon > 0$. The cases with $\epsilon = 0.1$ and $\epsilon = 0.5$ show no difference at large wavelengths, so that the perfect transmission for wavelengths $\lambda \gg d$ is a robust result.

However, some differences arise for larger wave numbers (i.e. smaller wavelengths). In particular, the bands of perfect transmission departing from $\omega = f$ weaken with increasing $\epsilon$ (especially the rightmost one). These bands corresponding to half multiple of the wavelength fitting perfectly inside a step (as explained in Sect. 3.4.4), we expect this mechanism to be less effective when the steps are no longer even, since these bands are due to resonances with waves that strongly depend on the step-size. Thus, this observation is perfectly in agreement with the physical interpretation of the bands near the inertial frequency.

<figure><img src="image_20.png" /><figcaption>Fig. 17. Transmission is enhanced for incident waves for which a multiple of the vertical semi-wavelength $\lambda_c/2 = \pi/k_c$ matches the vertical semi-wavelength of inertial waves that fits inside one step, i.e. $k_c d = n\pi$ with $n$ an integer. This is shown for $n = \{1, \dots, 3\}$.</figcaption></figure>

The rest of the figure does not differ much from the even step case until the value $\epsilon = 0.5$ is reached (right panel), for which it becomes obvious that there is less transmission near $\omega = \omega_+$. Some modulation of the placement and intensity of the perfect transmission bands between $f$ and $\omega_+$ can also be seen. This means that free modes of the staircase, given by the dispersion relation, are weakened when the steps are no longer even. This is expected to vary for each random realisation, which we have indeed verified.
<figure><img src="image_21.png" /><figcaption>Fig. 18. a) Same as Fig. 9c, overplotted with curves of equation $k_z d = n\pi$ for $n = \{1, \dots, 6\}$. A set of bands of perfect transmission departs from the frequency $f = 2\Omega \cos \Theta$, corresponding with a resonance with short-wavelength inertial waves. b) The case of 10 steps without rotation, overplotted with the branches of the dispersion relation given by Eq. (101). Transmission is enhanced along those branches, corresponding to resonances with free modes of the staircase.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. a) The even steps case and two uneven steps cases for one particular realisation with b) $\epsilon = 0.1$ and c) $\epsilon = 0.5$. These results have been obtained for 5 steps embedded in a stably stratified medium, with $\Omega = 0.4\bar{N}$ and $\Theta = \pi/4$. The white and red dashed arrows indicate $\omega = f$ and $\omega = 2\Omega$, respectively.</figcaption></figure>

## 3.5. Extension of the model: finite size interfaces

### 3.5.1. Physical set up

Now, we relax the assumption of having discrete interfaces. We denote by $l$ the vertical extent of the stably stratified interfaces, and define the aspect ratio

$$ \epsilon = \frac{l}{d}, \qquad (127) $$

where $d$ is still the vertical extent of the convective layers. The staircase thus consists of a succession of patterns containing a stably stratified and a convective layer (in red and orange respectively, see Fig. 20), of total size $L = l+d$. The total vertical extent of the staircase is $D = mL + l$. This model is of course more realistic since we naturally expect the stably stratified interfaces to have non-zero vertical extent.

The buoyancy frequency profile is defined step wise (see Fig. 20c) through

$$ N^2 = \begin{cases} N_i^2 & \text{for } -nL < z < -(nL + l), \\ n = \{0, \dots, m\}, & \\ 0 & \text{otherwise,} \end{cases} \qquad (128) $$

where $N_i^2$ has been defined in order to keep the mean stratification (and accordingly the mean density gradient, see Fig. 20) the same as in the previous model. This gives

$$ N_i^2 = \frac{\bar{N}^2}{\epsilon}. \qquad (129) $$

This definition of $N_i^2$ ensures that the integrated value of $N^2$ over one interface is independent of $\epsilon$ and equals $\bar{N}^2$. Also, for vanishingly small aspect ratios, we recover the reference model with discrete interfaces. This is illustrated in Fig. 21.
<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Summary of our physical model with finite size interfaces: *m* convective steps of constant size *d* and indexed by the integer *l* < *n* < *m*, are separated by interfaces of size *l*. *a*) General scheme: the incident, reflected and transmitted waves have amplitudes *A*₀, *B*₀ and *A*<sub>*m*+1</sub>, respectively. In the *n*th convective step, the ingoing wave has an amplitude *A*<sub>*n*</sub>, and the outgoing has an amplitude *B*<sub>*n*</sub>, while in the *n*th stably stratified interface they have amplitudes *C*<sub>*n*</sub> and *D*<sub>*n*</sub>, respectively. *b*) The corresponding density profile: in each interface, the density follows a gradient |dρ/dz| = Δρ/l, so that the density varies by an amount Δρ > 0. In the convective steps, it follows an adiabatic gradient. Thus, the mean density gradient |dρ/dz| = Δρ/L, where L = l + d. *c*) The corresponding buoyancy frequency profiles are N² = 0 in the convective layers and N<sub>i</sub>² = N̄²/ε in the interfaces to create the mean stratification N̄² = gΔρ/ρ₀L.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 21.** Buoyancy frequency profiles with various aspect ratios. As the aspect ratio decreases, $N_i$ increases to maintain the integrated value of N equal to $\bar{N}$. Our model converges towards the reference model (the limit being the Dirac distribution) when $\epsilon \to 0$.</figcaption></figure>

The corresponding density profile is shown in Fig. 20b. In convective layers (represented in orange in Fig. 20a), the density follows an adiabatic gradient, while in the interfaces the density gradient is

$$ \left| \frac{d\bar{\rho}}{dz} \right| = \frac{\Delta\rho}{l} = \frac{\Delta\rho}{\epsilon d}, \qquad (130) $$

where the density jump across one interface, Δρ, is maintained to be the same as in the reference model of Sect. 3.4, thanks to this procedure.

### 3.5.2. Analytic expression of the transmission coefficient

Since the density profile is continuous (there is no density jump), the boundary conditions between adjacent layers are simply

the continuity of the vertical velocity and its derivative (see Eq. (78)):

$$ w_n = w_{n+1}, \qquad (131) $$

$$ w'_{n} = w'_{n+1}. \qquad (132) $$

Now, the function $\hat{W}$ reads

$$ \hat{W}(z) = \begin{cases} \mathcal{A}_0 e^{k_{z,a}z} + \mathcal{B}_0 e^{-k_{z,a}z} & \text{for } z > 0, \\ C_0 e^{k_{z,i}z} + \mathcal{D}_0 e^{-k_{z,i}z} & \text{for } -l < z < 0, \\ \mathcal{A}_n e^{k_z(z+nL)} + \mathcal{B}_n e^{-k_z(z+nL)} & \text{for } -(n-1)L - l < z < -nL, n = \{1, \dots, m\}, \\ C_n e^{k_{z,i}(z+nL)} + \mathcal{D}_n e^{-k_{z,i}(z+nL)} & \text{for } z < -(mL + l). \end{cases} \qquad (133) $$

Using the boundary conditions given by (131)–(132), we get recurrence relations between coefficients at each interface between a convective layer and a stably stratified interface that, after little algebra, can be put into matrix form as in Sect. 3.2. We first get a relation between coefficients ($\mathcal{A}_n$, $\mathcal{B}_n$) and ($C_n$, $\mathcal{D}_n$),

$$ \begin{bmatrix} C_n \\ \mathcal{D}_n \end{bmatrix} = \mathbf{T}_1 \begin{bmatrix} \mathcal{A}_n \\ \mathcal{B}_n \end{bmatrix}, \qquad (134) $$

where

$$ \mathbf{T}_1 = \frac{1}{2} \begin{bmatrix} (1 + K_i^c) \Delta_{c,l}^{-1} \Delta_{i,l} & (1 - K_i^c) \Delta_{c,l} \Delta_{i,l} \\ (1 - K_i^c) \Delta_{c,l}^{-1} \Delta_{i,l}^{-1} & (1 + K_i^c) \Delta_{c,l} \Delta_{i,l}^{-1} \end{bmatrix}. \qquad (135) $$

Here, we have defined

$$ \Delta_{\alpha, \beta} = e^{k_{\alpha, \beta} \beta}. \qquad (136) $$
<figure><img src="image_25.png" /><figcaption>Fig. 22. Transmission coefficient for one step embedded in a convective medium with $\Omega = 0.4\bar{N}$ and $\Theta = \pi/4$ for different aspect ratios a) $\varepsilon = 0.01$; b) $\varepsilon = 0.1$ and c) $\varepsilon = 0.5$. The white and dashed red arrows indicate frequency $\omega = f$ and $\omega = 2\Omega$, respectively. A vertical band of perfect transmission appears near the inertial frequency $\omega = f$, this broadens with increasing $\varepsilon$.</figcaption></figure>

Here $\beta$ stands for either $l$ or $L$. Then, we get a second relation between coefficients ($\mathcal{A}_n, \mathcal{B}_n$) and ($C_{n+1}, \mathcal{D}_{n+1}$),

$$ \begin{bmatrix} \mathcal{A}_n \\ \mathcal{B}_n \end{bmatrix} = \bar{\mathbf{T}}_2 \begin{bmatrix} C_{n+1} \\ \mathcal{D}_{n+1} \end{bmatrix}, \qquad (137) $$

where

$$ \bar{\mathbf{T}}_2 = \frac{1}{2} e^{-i\bar{\varphi}} \begin{bmatrix} (1 + K_c^i)\Delta_{c,L} (1 - K_c^i)\Delta_{c,L} \\ (1 - K_c^i)\Delta_{c,L}^{-1} (1 + K_c^i)\Delta_{c,L}^{-1} \end{bmatrix}. \qquad (138) $$

Equations (134) and (137) can be combined in order to get a recurrence relation between coefficients ($C_n, \mathcal{D}_n$):

$$ \begin{bmatrix} C_n \\ \mathcal{D}_n \end{bmatrix} = \mathbf{T}_1 \bar{\mathbf{T}}_2 \begin{bmatrix} C_{n+1} \\ \mathcal{D}_{n+1} \end{bmatrix}. \qquad (139) $$

Finally, as we did in Sect. 3.5.2, we can then express $\mathcal{A}_0$ in term of solely $\mathcal{A}_{m+1}$,

$$ \mathcal{A}_0 = [b_a^T (\mathbf{T}_1 \bar{\mathbf{T}}_2)^m b_b] \mathcal{A}_{m+1}, \qquad (140) $$

where the left and right vectors are defined by

$$ b_a = \frac{1}{2} \begin{bmatrix} (1 + K_a^i) \\ (1 - K_a^i) \end{bmatrix}, \qquad (141) $$

$$ b_b = \frac{1}{2} \Delta_{b,l}^{-1} \begin{bmatrix} (1 + K_b^i)\Delta_{i,L}\Delta_{i,l} \\ (1 - K_b^i)\Delta_{i,L}^{-1}\Delta_{i,l}^{-1} \end{bmatrix}. \qquad (142) $$

Therefore, the transmission coefficient is given by

$$ T = \frac{k_{z,b}}{k_{z,a}} |b_a^T (\mathbf{T}_1 \mathbf{T}_2)^m b_b|^2, \qquad (143) $$

where $\mathbf{T}_2 \equiv e^{i\bar{\varphi}}\bar{\mathbf{T}}_2$ (we have $|e^{\pm i\bar{\varphi}}| = 1$), and the ratio $k_{z,b}/k_{z,a}$ is given by Eq. (116).

### 3.5.3. Comparison with the reference model

Now, we study the effect of having finite size stably stratified interfaces on the transmission coefficient. We maintain the angle $\Theta = \pi/4$, and choose a rotation rate $\Omega = 0.4\bar{N}$. The expression of the transmission coefficient is now given by Eq. (143). The results for different aspect ratios $\varepsilon = 0.01, 0.1$ and $0.5$ for a single step embedded in a stably stratified medium are displayed in Fig. 22. The panels can directly be compared to Fig. 9b, on which is displayed the single step case at the same colatitude and rotation rate, assuming infinitesimally thin interfaces.

As the aspect ratio is increased, a band of enhanced transmission around $\omega = f$ ($\approx 0.57\bar{N}$ here) appears, this being in agreement with Gerkema & Exarchou (2008). We note that this was already observed in the reference model when the transmission coefficient was plotted as a function of the vertical wave number, $k_z$. If we denote by

$$ \omega_{\pm}^{(i)} = \omega_{\pm}(N_i) \qquad (144) $$

the frequency limits between which the incident wave is propagative in the interfaces, there is a frequency range over which the wave is propagative in the convective layers (this happens for $0 < \omega < 2\bar{\Omega}$) and the stably stratified interfaces (this happens for $\omega_{-}^{(i)} < \omega < \omega_{+}^{(i)}$). In our case, the intersection of those two frequency domains is such that this happens for

$$ \omega_{-}^{(i)} < \omega < 2\bar{\Omega}. \qquad (145) $$

This is illustrated in Fig. 23, on which the frequency range above has been overplotted in dashed red lines. Therefore, there is enhanced transmission for waves that are propagative in both convective layers and stably stratified interfaces, as one would expect. The transmission is inhibited for $\omega < \omega_{-}^{(i)}$ for larger aspect ratios $\varepsilon$ – leading to less symmetric transmission about $\omega = f$ than our reference model (which corresponds to $\varepsilon \to 0$) – because the waves are evanescent in the interfaces for $\omega < \omega_{-}^{(i)}$. Thus, as we make $\varepsilon$ larger, the evanescent region becomes larger, which inhibits transmission.

Another difference with the reference model is that there is more transmission near $\omega_+$, and that the transmission peak is shifted to higher wave numbers when the aspect ratio is increased.
<figure><img src="image_26.png" /><figcaption>**Fig. 23.** Transmission coefficient for one step embedded in a stably stratified medium, with $\Omega = 0.4\bar{N}$, $\Theta = \pi/4$ and an aspect ratio $\varepsilon = 0.5$. Transmission is enhanced within the range of frequency $\omega_-^{(i)} < \omega < 2\bar{\Omega}$, for which the waves are propagative both in the convective steps and in the stably stratified interfaces.</figcaption></figure>

# 4. Conclusions and prospects

## 4.1. Summary of results

We have studied the transmission of internal waves through a staircase-like density profile that could be produced by oscillatory double-diffusive convection in giant planet interiors. First, we analysed the free modes of oscillation of such a density staircase in a rotating planet by deriving the dispersion relation that describes them, generalising Belyaev et al. (2015) to consider any latitude (they previously considered the effects of rotation at the pole and the equator). We then analysed the transmission of a wave through a density staircase consisting of one or multiple convective steps by extending Sutherland (2016) to include the complete Coriolis acceleration rather than neglecting its horizontal component, under the so-called ‘traditional approximation’ which is not appropriate when studying convective layers of giant planets (Ogilvie & Lin 2004).

We showed that the transmission of internal (inertial or gravito-inertial) waves is strongly affected by the presence of a density staircase in a frequency- and wavelength-dependent manner. This is true even if these waves are initially pure inertial waves. Large wavelength waves (with wavelengths $\lambda \gg D = \text{md}$, where $\lambda$ is either the vertical or horizontal wavelength and D is the total size of the staircase) are unaffected by the staircase. Low-frequency (inertial or gravito-inertial) waves of any wavelength are perfectly transmitted near the critical latitude $\theta_c = \sin^{-1}(\omega/2\Omega)$ (at which $\omega = f$), confirming a feature also obtained by Gerkema & Exarchou (2008). Otherwise, short-wavelength waves (with $\lambda \sim d \lesssim D$) are only efficiently transmitted if they are resonant with a free mode of the staircase (these are interfacial gravity or short-wavelength inertial modes, corresponding to the roots of the dispersion relation), and if not they are primarily reflected.

The frequency interval around $f = 2\Omega \cos\Theta$ on which transmission remains close to unity widens when rotation is increased (see Fig. 11). This means that for fast rotators, we expect a relatively broad frequency window for which layered

semi-convection has no particular effect upon the propagation of internal waves near the critical latitude. We note that this perfect transmission at $\omega = f$ can only occur at any latitude if $\omega < 2\Omega$, i.e. only for inertial and sub-inertial gravito-inertial waves. We do not expect this to happen for higher frequency gravity waves, since there is no colatitude $\Theta$ for which $\omega = f$.

We also found that non-traditional effects could be of major importance. First, they modify the frequency spectrum of propagation of internal waves, allowing the propagation of sub-inertial waves inside the convective layers of the staircase. Thus, the transmission of pure inertial waves can be strongly affected by the presence of a density staircase. Second, they modify significantly the transmission near the inertial frequency $\omega = f$, with bands of perfect transmission departing from this frequency. These have been physically interpreted as being resonances between pure inertial waves propagating inside a convective step whose wavelength fits exactly inside a step, and the incident wave.

We first modelled the staircase as having infinitesimally thin interfaces, which we also extended to consider interfaces of finite size. An analytical expression of the transmission coefficient has been derived in both set-ups, and its behaviour analysed. In addition to what was already described with the first model (that assumed infinitesimally thin interfaces), it has been found that transmission was significantly enhanced in the frequency range for which internal waves are propagative in both the convective steps and the stably stratified interfaces. As a result, perfect transmission was obtained near $\omega = f$ for any wavelength. Since this model is more realistic, this means that in reality almost any incident IWs or GIWs is perfectly transmitted at a special location, the critical latitude (if its frequency is smaller than $2\Omega$). Otherwise it is strongly affected by the staircase and is only transmitted if it has a large wavelength.

## 4.2. Consequences for the seismology of planets

A case from which we can draw relevant astrophysical conclusions is the case of multiple steps embedded in a convective medium, because gaseous giant planet envelopes are expected to be mostly convective. It is also relevant for the study of a portion of a more vertically extended staircase. In that case when $\omega \neq f$, and the mode is not resonant with a free mode of the staircase, transmission was found to be inhibited for larger wavelength incident waves as we increased the number of steps, being somewhat linearly with $1/m$ (this qualitative behaviour being the same for the other top and bottom layer properties). We find that the cut-off (indexed by the subscript “c”) wavelength below which transmission is less than a given threshold (when not resonant with a mode of the staircase) increases approximately linearly with the number of steps; or, alternatively, is approximately constant when normalised by the size of the staircase $D = \text{md}$. In Fig. 24, we show $\lambda_{z,c}/D$ as a function of the number of steps m, with a rotation rate $\Omega = 0.4\bar{N}$, for different transmission thresholds $T = 0.9, 0.5, 0.1$, at the particular frequency $\omega = 0.2\bar{N}$. We see that, for a given criteria on the transmission coefficient, the cut-off vertical wavelength only depends on the total size of the staircase, D, for a large enough number of steps. For instance, approximately perfect transmission is obtained in that case for $\lambda_z > 10^2 D$.

Assuming that the entire gaseous envelope is semi-convective, Leconte & Chabrier (2012) estimated that the size of the convective steps should lie in a range

$$10^{-9}-10^{-6} \lesssim \frac{d}{H_{\text{p}}} \lesssim 10^{-4}-10^{-2}, \quad (146)$$
<figure><img src="image_27.png" /><figcaption>**Fig. 24.** Vertical cut-off wave-length $λ_{z,c}$ normalised by the size of the staircase $D = md$, as a function of the number of steps $m$, measured for three criteria on the transmission: $T = 0.9$ (green), $T = 0.5$ (red), $T = 0.1$ (blue). These results were obtained in the case of a staircase embedded in a convective medium, with $Ω = 0.4N$ and at the particular frequency $ω = 0.2N$ which was not a resonant mode.</figcaption></figure>

where $H_p$ is the pressure scale-height, that is expected to be of the same order of magnitude as the planet's radius, $R$, in the deep interior of giant planets ($H_p \sim R$). We note that the hypothesis of a fully semi-convective envelope gives us an upper limit on the number of steps. In particular, Vazan et al. (2016) found it to be unlikely because of the efficiency of upward mixing by overturning convection in the upper layers of the envelope. In any case, it provides us with some order of magnitude estimates. We also note that local numerical simulations by for example Radko (2003) show merging of the double-diffusive layers over a finite time, until only one step remains in the local box. This differs from the hypothesis of assuming a large number of steps. However, we stress that those simulations are local and Boussinesq, so that the long-term evolution of double-diffusive convection in spherical geometry with a realistic density stratification is still unknown. Furthermore, the case of the Arctic Ocean on Earth seems to suggest that a staircase with a large number of steps can exist (Ghaemsaidi et al. 2016).

Following our reasoning from Eq. (146), this means that only incident waves coming from above the staircase with a wavelength

$$λ \geq λ_{z,c}, \qquad (147)$$

will be efficiently transmitted to deeper regions, where the cut-off frequency $λ_{z,c}$ is estimated to lie in a range

$$[10^{-7} - 10^{-4}]m \leq \frac{λ_{z,c}}{R} \leq [10^{-2} - 1]m. \qquad (148)$$

Incident waves with $λ_z < λ_{z,c}$ are expected to be primarily reflected at the top of the staircase, unless $ω \approx f$ or the incident wave is resonant with a free mode of the staircase. LC12 estimated that the number of steps in giant planet deep interiors lies in a range

$$[10^2 - 10^4] \leq m \leq [10^6 - 10^9]. \qquad (149)$$

For the bottom estimate $m = 10^2$ steps the estimate given by Eq. (148) yields

$$[10^{-5} - 10^{-2}] \leq \frac{λ_{z,c}}{R} \leq [1 - 10^2]. \qquad (150)$$

We note that the top range of the inequality above does not correspond to any wave that could propagate in giant planet interiors, since their wavelength is larger than the radius of the planet. More restrictive constraints on the layer sizes and their expected numbers would be needed to provide more precise estimates. However, the bottom range for $λ_{z,c}$ is of interest, since tidally excited inertial waves (in the presence of a core) are often thought to be of very short-wavelength (for example Ogilvie & Lin 2004), so their propagation may be strongly affected by a staircase.

On the other hand, Nettelmann et al. (2015) estimate that in the helium rain region, convective layers should have a size $d \sim 0.1-1$ km, corresponding to 10 000–20 000 layers. Following the discussion above and taking a typical giant planet’s radius $R \sim 10^5$ km, this suggests that the cut-off vertical wavelength is approximately given by

$$\frac{λ_{z,c}}{R} \sim [1-10] \qquad (151)$$

in the helium rain region. Because we expect tidal waves to have a vertical wavelength that is small compared to the planet's radius, Eq. (151) suggests that a region of layered semi-convection in the helium rain region (as modelled by Nettelmann et al. 2015) would act as a rigid wall for most internal waves propagating inside the planet's envelope, unless $ω \approx f$ or the incident wave is resonant with a free mode of the staircase.

Waves with shorter wavelengths than $λ_{z,c}$ will be strongly reflected and will not penetrate into deeper regions of the giant planet. This means that a tidally excited wave launched near the surface of a giant planet that has a region of layered semi-convection (for example just outside its core), could see the core's size artificially enlarged because it can be reflected at the top of the staircase region (which could have a radial extent of the order of $10^3-10^4$ km according to Nettelmann et al. 2015), rather than be reflected from the core itself (which would happen in the absence of a density staircase). This could have an impact on the nature of tidal modes that develop inside giant planets containing regions of layered semi-convection, and of course on tidal dissipation.

The (vertical and horizontal) wavelengths are related to the wave's frequency through the dispersion relation. Thus, an inertial wave that propagates in the convective envelope of a giant planet towards the centre has a frequency that is set by Eq. (15). This means that the cut-off wavelengths $λ_{⊥,c}$ and $λ_{z,c}$ define a frequency cut-off through the dispersion relation of pure inertial waves,

$$\omega_c = \left[ \frac{(\tilde{f}λ_{z,c} + fλ_{⊥,c})^2}{λ_{⊥,c}^2 + λ_{z,c}^2} \right]^{1/2}. \qquad (152)$$

In the presence of layered semi-convection, potentially observable modes of different frequencies will see a different size of the core. This observation could be a signature of layered semi-convection just outside the rocky/icy core of giant planets.

On the other hand, if density staircases are also present in the helium rain region, the spherical shell contained between that and the core could act as a trap for pure inertial waves in a certain frequency range, modifying the nature of the modes that can be sustained in that particular region. However, this could not be observed by seismology of the planet's surface.

Finally, we have seen that new modes can potentially be observed when invoking layered semi-convection compared to a fully convective envelope. These are the free modes of the staircase, first given by BQF15 and extended here to the case with
<figure><img src="image_28.png" /><figcaption>**Fig. 25.** Sketch presenting an illustrative summary of our motivations, results and future work. North-west quarter: how density staircases are possibly created (see Sect. 1). South-west quarter: main properties of the transmission of internal waves (see Sects. 3 and 4.2). South-east quarter: additional modes that could potentially be observed thanks to the seismology of planets (see Sect. 4.2). North-east quarter: the question we address in our second paper.</figcaption></figure>

rotation for any latitude through Eq. (101). This includes the gravito-inertial modes with a sufficiently large radial wavelength that the staircase is seen as a stably stratified medium characterised by the mean buoyancy frequency $\bar{N}$.

### 4.3. Perspectives

Our results help us to understand the impact of layered semi-convection on the propagation of internal waves at a qualitative level. In order to make quantitative predictions, one should extend this study to spherical geometry relevant for planets and stars. In particular, it could be interesting to analyse probabilities of excitation of free modes of the staircase. This could allow us to determine for example what fraction of tidal waves are perfectly transmitted through the staircase region, but this is beyond the scope of our paper.

In addition, we have considered a density staircase established by the double-diffusive instability, as described by e.g. Leconte & Chabrier (2012). We note, however, that it would be important in the near future to take into account the impact of rotation (Moll & Garaud 2017), differential rotation (Worthem et al. 1983) and magnetic fields on the formation of such a structure. Waves will also be affected by differential rotation (Mathis 2009; Mirouh et al. 2016) and magnetic fields (for example Mathis & de Brye 2011; Wei 2016). Moreover, nonlinear effects such as the feedback of waves on the background should be taken into account.

Finally, a synthetic sketch of our results is shown in Fig. 25. Based on this first study, we will explore in our second paper the tidal response and dissipation of a giant planet or star containing

a staircase-like density profile resulting from the presence of layered semi-convection.

**Acknowledgements.** We thank the referee for a constructive and helpful report. Q.A. was supported by ENS Paris-Saclay and CEA. Q.A. and S.M. acknowledge funding by the European Research Council through ERC SPIRE grant 647383. A.J.B. was supported by the Leverhulme Trust through the award of an Early Career Fellowship. Q.A. was also partly supported by this Leverhulme Trust award.

### References

Auclair-Desrotour, P., Le Poncin-Lafitte, C., & Mathis, S. 2014, A&A, 561, L7

Auclair Desrotour, P., Mathis, S., & Le Poncin-Lafitte, C. 2015, A&A, 581, A118

Baraffe, I., Chabrier, G., Barman, T. S., et al. 2005, A&A, 436, L47

Belyaev, M. A., Quataert, E., & Fuller, J. 2015, MNRAS, 452, 2700

Bodenheimer, P., Lin, D. N. C., & Mardling, R. A. 2001, ApJ, 548, 466

Bolmont, E., & Mathis, S. 2016, Celest. Mech. Dyn. Astron., 126, 275

Chabrier, G., & Baraffe, I. 2007, ApJ, 661, L81

Efroimsky, M., & Lainey, V. 2007, J. Geophys. Res. (Planets), 112, E12003

Fortney, J. J., Ikoma, M., Nettelmann, N., Guillot, T., & Marley, M. S. 2011, ApJ, 729, 32

Friedlander, S. 1987, Geophys. J., 89, 637

Fuller, J. 2014, Icarus, 242, 283

Garaud, P. 2013, in EAS Publ. Ser. 63, eds. G. Alecian, Y. Lebreton, O. Richard, & G. Vauclair, 285

Gaulme, P., Schmider, F.-X., Gay, J., Guillot, T., & Jacob, C. 2011, A&A, 531, A104

Gerkema, T., & Exarchou, E. 2008, J. Marine Res., 66, 617

Gerkema, T., & Shrira, V. I. 2005, J. Fluid Mech., 529, 195

Ghaemsaidi, S. J., Dosser, H. V., Rainville, L., & Peacock, T. 2016, J. Fluid Mech., 789, 617

González-Cataldo, F., Wilson, H. F., & Militzer, B. 2014, ApJ, 787, 79

Goodman, J., & Lackner, C. 2009, ApJ, 696, 2054

Guillot, T., & Showman, A. P. 2002, A&A, 385, 156

Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, The interior of Jupiter, eds. F. Bagenal, T. E. Dowling, & W. B. McKinnon, 35

Kurokawa, H., & Inutsuka, S.-I. 2015, ApJ, 815, 78

Lainey, V., Arlot, J.-E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957

Lainey, V., Karatekin, Ö., Desmars, J., et al. 2012, ApJ, 752, 14

Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286

Laskar, J., Boué, G., & Correia, A. C. M. 2012, A&A, 538, A105

Leconte, J., & Chabrier, G. 2012, A&A, 540, A20

Leconte, J., & Chabrier, G. 2013, Nature Geoscience, 6, 347

Ledoux, P. 1947, ApJ, 105, 305

Mathis, S. 2009, A&A, 506, 811

Mathis, S., & de Brye, N. 2011, A&A, 526, A65

Mathis, S., Neiner, C., & Tran Minh, N. 2014, A&A, 565, A47

Mathis, S., Auclair-Desrotour, P., Guenel, M., Gallet, F., & Le Poncin-Lafitte, C. 2016, A&A, 592, A33

Mayor, M., & Queloz, D. 1995, Nature, 378, 355

Mazevet, S., Tsuchiya, T., Taniuchi, T., Benuzzi-Mounaix, A., & Guyot, F. 2015, Phys. Rev. B, 92, 014105

Militzer, B., Soubiran, F., Wahl, S. M., & Hubbard, W. 2016, J. Geophys. Res. (Planets), 121, 1552

Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & Wood, T. S. 2012, ApJ, 750, 61

Mirouh, G. M., Baruteau, C., Rieutord, M., & Ballot, J. 2016, J. Fluid Mech., 800, 213

Moll, R., & Garaud, P. 2017, ApJ, 834, 44

Moll, R., Garaud, P., & Stellmach, S. 2016, ApJ, 823, 33

Nettelmann, N., Fortney, J. J., Moore, K., & Mankovich, C. 2015, MNRAS, 447, 3422

Ogilvie, G. I., & Lesur, G. 2012, MNRAS, 422, 1975

Ogilvie, G. I., & Lin, D. N. C. 2004, ApJ, 610, 477

Pollack, J. B., Grossman, A. S., Moore, R., & Graboske, Jr., H. C. 1977, Icarus, 30, 111

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62

Radko, T. 2003, J. Fluid Mech., 497, 365

Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2011, ApJ, 731, 66

Sakashita, S., & Hayashi, C. 1959, Progr. Theor. Phys., 22, 830

Salpeter, E. E. 1973, ApJ, 181, L83

Schwarzschild, M., & Härm, R. 1958, ApJ, 128, 348

Stellmach, S., Traxler, A., Garaud, P., Brumell, N., & Radko, T. 2011, J. Fluid Mech., 677, 554

Stern, M. 1960, Tellus A, 12, 172
Stevenson, D. J. 1975, Phys. Rev. B, 12, 3999

Stevenson, D. J. 1985, Icarus, 62, 4

Stevenson, D. J., & Salpeter, E. E. 1977, ApJS, 35, 239

Strugarek, A. 2016, ApJ, 833, 140

Sutherland, B. R. 2016, Phys. Rev. Fluids, 1, 013701

Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016, ApJ, 829, 118

Wahl, S. M., Wilson, H. F., & Militzer, B. 2013, ApJ, 773, 95

Wei, X. 2016, ApJ, 828, 30

Wilson, H. F., & Militzer, B. 2012a, Phys. Rev. Lett., 108, 111101

Wilson, H. F., & Militzer, B. 2012b, ApJ, 745, 54

Wood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157

Worthem, S., Mollo-Christensen, E., & Ostapoff, F. 1983, J. Fluid Mech., 133, 297

Zahn, J. P. 1966, Annales d'Astrophysique, 29, 313

Zahn, J.-P. 1989, A&A, 220, 112

## Appendix A: Derivation of the Poincaré equation

The linearised system we consider, when we adopt the Boussinesq approximation, is given by Eqs. (8)–(12). The aim of this appendix is to demonstrate the Poincaré equation, used in Sect. 2.3.

First, taking the combination $\partial_y(10) - \partial_z(9)$ gives

$$ \frac{\partial}{\partial t} \left( \frac{\partial w}{\partial y} - \frac{\partial v}{\partial z} \right) - (f \cdot \nabla)u = \frac{\partial b}{\partial y}, \quad (A.1) $$

where $f = (0, \tilde{f}, f)$. Then, taking the combination $\partial_z(8) - \partial_x(10)$ and using Eq. (11) gives

$$ \frac{\partial}{\partial t} \left( \frac{\partial u}{\partial z} - \frac{\partial w}{\partial x} \right) - (f \cdot \nabla)u = - \frac{\partial b}{\partial x}. \quad (A.2) $$

Then, taking the combination $\partial_x(9) - \partial_y(8)$ and using Eq. (11) gives

$$ \frac{\partial}{\partial t} \left( \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right) - (f \cdot \nabla)w = 0. \quad (A.3) $$

Then, by taking the combination $\partial_t(\partial_y(A.1) - \partial_x(A.2))$ and using Eqs. (11), (12) and (A.3), we finally obtain an equation on the vertical component of the velocity $w$

$$ \frac{\partial^2}{\partial t^2} \nabla^2 w + (f \cdot \nabla)^2 w + [N^2 \nabla_{\perp}^2] w = 0, \quad (A.4) $$

where $\nabla_{\perp}^2 \equiv \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$.