Elliptical instability is an instability of elliptical streamlines, which can be excited by large-scale tidal flows in rotating fluid bodies and excites inertial waves if the dimensionless tidal amplitude (*ε*) is sufficiently large. It operates in convection zones, but its interactions with turbulent convection have not been studied in this context. We perform an extensive suite of Cartesian hydrodynamical simulations in wide boxes to explore the interactions of elliptical instability and Rayleigh–Bénard convection. We find that geostrophic vortices generated by the elliptical instability dominate the flow, with energies far exceeding those of the inertial waves. Furthermore, we find that the elliptical instability can operate with convection, but it is suppressed for sufficiently strong convection, primarily by convectively driven large-scale vortices. We examine the flow in Fourier space, allowing us to determine the energetically dominant frequencies and wavenumbers. We find that power primarily concentrates in geostrophic vortices, in convectively unstable wavenumbers, and along the inertial wave dispersion relation, even in non-elliptically deformed convective flows. Examining linear growth rates on a convective background, we find that convective large-scale vortices suppress the elliptical instability in the same way as the geostrophic vortices created by the elliptical instability itself. Finally, convective motions act as an effective viscosity on large-scale tidal flows, providing a sustained energy transfer (scaling as $\epsilon 2$). Furthermore, we find that the energy transfer resulting from bursts of elliptical instability, when it operates, is consistent with the $\epsilon 3$ scaling found in prior work.

## I. INTRODUCTION

The elliptical instability arises when elliptically deformed streamlines excite pairs of inertial waves through parametric resonances.^{1–3} As long as retarding processes such as viscous damping can be overcome, an arbitrary small elliptical deformation can yield instability. The resulting inertial waves couple with the deformation,^{1} leading to the exponential growth of their amplitudes. This mechanism is in essence a triadic resonance interaction in which waves extract energy from the elliptical flow.

Its nonlinear evolution has been studied extensively. As the linear instability saturates, the inertial waves appear to collapse to rotating turbulence, which typically dissipates over time, leading to the flow becoming unstable again.^{4–10} This collapse to turbulence either occurs via weak inertial wave “turbulence,”^{4,6,9,10} or rotating turbulence involving large-scale geostrophic vortices or zonal flows.^{5,7–11} The inertial wave turbulence (involving a sea of weakly interacting inertial waves) may occur when the forcing amplitude is weak,^{10,12} or when geostrophic modes are suppressed, either by artificial frictional damping^{9} or via an external process such as the imposition of a magnetic field.^{6}

In recent years, elliptical instability also finds application as a tidal dissipation mechanism in stars and planets,^{5,8,13–18} extracting energy from a tidal flow. Tidal flows in stars or planets are usually split up into a large-scale equilibrium or non-wave-like tide, and a dynamical or wave-like tide.^{19,20} The equilibrium tide is the quasi-hydrostatic fluid bulge rotating around the body,^{19} while the dynamical tide consists of waves generated by resonant tidal forcing. The equilibrium tide is thought to be dissipated through its interaction with turbulence, usually of a convective nature,^{21,22} or by its own fluid instabilities,^{5,8,13–15,17,23} among which is the elliptical instability. This however requires careful consideration of the dynamics of the elliptical instability, particularly the properties of the turbulence which is expected, as well as its interaction with other processes in the system, such as magnetic fields or (stable or unstable) stratification.^{24}

The equilibrium tide deforms a body (body 1) into an ellipsoidal shape that follows an orbiting companion (body 2), and its deformation is represented by the ellipticity or tidal amplitude parameter,

with *m*_{1} and *m*_{2} being the masses of bodies 1 and 2 (e.g., a planet and its host star), *R*_{1} is the radius of body 1, and *a* is the orbital separation (semi-major axis). In an asynchronously rotating planet or star, the equilibrium tide in the frame rotating with the tidal bulge is an elliptical flow inside the planet. The rotation rate of this flow is the difference between the spin of the planet Ω and the orbital rotation rate *n* and is denoted by $\gamma \u2261\Omega \u2212n$. In this work, we model a small patch of an equilibrium tidal flow, which is treated as a background flow $U\u03030$ in the bulge frame (rotating at the rate *n* about the axis of rotation of the planet) following,^{5} such that:

where **x** represents the position vector from the center of the planet. In the frame rotating with the planet at the rate Ω, it can be written alternatively as the flow

where **x** now represents the position vector from the center of the planet in the frame rotating with the planet. This description represents the exact equilibrium tide of a uniformly rotating incompressible fluid body perturbed by an orbiting companion^{8,25} but approximates the main features of the equilibrium tide in more realistic models.^{20,26} We choose to work in the frame rotating with the planet at the rate Ω in this study. Larger deformations *ε* result in faster growth of the waves which means that they can overcome stronger viscous damping by either the viscosity of the fluid or by a turbulent viscosity.

The elliptical instability has been studied previously in simulations using a local Cartesian box model both with^{6} and without^{5} weak magnetic fields. The earlier study found that elliptical instability leads to bursty behavior, involving the interaction of instability-generated inertial waves with geostrophic columnar vortical flows produced by their nonlinear interactions. Irregular cyclic “predator**–**prey behavior” was obtained in which the elliptical instability first excited inertial waves, these interacted nonlinearly to produce vortices that inhibited further growth of waves until these vortices were damped sufficiently by viscosity, thereby enabling further growth of the waves. Similar behavior features in global hydrodynamical simulations of the elliptical instability,^{8} where zonal flows take the place of columnar vortices in the predator**–**prey dynamics. Upon taking magnetic fields into account in the local model, the behavior changed from bursts to a sustained energy input into the flow, as magnetic fields served to break up or prevent formation of strong vortices.^{6} Similar sustained behavior is observed if the vortices are damped by an artificial frictional force mimicking Ekman friction on no-slip boundaries.^{9}

These prior studies analyzed the elliptical instability in the convective regions of planetary (or stellar) interiors but did not incorporate convection explicitly (except perhaps by motivating a choice of viscosity if this is due to turbulence). However, convection can potentially interact with the elliptical instability in a number of ways. First, we might imagine that smaller-scale turbulent convective eddies could act like an effective viscosity in damping larger-scale inertial waves, thereby inhibiting or reducing the growth of the elliptical instability. Second, convection under the influence of rotation is known to generate mean flows (zonal flows or vortices), and these flows could also interact with those generated by the elliptical instability. In this study we wish to address the following questions: Can the elliptical instability operate in a turbulent convective background? How do convectively driven flows interact with the elliptical instability and modify (tidal) energy transfer rates?

The interaction of the elliptical instability with convection has been studied within linear theory,^{2,27} experimentally in cylindrical containers,^{28} and using idealized laminar global simulations in a triaxial ellipsoid.^{15} These studies illustrate that the elliptical instability can modify heat transport, though they did not focus on the dissipation of tidal flows, which is our primary focus here. The dimensionless heat transport is usually represented by the Nusselt number (Nu), which is a measure of the ratio of the total heat flux to the conductive flux (i.e., with no transport by fluid motions), as a function of the Rayleigh number (Ra), the dimensionless ratio of buoyancy driving to viscous and thermal damping, which measures the strength of convective driving. The Nusselt number was observed to be increased by the elliptical instability for Rayleigh numbers close to and below the value required for onset of convection. It was also observed to be larger than one even with stable stratification (Ra < 0), indicating that the elliptical instability can contribute to heat transport in this regime also.

In this paper, we build upon the local box model—which represents a small patch of a planet or star (see Fig. 1)—in the study by Barker and Lithwick,^{5,6} Le Reun *et al.*^{9} to study the interaction of the elliptical instability with convection with a focus on the resulting tidal dissipation. Our local model allows for higher resolution studies, which, in turn, allows us to reach more turbulent regimes than, e.g., Cébron, Maubert, and Le Bars^{15} and Lavorel and Le Bars.^{28} One of our aims is to study the behavior of the elliptical instability and to see if introducing convection could also lead to sustained energy injection in a similar fashion to a magnetic field (or frictional damping). As a secondary goal, we are interested in studying the modifications to heat transport by the elliptical instability, both in the weakly driven regime of convection and in stably stratified layers like those that may exist in giant planet interiors.

In Sec. II, we will discuss the linear properties of the elliptical instability and describe the model used, and in Sec. III, we describe the results obtained from our parameter survey in a qualitative manner. We investigate the sustained energy injection and analyze frequency and wavenumber spectra for the flow in Sec. IV. Then, we briefly discuss the scaling of the energy transfer from the background flow with *ε* in Sec. V. We discuss and conclude in Sec. VI.

## II. MODEL SETUP

### A. The elliptical instability

The linear properties of the elliptical instability have been reviewed by Kerswell.^{2} This instability operates when two inertial waves have frequencies that approximately add up to the tidal frequency $2\gamma $ (see Sec. I). In the short wavelength limit, this occurs for two waves with frequencies $\omega =\xb1\gamma $, which must each satisfy the dispersion relation for inertial waves, $\omega =\xb12\Omega kz/k$, where $kz/k=cos\u2009\theta $. The elliptical instability grows at a rate proportional to $\epsilon \gamma $.

Since we are investigating a small patch of a planet, we assume the tidal flow can be modeled locally as an unbounded strained vortex in the bulge frame.^{2} This approach yields a growth rate which depends on the angle the inertial waves make with respect to the rotation axis and the strain direction in the horizontal plane. For illustration, when the tidal bulge is stationary (*n *=* *0), $\gamma =\Omega $, and thus $kz/k=\xb11/2$ for the most unstable modes.^{2} Furthermore, the fastest growing waves have a phase aligned with the strain direction, by an angle of $\xb1\pi /4$ with respect to the elliptical deformation in the equatorial plane (i.e., plane containing the vortical flow).

The effects of viscosity, detuning, convection, and rotation of the elliptical bulge are also reviewed by Kerswell.^{2} Viscosity reduces the growth rate according to: $\sigma \u2212\nu |k|2$, where *σ* is the inviscid growth rate, *ν* is the viscosity, and **k** is the wavevector of the fastest growing mode. Detuning is a reduction of the growth rate as a result of the wave not satisfying the resonance conditions exactly, which also reduces the maximum growth rate. The rotation around the companion, and thus, the rotation of the elliptical bulge, modifies the growth rate depending on the rotation speed (*n*). The growth rate is decreased for most values of^{5} *n*, and it cannot operate in the interval $n=[\u22123/2\gamma ,\u22121/2\gamma ]$. In this interval, no inertial waves exist that satisfy the dispersion relation defined by

In the interval $n=[\u22121/2\gamma ,0]$, the growth rate is increased, though everywhere else it is decreased, over the case with *n *=* *0. The linear growth rate of an inviscid fluid at small *ε* without convection is given by^{2,29}

If an (un)stable stratification (aligned with the rotation axis) is present the dispersion relation is modified, as well as the above equation. The stratification introduced into the equation is represented by the Brunt–Väisälä (or buoyancy) frequency *N*. The modified dispersion relation is

Both the effects of unstable stratification (negative *N*^{2}) and stable stratification (positive *N*^{2}) can be computed using this equation. We observe that stable stratification typically inhibits elliptical instability, but that convection typically enhances growth.

For clarity of presentation $\gamma =\Omega $ is chosen, resulting in *n *=* *0, i.e., strictly representing the unphysical case where there is no rotation of the bulge. The body in question is not rotating around its companion which causes the tidal effects. However, it turns out that for simulations, the only linear effects of choosing a different value of Ω, and therefore a non-zero value of *n*, would be to modify the growth rate of the elliptical instability^{5} as well as the wavenumber of the most unstable mode.^{6}

### B. Governing equations and setup of the simulations

We model the convective instability using rotating Rayleigh–Bénard convection (RRBC) and the Boussinesq approximation. RRBC is chosen as it is the simplest model of rotating convection^{30} which allows us to study its interaction with elliptical instability (an even simpler model is “homogeneous convection” with periodic boundaries in the vertical, but this has unphysical nonlinear behavior). The Boussinesq approximation is justified when studying small-scale convective (and wavelike) flows, which satisfies the required conditions that flows are much slower than the sound speed, $u\u226acs$, and the vertical size of the domain *d* is much smaller than a pressure or density scale height, $d\u226aHp$.^{31} However, this neglects variations of the properties of a planet, i.e., density and temperature, and of course, any large-scale circulations cannot be modeled using this approximation.

The rotation axis is aligned with the *z*-direction, as is the temperature gradient, as indicated in Fig. 1. The box in the current setup thus represents a polar region, because the local rotation and gravity vectors are either aligned or anti-aligned (depending on the sign of Ω). The conduction state temperature profile *T*(*z*) between the hot plate at the bottom and the cold plate at the top, about which we will perturb, is

where *g* is the local gravitational acceleration (assumed constant) and *α* is the (constant) thermal expansion coefficient. Without loss of generality *T*_{0} is set to zero, so the temperature at the bottom (hot plate, typically) is $T(z=0)=0$, while the temperature at the top is $T(z=d)=N2/(\alpha g)$, such that the temperature drop is $\Delta T=\u2212N2/\alpha g$.

We non-dimensionalize by scaling lengths with the vertical domain size *d* (distance between the plates), scaling time using the corresponding thermal timescale $d2/\kappa $, thus scaling velocities with $\kappa /d$, pressures with $\rho 0\kappa 2/d2$, and finally scaling temperature with $T=\Delta T\theta $ (i.e., by the temperature difference between the plates). The governing equations for the dimensionless velocity and temperature perturbations **u** and *θ* to the background flow $U0$ and conduction state temperature *T*(*z*) in the Boussinesq approximation, in a frame rotating at the rate Ω about *z*, are then

where

with $u=(ux,uy,uz)$ and *p* being the perturbation to the pressure. The non-dimensional parameters describing the convection are the Rayleigh number,

where *ν* and *κ* are the constant kinematic viscosity and thermal diffusivity, the Ekman number (ratio of viscous to Coriolis terms)

and the Prandtl number $Pr=\nu /\kappa $. Note that we can relate the dimensional squared buoyancy frequency $N2=\u2212Ra\u2009Pr\u2009\kappa 2/(\alpha gd4)$, so that when $Pr=1$, the dimensionless value (in thermal time units) is $N2=\u2212Ra$. The tidal background flow also introduces the dimensionless ellipticity *ε* and the frequency *γ* in our chosen units.

Our simulations are executed in a small Cartesian box of dimensionless size $[Lx,Ly,1]$ with $Lx=Ly=L$. The boundary conditions in the horizontal directions are periodic, while in the vertical direction, they are impermeable, $uz(z=0)=uz(z=1)=0$, and stress-free, $\u2202zux(z=0)=\u2202zux(z=1)=\u2202zuy(z=0)=\u2202zuy(z=1)=0$. Stress-free boundary conditions are chosen both for numerical convenience and because they are probably more physically relevant than no-slip boundary conditions for modeling convection in the deep interior of a planet, far from boundaries. The convection in our box, thus, represents a single convection cell in the vertical. Boundary conditions in the vertical for the temperature perturbation are assumed to be perfectly conducting, $\theta (z=0)=\theta (z=1)=0$.

In the derivation of the elliptical instability, the choice is often made to work with so-called shearing waves.^{1,2} Shearing waves have time-dependent wavevectors, which allows us to account for the stretching and rotation of waves due to a background flow, such as the equilibrium tide in our simulation. A single shearing wave (sometimes also referred to as a Kelvin wave,^{2} although strictly different to a coastal Kelvin wave) is represented as

where $k\u0307\u22a5=\u2212ATk\u22a5$ and $k\u22a5=(kx,ky,0)$, and we use a basis of these waves in our simulations following Barker and Lithwick,^{5} except that we use a sine–cosine decomposition in *z* similar to Duguid, Barker, and Jones.^{32}

The simulations are executed using the Snoopy code.^{33} The Snoopy code implements a Fourier pseudo-spectral method using FFTW3 in a Cartesian box. We use a sine–cosine decomposition in *z*, as in Eq. (15), and shearing waves (i.e., Fourier modes) in *x* and *y*. A third-order Runge-Kutta scheme is used for the time-stepping, together with a Courant–Friedrichs–Lewy (CFL) safety factor to ensure the timesteps are small enough to accurately capture non-linear effects, usually set to 1.5 (which is smaller than the stability limit of $3$). The anti-aliasing in the code uses the standard 2/3 rule.^{34} A variety of different Rayleigh numbers were analyzed using the simulations. In addition, some simulations were performed with Rayleigh numbers in the stably stratified regime, i.e., with $Ra<0$. The values of the Rayleigh number are typically reported as $Ra/Rac$ for clarity, where $Rac$ is the onset Rayleigh number (determined numerically), and the range of this ratio studied at $Ek=5\xd710\u22125.5$ is from 2 to 20 and from –10 to 0.8, in the convectively unstable and stable regimes, respectively. We vary *ε* from 0.01 to 0.20.

In simulations of RRBC in a local box model, large-scale vortex (LSV) structures emerge in the flow when rotation dominates.^{35–37} This LSV emerges with our chosen boundary conditions (more details below) and grows to the size of the box in the horizontal. One of the effects of the LSV is to reduce heat transport, as the vertical motions are suppressed by such a vortex.^{35} An additional reason to study convective LSVs is that the elliptical instability can be suppressed by the presence of strong vortices.^{5} The convective LSV might suppress the elliptical instability as well, potentially preventing it from operating efficiently. Since our flow is likely to be rotationally dominated due to our choice of Ekman number and computationally feasible Rayleigh numbers, a horizontal box size was chosen that would capture this LSV. The question remains, however, what effect changing the aspect ratio of the box would have on the effects presented in this report, as the aspect ratio *L*/*d* (the ratio of horizontal length of the box to its vertical length) influences the ratio of the vertical to total kinetic energy.^{35}

### C. Energetic analysis of simulations

To analyze the energy of the flow, we derive a kinetic energy equation by taking the scalar product with **u** of Eq. (9) and then averaging over the box. We define our averaging operation on a quantity *X* as $\u27e8X\u27e9=1L2d\u222bVX\u2009dV$. We obtain

where we have defined the mean kinetic energy,

the mean viscous dissipation rate,

and the energy injection rate (more generally, energy transfer rate) from the tidal to convective flows (or vice versa),

To obtain an equation for the thermal (potential) energy, we multiply Eq. (11) by –RaPr*θ* and average over the box to obtain

where we have defined the mean thermal energy as

and the thermal dissipation rate as

The total energy is $E=K+P$, which, thus, obeys

In a statistically steady state, it is expected that the (time-averaged value of the) energy injected balances the total dissipation, i.e., $I\u2248D\u2261D\nu +D\kappa $ (on average). This sum of the two dissipation rates then represents the tidal energy dissipation rate resulting from the tidal energy injected. Therefore, to interpret the tidal energy dissipation rate, we examine the tidal energy injection rate *I*.

Arguments to describe scaling laws for the dissipation due to the elliptical instability were first described by picturing the instability saturation as involving a single most unstable mode whose amplitude saturates when its growth rate (*σ*) balances its nonlinear cascade rate.^{5} Thus, if the most important mode of the elliptical instability satisfies $\sigma \u223cku$, where *k* is its wavenumber magnitude and *u* is its velocity amplitude, then we find $u\u223c\epsilon \gamma /k$. The total dissipation rate *D*, therefore, scales as $D\u223cu2\sigma \u223c\epsilon 3\gamma 3/k2$. Thus, in such a statistically steady state, the dissipation and energy injection rate are expected to scale as

and this is consistent with some local and global simulations^{5,8} as well as the scaling found for related instabilities like the precessional instability.^{38,39} We are interested in exploring whether convection could lead to a different result and potentially reduce this steep *ε* scaling.

Since we know both the elliptical instability^{5} and convection^{35} in isolation can produce geostrophic flows such as vortices, we introduce further diagnostics to analyze these flows and the roles they play. To do this, we decompose the total energy injection from the background flow into

where we have defined $I2D=\u2212\u27e8u2DAu2D\u27e9$ and $I3D=\u2212\u27e8u3DAu3D\u27e9$. $I2D$ and $u2D$ are defined to include all (geostrophic) modes where the wavevector has only non-vanishing *x* and *y* components, with *k _{z}* = 0, and $I3D$ and $u3D$ includes all the modes with $kz\u22600$. Thus, we have decomposed the total energy injection rate into energy injection into the barotropic (

*k*= 0) and baroclinic ($kz\u22600$) flow. A pure inertial wave with

_{z}*k*= 0 would have zero frequency, while in convectively unstable simulations, which is the main focus of this work, no gravity waves exist which could have $\omega \u22600$ even when

_{z}*k*= 0, and therefore, one can crudely think of this decomposition as one into geostrophic vortex modes ($I2D$) and waves ($I3D$). We have found that the time-averaged energy input into the vortical motions $I2D$ is approximately zero,

_{z}^{5}but that the input into the waves $I3D$ is on average non-zero (which it must be when the elliptical instability operates) and clearly demonstrates any bursty behavior observed. Based on this observation, only results derived from $I3D$ will be plotted in this paper. The total kinetic energy

*K*is also split up into a 2D and 3D component in a similar manner by defining $K2D$ and $K3D$, so as to allow us to determine which components contribute the most and dominate the flow.

To further analyze the energy transfer rates *I* and $I3D$, we also convert these to an effective viscosity *ν _{eff}* and $\nu eff,3D$, respectively. The effective viscosity represents the energy dissipation that would result from a constant kinematic viscosity with the value $\nu =\nu eff$, and this quantity allows us to interpret the value of

*I*. In particular, this is a useful comparison to quantify the rate at which turbulent convection could damp our tidal flow, if this interaction behaves like a turbulent viscosity. To define the effective viscosity, we equate the work done by the tidal flow on the convective flow with the viscous dissipation rate of the tidal flow, assuming that this is due to a constant kinematic viscosity

*ν*, following Duguid, Barker, and Jones,

_{eff}^{32}Goodman and Oh,

^{40}Ogilvie and Lesur,

^{41}Vidal and Barker.

^{42}We note that

We define the strain rate tensor for the tidal flow as $eij0\u226112(\u2202iU0,j+\u2202jU0,i)$, such that the rate at which energy is dissipated is given by

The effective viscosity is then defined by

The injection terms represent energy being transferred from the background flow to the perturbations or vice versa. This, by definition, impacts the energy in these flows. The evolution of the tidal flow $U0$, however, is not explicitly accounted for in our model, which we treat as a fixed (but time-dependent) background flow. The idea is that it has much larger energy than the perturbations, so it is treated as an infinite reservoir in our simulations. These results therefore give us only a snapshot at a certain point in time of the evolution of the system, which is reasonable because tidal evolutionary processes usually occur slowly relative to convective or rotational timescales.

Finally, we compute the vertical heat transport in our simulations, represented by the Nusselt number, which we define as

This gives the ratio of the total heat flux to the conductive flux and would take the value one in the absence of flows (i.e., heat transported purely by conduction).

### D. Linear growth rates and numerical validation

The predicted growth rate of the elliptical instability is given in Eq. (7), however, this was derived for an unbounded flow in *z*. Since we have adopted the RRBC setup with impermeable walls in *z*, we must determine how this affects the growth rate of instability, although we expect it remains unchanged (and we have also demonstrated this analytically, though we omit this derivation). Hence, we performed multiple test simulations analyzing the linear growth rates of both the elliptical and convective instabilities. We initialized them with random noise and the non-linearities were switched off, thus leading to a continuous exponential growth, allowing for easy extraction of the growth rate. Fits were performed to the mean kinetic energy on a log-scale to determine *σ*, by noting that if velocity components grow as $exp\u2009(\sigma t)$, then $K\u221d\u2009exp\u2009(2\sigma t)$.

The results, along with the theoretical growth rate predictions for both instabilities, are plotted in Fig. 2. The top panel shows the growth rate of the elliptical instability as a function of *ε* (when *n *=* *0), where we have adopted a time unit $\gamma \u22121$, equivalent to using *γ* = 1, for the purposes of this figure. It can be concluded that the modeling of the inviscid growth rate is approximately correct. We also accounted for viscous damping by including the viscous decay rate so that the total growth rate is $\sigma \u2212\nu k2$ (where *k* is the wavevector magnitude of the mode), which is in excellent agreement with our simulations. We obtain values that are very slightly smaller than the theoretical prediction though, even when taking into account viscosity, which is likely due to the detuning effect discussed previously (i.e., that the mode does not precisely satisfy $kz/k=1/2$). In numerical simulations, this detuning arises because of the finite number of grid points, which, in addition to a chosen aspect ratio, prohibits the waves from precisely satisfying the aforementioned condition. However, because this condition does not stipulate the size of the wavenumbers, instead stipulating their ratio, and as such direction, the fastest growing mode that dominates the volume-averaged energy will be as large-scale as possible while still adequately satisfying the resonance condition in order to reduce the viscosity correction. This means that it should be unaffected by resolution, instead being controlled by the aspect ratio of the box. The effect on the growth rate of this detuning, which here corresponds to the difference between the markers and the viscosity corrected growth rate, at this aspect ratio and $\gamma =\Omega =1$, corresponding to the top panel of Fig. 2 and the other simulations in this work, is determined numerically to be $\u22480.002$.

The growth rate of the elliptical instability as a function of Ω, keeping *γ* = 1 is shown in the middle panel of Fig. 2 and also follows the theoretical prediction well but is again slightly lower for the same reasons. The growth rate of the convective instability for $Ek=5\xd710\u22125.5$ is shown in the bottom panel of Fig. 2, now using thermal time units, and this is also in very good agreement with the linear convective growth rate for the fastest growing mode (which scales as $Ek\u22121/3$ as obtained from solving the relevant cubic dispersion relation numerically) as expected. We can, therefore, be confident that both instabilities have been captured correctly.

## III. NUMERICAL RESULTS

### A. Qualitative analysis of illustrative simulations

We begin our discussion of simulation results by presenting the *z*-averaged vertical vorticity $\u27e8\omega z\u27e9z$ of the flow in Fig. 3, taken from a snapshot at *t *=* *0.08 in a simulation in which only the elliptical instability and its associated nonlinear dynamics are present with $\epsilon =0.1,\u2009Ra=0$ (left) and a simulation in which both the elliptical instability and the convective instability, as well as both their associated nonlinear dynamics, are present with $\epsilon =0.1,\u2009Ra=6Rac$ (right). This time is after the initial saturation of instability, in which an LSV has formed in the flow in both cases. These LSVs are an important feature produced by nonlinear evolution of both elliptical instability and convection. Note that with our chosen aspect ratio and Ekman number, the convective LSV emerges when $Ra\u22733Rac$. Thus, they emerge at similar values of the Rayleigh number as the convective LSV in Guervilly, Hughes, and Jones,^{35} which can be readily seen upon taking into account the factor of $\u22488.7$ included in the definition of $Rac$, as a result of which the parameter $Ra\u0303=RaEk4/3\u227320$ from their work^{35} translates to $Ra\u227320Ek\u22124/3\u22482.5Rac$. The vortices are centered in these images for clarity. The nonlinear evolution of the elliptical instability in the left panel creates a cyclonic vortex and a smaller and weaker anticyclonic vortex (at the corners, noting the periodic boundaries). In the right panel, the convection dominates the flow, and as a result, the convective LSV dominates and results in a single cyclonic vortex, with a primarily anticyclonic background. The vorticity in the center of this vortex is larger than the vorticity in the center of the one in the left panel. In the right panel, small-scale convective eddies are present throughout the box, making the flow appear much noisier compared to the elliptical instability in isolation.

Barker and Lithwick^{5} found that these vortices are produced by the nonlinear saturation of the elliptical instability and play a key role in the predator–prey behavior of the inertial waves and vortices that they observed. They simulated cubic boxes ($Lx=d=1$) and tall thin boxes $(Lx/d<1)$, but the dynamics in wider boxes $(Lx/d>1)$, such as those that are typically used to study convection, were not examined there.

The top panel of Fig. 4 shows the volume-averaged kinetic energy as a function of time in a simulation of the elliptical instability in a wide box with $\epsilon =0.1,\u2009Ra=0$. The vortex observed previously in a 1 × 1 × 1 box by Barker and Lithwick^{5} dominates the flow to an even greater extent in the 4 × 4 × 1 box, as we can see by the dominance of the energy in the 2D component of the flow ($K2D$) at all times after the initial saturation. The 2D, or geostrophic, modes have energies ($K2D$) much larger than the inertial waves (quantified by the energy in the 3D modes, $K3D$), but the inertial waves undergo transient bursts temporarily increasing their energy, though $K2D$ remains dominant unlike in the 1 × 1 × 1 case in Barker and Lithwick.^{5} Each burst in the 3D energy later results in an increase in the 2D energy, indicating energy transfers from inertial waves to vortices. The vortex slowly decays viscously, however, the bursts of inertial wave energy are sufficient to compensate this lost energy, enhancing it further until a quasi-steady state is reached after $t\u223c0.1$. The corresponding energy injection ($I3D$) in the bottom panel of Fig. 5 in black shows that the 3D energy increase is a result of a direct energy injection into those modes. The 2D injection ($I2D$, not shown) is oscillatory in sign and has a small value consistent with 0. Meanwhile, the bottom panel of Fig. 4 shows the clear dominance of the LSV in a purely convective simulation with *ε* = 0, $Ra=4Rac$, as the 2D energy of the LSV, and by extension the LSV itself, continuously grows for all times plotted. The LSV will continue to grow until it reaches either the horizontal box-scale or its growth is balanced by viscous dissipation. A steady level of 3D energy is present in this simulation after the initial saturation, representing the energy in the convective eddies.

The interaction of convection and the elliptical instability varies according to the parameters chosen. First, we present a simulation with weak convection but strong ellipticity in Fig. 5, with $Ra=4Rac$ and $\epsilon =0.1$. The convection in this simulation leads to an LSV, which results in continuous growth of the 2D modes. This, however, does not inhibit the elliptical instability, and a multitude of bursts is observed. The elliptical instability in fact enhances the energy in the 2D modes by at least one order of magnitude compared to the purely convective simulation in the bottom panel of Fig. 4, as the bursts input more energy into the LSV. There is a continuous decrease in the 2D energy, from *t *=* *0.1 to *t *=* *0.17. In this period of time neither the bursts, weakened by the strong vortex, nor the convective eddies provide enough energy to compensate the viscous dissipation. The corresponding energy injection in the bottom panel of Fig. 5 in green roughly matches the energy injection of the purely elliptical simulation, although it is initially maintained at a higher value. From *t *=* *0.11 onwards, when the LSV has reached its strongest value, the behavior gradually changes from bursts to an almost continuous energy injection. From *t *=* *0.21, once the LSV has been sufficiently weakened, bursty behavior is again observed.

### B. Varying strength of convective driving and ellipticity

In Fig. 6, we present results for a range of values of $Ra/Rac$ and *ε* with $Ek=5\xd710\u22125.5$. The figures on the left show the time evolution of the kinetic energy components, and those on the right the energy injection term $I3D$. Figures 6(a) and 6(b) show simulations with $Ra=6Rac$ and $\epsilon =0.2$. Barker and Lithwick^{5} observed a change in behavior at $\epsilon \u22730.15$, seeing a sharp increase in the frequency and strength of the elliptical instability bursts. The 3D component of the kinetic energy is maintained at a higher level in this case, but is still lower than the energy in the 2D component. The energy injection features many bursts in a short time frame, with multiple bursts injecting energy at the same rate as the initial burst in linear growth phase. The increased burst frequency also leads to a sustained energy injection throughout the simulation. There appears to be a secondary transition around *t *=* *0.045 where the energy injection increases steeply and maintains a significant non-zero energy injection, much larger than the initial burst. We observe a correspondingly higher minimum level of the 3D component of the energy during this simulation.

The kinetic energy in Fig. 6(c) shows that increasing the Rayleigh number, i.e., making the convection stronger compared to the elliptical instability, results in fewer visible bursts into the 3D component in the first half of the simulation compared with the top panel of Fig. 5, and the total kinetic energy is further dominated by the 2D component. The increased convection strength, and therefore (for our parameters) stronger LSV compared with Fig. 5, drowns out most of the bursts from the elliptical instability. The reduced presence of the elliptical instability is also clearly visible from the $I3D$ term in Fig. 6(d), showing considerably fewer bursts in the first half of the simulation, decreasing the tidal dissipation. On the other hand, a “floor value” corresponding to a non-zero continuous energy injection arises, which is most clearly visible in between bursts. This sustained energy injection arises from the interaction between the convection and the equilibrium tidal flow.

Figure 6(e) and especially Fig. 6(f) confirm this sustained injection occurs as the convection is strengthened relative to the elliptical instability. The ellipticity has been reduced to $\epsilon =0.05$, thus weakening the elliptical instability (whose growth rate is proportional to $\epsilon )$. As a result, the bursts from the elliptical instability have vanished, with only a short initial burst remaining, after which a continuous injection arises. These simulations suggest there is a point at which the convection, with both its LSVs and its resulting effective viscosity acting to damp the inertial modes, overpowers the elliptical instability such that the bursts are completely suppressed.

Increasing the Rayleigh number to $Ra=20Rac$ and maintaining $\epsilon =0.1$ instead of decreasing *ε* in Figs. 6(g) and 6(h) leads to similar behavior, with no bursty behavior for the elliptical instability and instead a sustained energy injection. Thus, increasing Ra inhibits bursts even if they were present at lower values of the Rayleigh number. Additionally, the sustained injection term has increased by a factor of about 20 compared to Fig. 6(f). The sustained energy injection then increases with Ra and *ε*. Thus, introducing convection has two effects: (1) the bursts of elliptical instability are suppressed to a greater extent as the strength of convection increases relative to the strength of the elliptical instability, and (2) a sustained energy injection arises from the interaction between convection and the background flow.

### C. Heat transport modification by elliptical instability

A further effect of the elliptical instability is to modify heat transport. The inertial waves excited by the elliptical instability are capable of transporting heat in the system.^{15,28} The elliptical instability can occur and affect heat transport in both stably stratified and convectively unstable fluids [see Eq. (7)]. We observe that the heat transport is tied to the bursts of elliptical instability, and is similarly bursty in its operation. The Nusselt number represents the heat transport in the simulation and is shown as a function of time for $\epsilon =0.1,\u2009Ra=4Rac$ in the top panel of Fig. 7. This shows the two-sided effects of the elliptical instability in this simulation. The bursts increase heat transport, temporarily increasing the Nusselt number. Then, the enhanced cyclonic vortex as a result of the elliptical instability slightly decreases the heat transport, similar to the reduction observed in the presence of convective LSVs by Refs. 35 and 36, which they proposed occurs because cyclonic vortices act to effectively increase the rotation, thus further constraining the vertical motions, and as a consequence the heat transport, according to the Taylor–Proudman theorem.^{43}

In stably stratified fluids ($Ra<0$), or in convectively stable but unstably stratified fluids ($Ra<Rac$), in the absence of the elliptical instability, there are no sustained vortices or vertical motions. The heat transport in this regime would then be purely conductive (at long times, after decay of transients) with $Nu=1$. The vertical motions introduced by the elliptical instability in either regime may though transport heat during the bursts, even for convectively stable fluids.

The effects of the elliptical instability on heat transport have been quantified using time averages of the Nusselt number in the bottom panel of Fig. 7, for simulations performed with *ε* = 0 and $\epsilon =0.1$. Due to the absence of bursts of the elliptical instability, there is no observable enhancement to the heat transport when the convection is strong. In fact, since the ellipticity of the equilibrium flow slightly enhances the 2D energy of the vortex it is likely to result in a slightly diminished total heat transport, possibly supported by the reduced Nusselt number in the inset at $Ra=15Rac$ and $Ra=20Rac$. At low positive Rayleigh numbers, the elliptical instability plays a major role in enhancing or hindering the heat transport as a function of time. For very weak convection, in which there is no convectively generated LSV, the elliptical instability enhances the net heat transport strongly, while at higher convection strengths it cancels or slightly decreases the heat transport. Finally, the elliptical instability does indeed produce heat transport in stably stratified regimes, although the additional heat transport is highly variable in time, only occurring during a burst in energy injection of the elliptical instability, and decreases as the stratification increases (i.e., $\u2212Ra$ increases). This is represented by the Nusselt number tending toward one on average as the fluid becomes more stably stratified, where we note that all cases plotted are linearly unstable.

### D. Growth rate of 2D convective vortices

In Fig. 8, we examine the growth rate of the 2D energy $K2D$ produced by convection during the initial burst as a function of the Rayleigh number [i.e., we determine $\sigma 2D$ defined by $K2D\u221d\u2009exp\u2009(\sigma 2Dt)$ in the initial phases of the simulation]. We examine here simulations with *ε* = 0 and ones with higher *ε*, but all in the sustained energy injection regime. We find that non-zero ellipticities have little effect on this growth rate, modifying it slightly but not significantly. The growth rate $\sigma 2D$ is observed to increase as the Rayleigh number increases. Deviations in this growth behavior are observed at higher Rayleigh numbers ($>10Rac$). A potential explanation could lie in the slow increase in the growth rate of the 2D component as the 3D component of the kinetic energy grows. Because the growth of the 3D component is rapid at these values of the Rayleigh number, the non-linear breakdown happens before the 2D component reaches its maximum growth rate. Comparing the scaled 2D growth rate with the convective velocity in each simulation (plotted as the blue diamonds), and by extension the convective Rossby number of our simulations [$Roconv=\u27e8uz2\u27e9/(2\Omega d$)], we find good agreement. This implies that the growth of the 2D modes scales with the Rossby number in the initial regime, i.e., $\sigma 2D\u221dRoconv$. Note that this is not consistent with the predicted scaling due to weak nonlinear interactions between pairs of inertial waves (in the absence of convection) of, e.g., Kerswell,^{44} which would predict $\sigma 2D\u221dRo2$. It however agrees with the prediction for generation of an LSV from interactions between inertial waves and geostrophic modes at sufficiently large Rossby number,^{45} though a theory for convective LSVs has not yet been presented.

## IV. ANALYSIS OF THE SUSTAINED ENERGY INJECTION

### A. Origin and parameter regime for sustained tidal energy injection

We have performed a range of simulations in which *ε* and Ra are varied for $Ek=5\xd710\u22125.5$ to determine when sustained energy injection (vs burstiness) is obtained. First, a “phase diagram” was created based on these simulations, plotted in Fig. 9, which indicates in which simulations we observe bursts of the elliptical instability, in which we observe sustained energy injection. Simulations containing any bursts of the elliptical instability, even just at the earliest times have been labeled bursty (orange markers), while simulations containing no such bursts have been labeled sustained (blue markers). Bursty simulations may still feature a sustained energy injection, however in the interest of determining where the bursts of the elliptical instability still exist we have classified them as bursty here. Some of the simulations are very difficult to tell by eye whether they are bursty or sustained, and have therefore been labeled as uncertain (yellow markers). The transition between the two “phases” is likely situated close to these markers. Finally, one point has been labeled “neither,” namely, the point corresponding to the simulation with $Ra=2Rac,\u2009\epsilon =0.02$. This simulation features no bursts as the elliptical instability is too weak to operate at this level of convection. However, it also does not display any sustained energy injection. Of further note is that at this supercriticality ($Ra/Rac$), the LSV is absent, which was also observed by Guervilly, Hughes, and Jones^{35} to be independent of Ekman number.

To explain the absence of bursts of elliptical instability, we investigate the 3D modes in the flow in both real space and Fourier space. In real space, we use a method akin to the one used in Favier, Guervilly, and Knobloch^{46} to determine the 3D flow. The 3D velocity components are obtained by taking the difference between the total velocity and the *z*-averaged horizontal (or 2D) velocity components

Here, $\u27e8ux\u27e9z,\u2009\u27e8uy\u27e9z$ are the depth averaged *x* and *y* components of the velocity, respectively, i.e., the horizontal velocity components. From this, the magnitude of the 3D velocity is calculated. Favier, Guervilly, and Knobloch^{46} showed that the convective LSV suppresses 3D motions, resulting in an area with lower 3D velocities inside the vortex. Since we observe similar LSVs, such a suppression of 3D modes might contribute to the suppression of the elliptical instability.

First, we examine a case of the elliptical instability in isolation using this method in Fig. 10(a), with parameters $\epsilon =0.1,\u2009Ra=0$. We show results for the total velocity magnitude $u3D$ at *t *=* *0.05 at the mid-plane ($z=d/2$) during a burst of instability, after vortices have formed following initial saturation. In Fig. 10(a), we see that the power during a burst is concentrated inside the vortex, particularly in its center.

On the right panel, in Fig. 10(b), we show a similar result from a simulation with the elliptical instability and convection, using the parameters $\epsilon =0.1,\u2009Ra=6Rac$ during a burst of the elliptical instability at *t *=* *0.1. We again observe that the burst is concentrated in the center of the vortex. Examining the same simulation in the absence of a burst of the elliptical instability as well as a simulation with $Ra=6Rac,\u2009\epsilon =0.05$ (not pictured) reveals that there is no such concentration of power in the center. Instead the same picture of suppression of convective 3D motions is obtained as found by Favier, Guervilly, and Knobloch.^{46}

Bursts of the elliptical instability are thus primarily concentrated in the center of LSVs according to these results. The LSV is, therefore, expected to have a strong effect on the growth of the elliptical instability, as the free inertial waves existing within these vortices will differ from those of the original flow, thus acting as a constraint to detune the elliptical instability.

### B. Frequency-wavenumber spectrum analysis

We now present further analysis of our simulations using the approach devised by Le Reun *et al.*^{9} by computing the frequency-wavenumber spectrum of the flow to identify the inertial waves and the convection. This Fourier space analysis shown in Fig. 11 uses two properties of the elliptical instability to identify it in the spectrum: (1) it has a preferred direction (wavevector orientation) and (2) the dispersion relation of the inertial waves relate each direction to a particular frequency.

The direction of the flow, i.e., angle the wavevector makes with the rotation axis *θ*, is obtained from the ratio $kz/k=cos\u2009(\theta )$. Each velocity component, i.e., $ux,uy,uz$, associated with a specific wavevector can be put into a bin corresponding to its angle and wavenumber (wavevector magnitude) at every time step in a simulation. We use 60 bins for the angle *θ*, with $\theta =[0,\pi /2]$; likewise we have chosen bins of size $\pi /2$ for the *k*-bins and enough of these to cover all values allowed by the spatial resolution of the simulation. To obtain a spectrum as a function of the frequency and angle we sum over the *k*-bins, resulting in the total velocity component in each wavevector angle bin, at each time step. Equal timesteps of size $10\u22126$ are used, such that the fastest inertial waves, with periods $\pi \xd710\u22124.5$ can be properly captured. The Fourier transform in time of all three velocity components gives the frequency spectrum of each velocity component. We then multiply the transformed velocities with their complex conjugates and add all three components to obtain the energy in each *ω* and *θ* bin. We consider the interval of wavenumber bins $k=[2,50]$ to avoid the contribution of small-scale motions, which contain little energy. The geostrophic modes with $\theta =\pi /2$ strongly dominate the energy, so for clarity, we set the rightmost column at $\theta =\pi /2$ to zero on these plots since we wish to analyze the waves.

We plot the dispersion relation of inertial waves (in the absence of vortices and stratification) as a solid black line on all $\theta \u2212\omega $ energy spectra. Thus, the dispersion relation given in Eq. (4) is plotted. This choice is suitable for the Rayleigh numbers plotted, because convection tends to reduce the magnitude of the buoyancy frequency in the bulk of the box, leading to an effective buoyancy frequency $Neff2>N2$ (keeping in mind that *N*^{2} is negative). At the plotted Rayleigh numbers of $Ra=4Rac,\u2009Ra=6Rac$, and $Ra=8Rac$ the respective effective buoyancy frequencies are: $Neff2\u2248\u22121.5Rac,\u2009Neff2\u2248\u22122.5Rac,\u2009Neff2\u2248\u22123Rac$. Therefore, implying $Neff2/\Omega 2\u223cO(10\u22122)$ at the Rayleigh numbers used in these simulations; thus, the second term in Eq. (6) is close to zero, only affecting the dispersion relation around $\theta \u2248\pi /2$. The initial bursts of the elliptical instability are expected to be located at their preferred angle of $\theta =arccos(1/2)=\pi /3$ (when *n *=* *0). Combined with their dispersion relation, the fastest growing mode of the elliptical instability is, thus, expected to be located at $\theta =\pi /3,\u2009\omega =\Omega $ on these figures. We can, thus, very easily identify the elliptical instability in such a Fourier spectrum. An example where we can clearly identify the elliptical instability is given in Fig. 11(a), computed from the linear growth phase of the simulation with $Ra=0,\u2009\epsilon =0.05$. All modes with non-negligible energy during the linear growth phase are shown to be centered on this point, as well as at $\theta =\pi /3,\u2009\omega =3\Omega $, where the latter result from “nonlinear” interactions between the background tidal flow with frequency $2\Omega $ (and wavenumber zero) and the dominant modes at $\omega =\Omega ,\theta =\pi /3$.

Beyond the initial linear growth phase, inertial wave breakdown is observed, resulting in power concentrated around the initial fastest growing mode, but with a distribution primarily following the inertial wave dispersion relation. There is also energy in the geostrophic modes ($\theta \u223c\pi /2$, not shown) as well as the “mirrored dispersion relation” from secondary non-resonant interactions of the waves with the tidal flow.^{9} Figure 11(b) shows the same simulation as Fig. 11(a) but after the linear growth phase. Most of the power is concentrated around the initial fastest growing modes, however the energy is also spread throughout the figure, away from the dispersion relation, i.e., the resulting energy is no longer solely contained within the set of inertial waves.

Figure 11(c) shows the spectra for a simulation of rotating convection without the elliptical instability, with $Ra=6Rac$, *ε* = 0 from *t *=* *0.11 to *t *=* *0.13. Convection is shown to introduce modes at high values of *θ*. This can be understood from linear growth rate predictions, where we can show that convective instability of the conduction state requires $\theta \u2248[1.4,\pi /2]$ for *n *=* *1 modes at this Rayleigh number. The dominant modes are indeed concentrated in convective modes at $\theta \u2248[1.4,\pi /2]$ in this figure. The power away from these modes broadly follows the dispersion relation. This is particularly interesting because convective modes at onset are steady, so they should have $\omega \u223c0$ and be concentrated at the bottom of this figure. We might, therefore, speculate that the turbulence generated by convection is swept up by rotation into inertial waves, explaining the frequency of these modes. Inertial waves in rotating convection are expected to arise from oscillatory convective modes if $Pr<1$ (technically for Pr < 0.677)^{30} and have previously been observed in simulations^{47} at $Pr<1$. Supporting our argument for inertial waves arising due to rotating convection at $Pr\u22651$ is the detection of inertial waves at $Pr>1$ in spherical shell simulations of convection.^{48}

Spectra featuring both the elliptical instability and convection would be expected to look like a combination of these features, although it is likely to be difficult to distinguish the turbulent phases of the elliptical instability from the convective motions. However, the location of the elliptical instability bursts should shed some light on whether the sustained energy injection contains a weak (overshadowed) burst or whether the elliptical instability is absent entirely. The simulation with $Ra=6Rac,\u2009\epsilon =0.05$ analyzed from *t *=* *0.11 to *t *=* *0.13 is shown in Fig. 11(d). This shows the expected convective modes, but no enhanced power at the expected location of the elliptical instability. Thus we conclude that the operation of the elliptical instability has been inhibited by convection at these parameters.

The convective motions are expected to be small-scale motions, based on the visible fluctuations in Fig. 10(b), three and on the linear theoretical analysis of convection, which at these parameters predicts unstable modes in the range of wavenumber bins $k=[30,50]$. Meanwhile, the energetically dominant inertial waves are likely to be large-scale. This is a direct consequence of the elliptical instability being directional and, thus, choosing the mode with the smallest viscous effects (the largest possible modes, which also have the longest nonlinear cascade times). Therefore, we have reproduced the plots in Fig. 11 with the limited wavenumber range of $k=[2,12]$ in Fig. 22 in Appendix C.

Using the Fourier spectrum, we can determine the wavenumbers that contain the most energy in the simulation as a function of *θ* at a given *ω*. To this end, we do not sum over all wavenumbers *k*, but instead construct a $\theta \u2212k$ spectrum in Fig. 12. We are interested in the wavenumber magnitudes that are active on the dispersion relation. Therefore, we have done a Fourier transform on the wavenumbers, resulting in a $\omega \u2212k\u2212\theta $ matrix. Slices were then taken of the dispersion relation by taking the energy in all *k*-values at a combination of *ω* and *θ* that lies on the dispersion relation. The left panel of Fig. 12 shows the same simulation as Fig. 11(a), i.e., pure elliptical instability during its linear growth. The power here is concentrated along various black solid curves defined by the relation between *θ*, the vertical wavenumber, $kz=nz\pi $, and total wavenumber, *k*: $\theta =arccos(kz/k)$. Each curve has a different integer vertical wavenumber *n _{z}*, with

*n*= 1 the lowest curve,

_{z}*n*= 2 the one above, etc. The curves with $nz=1,2,3,4,5,6$ are plotted. In this simulation the energy of the elliptical instability burst is concentrated in modes with

_{z}*n*= 2, $k=4\pi $. The mode corresponding to this with horizontal wavenumber integers

_{z}*n*=

_{x}*n*is the (5, 5, 2) mode. This method of analysis is powerful as it clearly shows which inertial modes are growing in the simulation. In the right panel of Fig. 12, this method is applied to the same simulation as in Fig. 11(d) of both convection and the elliptical instability in the sustained regime. This panel shows that there is less power on the (5, 5, 2) mode. Instead, we see power concentrated on the

_{y}*n*= 1 curve, and concentrated toward higher wavenumbers and higher angles, as we would expect of convective motions. Indeed, we again observe no clear sign of the elliptical instability during this simulation.

_{z}Finally, we are interested in analyzing further the energy injection term $I3D$. The 3D component is any component, that is, not in the rightmost column of these $\theta \u2212\omega $ spectra (since that column has $\theta =\pi $, implying that *n _{z}* = 0). Therefore, to study $I3D$ we just set the rightmost column to zero, which was already done for visibility in the above plots. $I3D$ is calculated from the Reynolds stress components $uxuy,\u2009ux2$ and $uy2$. After calculation, this quantity is Fourier transformed, and its real part is then plotted here. The colorbar minimum has been increased compared to previous figures to reduce the impact of convective noise on the figure. On the left of Fig. 13, we show just the elliptical instability, for the same simulation as Fig. 11(a). We see that the energy injection is into the resonant inertial waves during this initial burst. It seems two frequency bins, in particular, contain the majority of the energy injection, one on the dispersion relation and one on the mirrored dispersion relation. To compare, $I3D$ of the same simulation as Fig. 11(d) is plotted on the right panel of Fig. 13. Energy injection is present predominantly on the right-hand side of the figure and is no longer concentrated along the dispersion relation, suggesting the energy injection is primarily into convective motions, instead of the inertial waves of the elliptical instability. This would be consistent with treating the energy transfer between tidal and convective flows as being due to a turbulent effective viscosity from the convective motions.

### C. Linear growth of elliptical instability on a convective background with an LSV

Based on the real space analysis in Sec. IV A, we concluded that the convective flow and its resulting LSV modifies subsequent growth of the elliptical instability, similar to the modification of the LSV resulting from the elliptical instability. Furthermore, based on the Fourier space results in Sec. IV B, we found that the elliptical instability is largely inhibited by the convective flow (and its LSV) for lower values of *ε*.

To further examine the effects of convection on the elliptical instability, we analyze the growth rate of bursts of the elliptical instability on a convective background. First, we measured the growth rates in simulations without convection ($Ra=0$) as a reference. We split the results of each simulation into the initial burst and any further bursts. The initial burst should then be close to the linear prediction, while any subsequent bursts are affected by non-linear effects, such as the LSV. The growth rate of the bursts of the elliptical instability is obtained from (half) the growth rate of $K3D$, and our results normalized by $\sigma \epsilon $ are shown in Fig. 14. The growth rates of the initial burst without convection (blue diamonds) lie close to the linear theoretical prediction, $(9/16)\epsilon \gamma $, plotted as the solid black line. Further bursts of these simulations (orange diamonds), however, substantially deviate from this prediction. A large reduction of the growth rate is found, likely due primarily to the LSV created by the initial burst. Note that these growth rates reduce further as the simulation continues and the energy in the LSV grows.

To compare these with similar results in the presence of convection, we ran new simulations which have been initialized with the flow and temperature fields from a purely convective simulation long after saturation of initial instability. We started several simulations with various ellipticities $\epsilon >0.05$ from our simulation with $Ra=4Rac$ (*ε* = 0), and several with $\epsilon >0.1$ from our simulation with $Ra=8Rac$ (*ε* = 0). These results are shown in Fig. 14 using yellow circles for the initial burst and purple circles for further bursts at $Ra=4Rac$ and green and cyan squares at $Ra=8Rac$, respectively. Focusing first on $Ra=4Rac$, we see that the suppression of the initial burst of the elliptical instability occurs for $\epsilon \u22720.07$. Our previous results and phase diagram (Fig. 9) indicated bursts of instability for $\epsilon \u22730.05$. These results also show that the growth rate is strongly affected by the convection, as the markers are substantially below the theoretical growth rate. The further bursts show a wide spread of growth rates. The highest measured growth rates overlap with those of the simulations of the pure elliptical instability. This implies that the convective LSV inhibits the elliptical instability in the same way as the LSV generated by the elliptical instability itself, but can inhibit it more strongly.

Our results for the case of more turbulent convection with $Ra=8Rac$ show higher growth rates than at $Ra=4Rac$, particularly at $\epsilon >0.14$. This is possibly indicative of a reduced suppression of the elliptical instability or an enhancement of the growth rate as the convection becomes stronger. A possible explanation for the enhanced growth rate could lie in Eq. (7). Increasing the Rayleigh number increases $\u2212N2$ for the conduction state. However, at $Ra=8Rac$ the growth rate is only increased by a factor of $\u22481.13$ compared to $Ra=0$, and even at $Ra=20Rac$, the increase is only a factor of $\u22481.3$ compared to $Ra=0$. Furthermore, this factor is likely to be less important than this would predict, as convection acts to reduce $\u2212N2$, and the efficiency of rotating convection increases with^{43} Ra.

We compare these results with some theoretical arguments in the figure. We use the energy injection $I3D$ of each simulation as a function of the Rayleigh number to obtain an effective viscosity *ν _{eff}* that corresponds with this sustained rate of injection (strictly acting on the tidal flow rather than the waves). We add this effective viscosity to the fluid viscosity

*ν*to obtain the “total viscosity,” which is used to compute a viscous damping rate $\u2212(\nu +\nu eff)k2$. Predictions for the growth rate after introduction of an effective viscous damping rate corresponding with $Ra=4Rac$ and $Ra=8Rac$ are plotted in purple and cyan lines, respectively, assuming the dominant wavenumbers and resonance conditions are unchanged. Incorporating the microscopic viscosity and/or effective viscosity decay rates are both inconsistent with the numerically-obtained growth rates. Indeed, results from simulations with just the elliptical instability imply that the suppression by an LSV is much stronger than would be predicted by such an effective viscosity. Furthermore, the slope of the growth rate as a function of

*ε*for simulations initialized on a convective background deviates from the 9/16 prediction in a similar manner with and without the convective background. To modify this 9/16 value there must be some change in the resonance conditions and the dominant wavenumbers, due to either detuning as previously remarked upon in this work, or the phases of the inertial waves, to explain the burst behavior.

We investigate the dominant wavenumber in each simulation, to examine if weakening of the elliptical instability occurs because the LSV changes this dominant wavenumber, using the approach described in Sec. IV B. The $\theta \u2212k$ spectrum for inertial modes following the dispersion relation is shown with $\epsilon =0.15$ for two $Ra$ values in Fig. 15. This shows that there is indeed a modification of the dominant wavenumber of the initial burst when initializing on a convective background. The $\theta \u2212k$ spectrum in the left panel shows the first elliptical instability burst in the simulation with $Ra=4Rac,\u2009\epsilon =0.15$. The energetically dominant wavenumber is no longer the (5, 5, 2) mode satisfying the ideal resonance condition without convection, and instead the power is concentrated at *n _{z}* = 1 with

*θ*close to but larger than the ideal value $\theta =\pi /3$. In the right panel of Fig. 15 we show the same for the first elliptical instability burst, but with $Ra=8Rac,\u2009\epsilon =0.15$. At both values of the Rayleigh number the subsequent inertial wave breakdown results in power in larger wavenumbers, including the (5, 5, 2) mode, which is viscously dissipated but otherwise maintained until the next burst. The availability of the (5, 5, 2) mode in the subsequent bursts results in higher growth rates compared to the initial burst, but still far below the ideal linear prediction. However, we observe power at the expected $kz/k=1/2$ and $\omega =\gamma $, therefore apparently arguing against the hypothesis that detuning the dominant resonance is responsible for most of the reduction in growth rate. This leaves the perturbed phase argument originally proposed by Barker and Lithwick

^{5}to potentially explain the observed change of the growth rate prefactor.

## V. SCALING LAWS OF THE ENERGY INJECTION

The energy injection rate ($I3D$) due to the elliptical instability on its own scales consistently with $\epsilon 3$ when the flow is sufficiently turbulent.^{5,6} However, the sustained energy injection in our simulations is not the result of the elliptical instability in isolation. We plot the energy injection $I3D$ as a function of *ε* at various values of Ra at fixed $Ek=5\xd710\u22125.5$ in the top panel of Fig. 16, which we divide into two regimes by a vertical dashed line located at $\epsilon =0.08$, in accordance with our discussion of the simulations initialized on a convective background.

To the left of this line the simulations show sustained energy injection without obvious bursts for $Ra\u22732Rac$ whereas to the right of the line the simulations show clear bursts. We fit both sides separately using the $Ra=6Rac$ data. The black line is the fit to the sustained energy injection using points on the left side, which scales like $\epsilon 2$. This is predicted if the convection acts like an effective viscosity in Eq. (28).

The bursts of elliptical instability on the right side of the figure contribute on top of this sustained energy injection. We fit using both the (naive) theoretically predicted $\epsilon 3$ scaling and one like $\epsilon 6$ as previously observed.^{5} Both are consistent with the data on the right-hand side and are inconsistent with data on the left. Furthermore, the fits are consistent with data from simulations at all values of the Rayleigh number, indicating that this scaling may be independent of Ra. In this regime, the elliptical instability is much more efficient than the effective viscosity of convection and would only be surpassed by the latter when $\epsilon \u22480.01$ if we extrapolate the former with an $\epsilon 3$ scaling.

## VI. DISCUSSION AND CONCLUSION

### A. Comparison with previous work

As mentioned previously in Sec. II B, in the linear study of elliptical instability in an unbounded strained vortex, convection enhances the growth rate of elliptical instability.^{2} In our nonlinear simulations, however, we observe the opposite, as the suppression of elliptical instability increases with the increasing Rayleigh number. Linear analysis of elliptical instability in a heated cylindrical annulus on the other hand finds that the growth rate can decrease with the increasing Rayleigh number.^{27} Expanding on this linear analysis are the experiments of the elliptical instability with convection.^{28} These experiments also measured the growth rate of the elliptical instability with convection and obtained the same result. In addition, they found that smaller Ekman number leads to a faster growth rate of the instability. This result is also supported by previous numerical simulations in ellipsoids.^{15} Although none of these experiments or simulations have the same geometry as our simulations or clearly feature an LSV, a suppression of elliptical instability due to convection is also observed. We conclude that this result is likely to be universal that elliptical instability is weakened by convection. We find similar heat transport as a result of the elliptical instability as Cébron, Maubert, and Le Bars^{15} Lavorel and Le Bars,^{28} in the stably stratified regime, as well as the enhanced heat transport in simulations with the Rayleigh number below or just above the critical Rayleigh number. However, we do not observe a constant Nusselt number but rather observe a weakening of this effect as the stratification increases. We also find decreased heat transport at $Ra\u226bRac$ compared to the purely convective simulations, likely due to the stronger vortices formed.

### B. Future work

In our current setup, the local box is appropriate to model the poles of a planet with the gravity and rotation axis both pointing in the *z*-direction. The location of the box, and thus, accounting for misaligned gravity and rotation could have important effects on the interactions of these instabilities. Rotating convection has been studied with misaligned gravity and rotation by Currie *et al.*^{43} They found that convective plumes (for rapid rotation) align with the rotation axis, and that strong zonal flows tend to predominantly form rather than LSVs at non-polar latitudes. These zonal shear flows have an important effect on heat transport, but could also modify the excitation and saturation of inertial waves due to the elliptical instability that should be explored in further work. Global simulations that are sufficiently turbulent and rapidly rotating to capture regimes similar to those we have explored would also be worthwhile, somewhat along the lines of the previous laminar simulations presented in, e.g., Cébron, Maubert, and Le Bars.^{15} Following Barker,^{8} one could also study the interaction of the bursty non-linear dynamics of elliptical instability with convection. One advantage of global simulations (in full ellipsoids) is that the linearly-excited inertial waves are no longer constrained by the (artificial) aspect ratio of the box.

It would also be of interest to further explore the parameter regime in our simulations, particularly by varying the Prandtl number. In particular, the low Prandtl number ($Pr<0.67$) rotating convection itself excites inertial waves.^{47} These convectively created inertial waves might also be unstable to elliptical instability and could, due to their constant generation by the oscillatory convection, result in another source of potentially continuous tidal dissipation. Finally, Hot Jupiters, like Jupiter itself, tend to have strong magnetic fields.^{49} Therefore, the inclusion of MHD is likely to be important and can have significant effects on tidal dissipation. In these simulations the LSVs of convection and the elliptical instability are likely to be suppressed, as magnetic fields inhibit the formation of large-scale structures. This should allow for a continuous operation of elliptical instability.^{6} It is, however, unclear what the influence of convection will be in this interaction.

### C. Conclusion

We have investigated the interactions of elliptical instability and rotating Rayleigh–Bénard convection in a Cartesian model using psuedo-spectral hydrodynamical numerical simulations involving horizontal shearing waves. First, we simulated elliptical instability without convection in wide boxes (with stress-free impenetrable boundaries in the vertical) for the first time and found the nonlinear evolution of the instability to produce geostrophic vortices that dominate the flow to an even greater extent than in cubical boxes. The introduction of convection leads to a suppression of elliptical instability that we argue is primarily due to the convectively generated LSV. It also gives rise to a sustained energy injection into the flow (i.e., transfer from the elliptical/tidal flow) that scales as $\epsilon 2$, which can be interpreted such that the convection operates as an effective viscosity (independent of *ε*) in damping the tidal flow.

The suppression of elliptical instability by convection was investigated in detail using numerous approaches. Measuring the 3D motions, which are weakened by the LSV, we showed that during a burst of elliptical instability, the power is concentrated in the center of the vortex. We also presented a detailed analysis of the frequency and wavenumber Fourier spectra of the energy in our simulations to clearly identify inertial modes and convective flows. We observed that the elliptical instability (and energy injection into inertial modes more generally) is indeed inhibited by convective flows. Rotating convection also weakly excites inertial modes, which are identified as power in modes along the dispersion relation, in the absence of elliptical instability.

When initializing simulations of elliptical instability from a convective turbulent state including an LSV, it was found that this LSV reduces the growth rate of elliptical instability compared with the inviscid or viscous growth rate prediction. It is also reduced compared with the prediction modified by crudely adopting the aforementioned effective turbulent viscosity. The reduction of the growth rate by the LSV indicates that the dominant resonances are de-tuned by it or that there are significant perturbations in the phases of the waves by the LSV. Our Fourier space analysis showed that the fastest growing mode with an LSV is the same as the one found in the absence of an LSV for all bursts of elliptical instability. This indicates that the latter argument may be more applicable.

We also found that the inertial waves excited by elliptical instability can transport heat; when the elliptical instability is weak relative to convection or suppressed this has little effect on the Nusselt number, but when the elliptical instability is comparable in strength to the convection, it can significantly enhance transport. The elliptical instability can also result in heat transport in stably stratified regimes, but this weakens as the stratification becomes stronger (within the linearly unstable regime).

The elliptical instability leads to an energy transfer rate from the tidal/elliptical flow (and hence dissipation rate), that is, approximately proportional to $\epsilon 3$, as previously found in the absence of convection.^{5,8} This scaling is similar to results obtained for related instabilities such as the precessional instability.^{38,39} Indeed, when the elliptical instability operates, the energy transfer rates are quantitatively similar to those found in prior work.^{6} This implies that when it is not suppressed by convection, the astrophysical energy transfer rates (e.g., in hot Jupiters) from the elliptical instability are negligibly affected by convection. However, we should point out that in the narrow range of simulations where the $\epsilon 3$ scaling is obtained without being suppressed by convection, the data are also consistent with a stronger $\epsilon 6$ scaling observed previously. Further work exploring more turbulent regimes at higher Ra and smaller Ek and Pr would be beneficial to further explore the scaling laws to allow robust extrapolation to stars and planets.

## ACKNOWLEDGMENTS

We thank T. Le Reun for kindly supplying his routines in Snoopy that were modified and used to produce the Fourier spectra in Sec. IV. We also thank the three referees for their positive and constructive comments that helped us to improve the article. N.B.V. was supported by EPSRC studentship No. 2528559. A.J.B. and R.H. were supported by STFC Grant Nos. ST/S000275/1 and ST/W000873/1. R.H. would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme “Frontiers in dynamo theory: from the Earth to the stars” where work on this paper was undertaken. This work was supported by EPSRC Grant No. EP/R014604/1. R.H.'s visit to the Newton Institute was supported by a grant from the Heilbronn Institute. This work was undertaken on ARC4, part of the High Performance Computing facilities at the University of Leeds, and using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital Grant Nos. ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations Grant No. ST/R001014/1. DiRAC is part of the National e-Infrastructure.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Nils Beijing de Vries:** Conceptualization (equal); Formal analysis (lead); Software (equal); Writing – original draft (lead); Writing – review & editing (equal). **Adrian John Barker:** Conceptualization (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). **Rainer Hollerbach:** Conceptualization (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: RESOLUTION

Multiple tests were performed to ensure that our simulations were properly resolved. We ensure that heat transport is well resolved by testing the Nusselt number. If upon increasing the vertical (or horizontal) resolutions (*n _{x}*,

*n*,

_{y}*n*are the numbers of grid points in each direction), the Nusselt number was negligibly altered, i.e., the previous resolution was suitable for resolving the convection. In Fig. 17, the Nusselt number for different vertical resolutions for one of our most demanding simulations with $Ra=15Rac$ and $\epsilon =0.1$ ($Ek=5\xd710\u22125.5$) is plotted. It is clear that too small vertical resolution influences Nu and that we get convergence numerically in this case when $nz\u2265160$. Similar test simulations were done for all Ra to ensure a good vertical resolution.

_{z}A horizontal wavenumber power spectrum of the kinetic energy showing simulations with various horizontal box sizes is shown in Fig. 18 for a demanding case with $Ra=20Rac$ and $\epsilon =0.2$ ($Ek=5\xd710\u22125.5$). The LSV caused by convection leads to the smallest wavenumber modes becoming dominant. Choosing a larger box only serves to let the vortex grow larger over time. However, during the initial burst phase at the start of the simulation, around *t *=* *0.005, the LSV is still forming and power is not yet contained in the largest scales. Therefore, we can check these early phases to see whether the box size is appropriate to accommodate the convection and elliptical instability. Additionally, we can examine the power at the anti-aliasing scale, as this will reveal whether the flow is well resolved. We desire that the power here is at least a factor of 10^{3} lower than that in the peak to consider a simulation to be “well resolved.” Based on the horizontal power spectra in Fig. 18, we conclude that a horizontal resolution of 256 × 256 with a box size of 4 × 4 is suitable. In addition, two power laws are plotted, in dashed-black, the $k\u22a5\u22125/3$ associated with the Kolmogorov spectrum of turbulence, which matches the middle or inertial subrange of all spectra quite well, except for 8 × 8 × 1. Furthermore, the dashed-red $k\u22a5\u22123$ power law associated with wave turbulence is plotted, possibly matching the smaller wavenumbers.

Increasing Ra leads to more turbulent simulations, requiring higher resolutions to accurately capture small-scale effects, which becomes computationally expensive. To minimize computational expense, we desire to minimize resolution subject to the simulation being well resolved. One further check is that the most important quantity we study, the energy injection term *I*, is numerically converged. To this end, an additional test was performed with a simulation with $Ra=20Rac,\u2009\epsilon =0.1$ and a resolution of $512\xd7512\xd7224$. The energy injection term $I3D$ of this simulation is compared with one with resolution $256\xd7256\xd7224$ in Fig. 19. Aside from small fluctuations the two simulations are in agreement, indicating that horizontal resolutions of 256 × 256 are appropriate to study the energy injection accurately. The resolutions used for all Ra at fixed Ek are given in Table I.

$Ek=5\xd710\u22125.5$ . | $nx\xd7ny$ . | $nz$ . |
---|---|---|

$Ra/Rac=\u22126,\u22124,\u22123,\u22121,\u22120.8,0.3,0.8,1.99,3,4,6$ | 256 × 256 | 96 |

$Ra/Rac=7,8$ | 256 × 256 | 128 |

$Ra/Rac=\u221210,9,10,15$ | 256 × 256 | 160 |

$Ra/Rac=20$ | 256 × 256 | 224 |

$Ek=5\xd710\u22125.5$ . | $nx\xd7ny$ . | $nz$ . |
---|---|---|

$Ra/Rac=\u22126,\u22124,\u22123,\u22121,\u22120.8,0.3,0.8,1.99,3,4,6$ | 256 × 256 | 96 |

$Ra/Rac=7,8$ | 256 × 256 | 128 |

$Ra/Rac=\u221210,9,10,15$ | 256 × 256 | 160 |

$Ra/Rac=20$ | 256 × 256 | 224 |

### APPENDIX B: DIFFERENT BOX SIZES

To test the effect of different box sizes on energy transfers and Fourier spectra, we ran multiple simulations with $Ra=6Rac$ and different *ε* ($Ek=5\xd710\u22125.5$). We have plotted $I3D$ as a function of *ε* with different markers denoting different box sizes in Fig. 20. The blue markers represent results at 1 × 1 × 1, the green markers are 2 × 2 × 1, the yellow markers 3 × 3 × 1 and the burgundy markers 4 × 4 × 1. The mean energy injection in the sustained regime is shown to be independent of box size, and all markers follow the same $\epsilon 2$ scaling. The values on the right in the bursty regime do change, as the 4 × 4 × 1 results becomes bursty for $\epsilon \u22650.05$; whereas simulations with smaller horizontal box sizes do not. Results at 3 × 3 × 1 only indicate burstiness at $\epsilon \u22650.1$, while 2 × 2 × 1 and 1 × 1 × 1 remain in the sustained regime beyond $\epsilon =0.1$.

The explanation for this behavior lies in the allowed values of $k\u22a5$ as the box size is varied. For smaller box sizes, the values of $k\u22a5$ and, hence, *k* increase (for the same *k _{z}*), so for a resonance with $kz/k=\xb11/2$,

*k*(and hence

_{z}*k*) must also be larger. Figure 21 shows the $\theta \u2212k$ spectrum on the dispersion relation for 1 × 1 × 1, $Ek=5\xd710\u22125.5,\u2009Ra=0,\u2009\epsilon =0.1$, during the initial burst of elliptical instability; these are the linearly most unstable modes. The dominant

*k*modes lie on lines of

*n*=

*4 and*

*n*=

*5. These larger values of*

*k*imply larger decay rates $\u2212\nu k2$, and therefore, a decreased growth rate due to viscosity. This suppresses the elliptical instability for larger

*ε*when the box is smaller. The suppression of the elliptical instability is, thus, artificially enhanced (reduced) by the choice of a smaller (larger) box. Upon extrapolating this effect to a full planet, it is expected that the viscosity suppression of the elliptical instability is weak due to the large scales available to the system.

### APPENDIX C: Limited wavenumber range frequency spectra

We present the same plots as in Fig. 11 with a limited wavenumber range of $k=[2,12]$. One effect of the limited wavenumber range, combined with our finite grid, is that a number of columns on the right will contain no energy. Only higher wavenumbers can have these angles. The limited *k* range does not affect the linear growth spectrum in Fig. 22(a) as the power is concentrated in wavenumbers within our adopted range. During the inertial wave breakdown in Fig. 22(b), we can clearly see the power concentrated along the inertial wave dispersion relation, as well as the mirrored dispersion relation representing the secondary non-resonant interactions between the waves and the background tidal flow.^{9} In the convective simulations in the bottom two panels, there is indeed power along the dispersion relation where inertial waves are expected, providing another tentative hint for inertial waves in rotating convection.