© ESO, 2014

1. Introduction

The study of tidal friction has a long history in geo- and astrophysics (e.g., Darwin 1879). The basic picture of a tidally deformed fluid body is that of an ellipsoid with the long axis pointing towards the tidal companion, and in which the fluid rotates with respect to the major axes of the ellipsoid, so that each fluid particle follows an elliptical streamline. It was discovered in the mid 1970s that flows of constant vorticity with elliptical streamlines are prone to a hydrodynamic instability, the so-called elliptical instability (for a review, see Kerswell 2002). One expects to find vastly different dissipation rates in laminar and in unstable turbulent flows. A determination of the stability limit of the basic tidal flow is thus an indispensable prerequisite to any computation of tidal friction. In this paper, we compute the stability limit within a simple model of tidal flow taking into account the compressibility of a stellar or planetary envelope.

Fluids in solid body rotation (i.e., flows with circular streamlines) support a type of waves known as inertial waves in which the Coriolis force acts as a restoring force. In elliptical instability, the finite ellipticity couples inertial modes in pairs and can lead to instability if certain resonance conditions are met. The same phenomenon is often called a triad resonance, the triad being formed by the two inertial waves and the basic, elliptically deformed flow which is itself counted as third wave.

In this paper, we calculate the growth rate of the elliptical instability in a slightly deformed sphere. Incompressible flows have already been investigated in this geometry (Lacaze et al. 2004; Le Bars et al. 2010). Elliptical instability can be studied in its simplest form in an unbounded domain (Bayly 1986; Landman & Saffman 1987; Waleffe 1990; Miyazaki & Fukumoto 1992; Miyazaki 1993; Miyazaki et al. 1995). Another convenient geometry is an elliptically deformed cylinder (Gledzer et al. 1975; Malkus 1989; Eloy et al. 2003), and other studies exist for deformed spheroids (Gledzer & Ponomarev 1977; Kerswell 1994; Cébron et al. 2010a) and spherical shells (Aldridge et al. 1997; Seyed-Mahmoud et al. 2000; Cébron et al. 2012b).

Several additional ingredients have been added to the problem, such as stratification (Miyazaki & Fukumoto 1992; Miyazaki 1993; Guimbard et al. 2010), magnetic fields (Kerswell 1994; Lacaze et al. 2006; Herreman et al. 2009; Cébron et al. 2012a), rotation of the elliptical perturbation (Craik 1989; Gledzer & Ponomarev 1992; Miyazaki 1993; Miyazaki et al. 1995; Seyed-Mahmoud et al. 2000; Le Bars et al. 2010), and viscosity (Landman & Saffman 1987; Kerswell 1994; Lacaze et al. 2004; Le Bars et al. 2010). The importance of the elliptical instability for tidal dissipation is examined in a recent paper by Barker & Lithwick (2013). They studied the non-linear evolution of the elliptical instability. They performed three-dimensional hydrodynamical simulations of a box with periodic boundary conditions with a base flow such that this box can be considered as a small patch of a tidally deformed fluid in a planet or a star. They found that for the astrophysically relevant values of the ellipticity the wave driving mechanism is not sustained permanently because of the presence of strong columnar vortices whose presence effectively suppresses the driving mechanism. This leads to less dissipation compared to the case with a sustained wave driving mechanism, but it is certainly higher compared to case with a stable flow. A strong stratification can prevent the process of re-laminarization as found in an experiment with a stratified fluid in a rotating cylinder by Guimbard et al. (2010). The same holds for the instability in the presence of a weak magnetic field (see Barker & Lithwick 2014, a companion paper to Barker & Lithwick 2013). Barker & Lithwick (2014) found that this field prevents the vortices from forming. They calculated the dissipation and conclude that the inferred tidal dissipation is potentially important at short orbital periods. In this companion paper they also neglect the effects of a realistic geometry and the additional presence of turbulent convection, both of which possibly enhance the dissipation such that the instability becomes important more generally.

The influence of a temperature-gradient on the elliptical instability in a triaxial ellipsoid was examined by Cébron et al. (2010b) with numerical methods. They found that the growth rate of the elliptical instability is significantly enhanced by a thermal stratification and that in a convective flow the elliptical instability can still grow, but with a reduced growth rate. They were not able to reach the regime of very large Reynolds numbers, where the results of Fabijonas & Holm (2003) predict an increased tidal destabilization. For stars the case of a bi-layer flow with a stable stratified and a convective region is interesting. In the paper by Cébron et al. (2010b) it was shown that even in such a flow the instability can grow over the whole fluid.

In this paper, we will include the rotation of the perturbation, viscous dissipation, and most importantly, compressibility. The last has been neglected in nearly all other publications on the subject with the exception of the work by Cébron et al. (2013) who numerically simulated one particular set of parameters. In contrast, we use a semi-analytical method to obtain a broad overview. The most important issue this study will address is how the compressibility of the flow influences growth rates of the elliptical instability. Surprisingly, we will find that the structure of the most unstable modes in ideal fluids depend on the compressibility, but the maximum growth rate is independent of compressibility. Section 2 collects all the formulas describing the model and Sect. 3 explains their numerical implementation. The results are presented in Sect. 4.

2. Mathematical formulation of the model

We consider the equatorial tide raised on a central body by a tidal perturber or perturbing body. We choose a frame of reference in which the perturbing body is at rest and use a Cartesian coordinate system x,y,z with its origin at the center of the central body, its z-axis directed along the rotation axis of the reference system, and the x-axis is pointing towards the perturbing body. This reference frame rotates at rate Ω_P = Ω_Pẑ (hats denote unit vectors) relative to inertial space. In this frame, the central body rotates about the z-axis with angular velocity $\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\frac{\mathit{b}}{\mathit{a}}\mathrm{+}\frac{\mathit{a}}{\mathit{b}}\mathrm{)}{\mathrm{\Omega }}_{\mathrm{F}}$ within an ellipsoid with semi-major axes a,b,c and the surface ${\left(\frac{\mathit{x}}{\mathit{a}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{y}}{\mathit{b}}\right)}^{\mathrm{2}}\mathrm{+}{\left(\frac{\mathit{z}}{\mathit{c}}\right)}^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathit{,}$(1)so that the motion of the central body within the chosen frame of reference is given by the velocity field u₀: ${{u}}_{\mathrm{0}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{F}}\left(\mathrm{-}\frac{\mathit{a}}{\mathit{b}}\mathit{y}{x̂}\mathrm{+}\frac{\mathit{b}}{\mathit{a}}\mathit{x}{ŷ}\right)\mathit{.}$(2)In order to determine the stability of this flow, we start from the full Euler equation

$\begin{array}{ccc}& & {\mathit{\partial }}_{\mathit{t}}{v}\mathrm{+}\mathrm{(}{v}\mathrm{·}\mathrm{\nabla }\mathrm{)}{v}\mathrm{+}\mathrm{2}{{\Omega }}_{\mathrm{P}}\mathrm{\times }{v}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{\varrho }}\mathrm{\nabla }\mathit{P}\mathrm{+}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{self}}\mathrm{+}\mathrm{\nabla }{\mathrm{\Phi }}_{\mathrm{P}}\\ & & {\mathit{\partial }}_{\mathit{t}}\mathit{\varrho }\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{\left(}\mathit{\varrho }{v}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$where v stands for the velocity, ϱ the density, P the pressure, and Φ_self for the potential terms created by the central body itself (gravitational and centrifugal), whereas Φ_P is the perturbing tidal potential. We will assume u₀ to be a stationary solution of the above equations for suitable density profiles and potentials ρ₀ and Φ₀: $\begin{array}{ccc}& & \\ & & \mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}}{{u}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{0.}\end{array}$

In order to end with a tractable problem, we will have to choose a density profile such that the eigenmode calculation below leads to a separable equation. This is achieved by setting ${\mathit{\rho }}_{\mathrm{0}}\mathrm{=}\begin{array}{c}\mathrm{˜}\\ {\mathit{\rho }}_{\mathrm{0}}\end{array}{\left(\mathrm{1}\mathrm{-}{\left(\frac{\mathit{x}}{\mathit{a}}\right)}^{\mathrm{2}}\mathrm{-}{\left(\frac{\mathit{y}}{\mathit{b}}\right)}^{\mathrm{2}}\mathrm{-}{\left(\frac{\mathit{z}}{\mathit{c}}\right)}^{\mathrm{2}}\right)}^{\mathit{\beta }}$(7)for arbitrary prefactors $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{\rho }}_{\mathrm{0}}\end{array}$ and exponents β (Wu 2005a).

The question arises which beta one should use for the calculations. We must first determine which polytropic index n in the polytropic relation $\mathit{P}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{K}{\mathit{\rho }}^{\mathrm{1}\hspace{0.17em}\mathrm{+}\hspace{0.17em}\mathrm{1}\mathit{/}\mathit{n}}\mathrm{\left(}\mathit{r}\mathrm{\right)}$(8)is appropriate for the central body. We choose n = 3, n = 3/2, and n = 0 for our calculations; n = 3 is suitable for stars which are well modeled by a relativistic completely degenerate electron gas, like relativistic white dwarfs. The same polytropic index also describes main sequence stars with M ≳ M_⊙, (Kippenhahn & Weigert 1990). n = 3/2 is appropriate for objects which are well modeled by a non relativistic completely degenerate electron gas, like non-relativistic white dwarfs and brown dwarfs. And n = 3/2 is also relevant for main sequence stars with a mass below M ~ 0.4  M_⊙, these are fully convective stars (Kippenhahn & Weigert 1990; Chabrier et al. 2009). For planets, the range of n is 0 ≤ n ≲ 1.5, depending on how massive they are (Horedt 2004). For Jupiter mass objects, n = 1 is a good value (Chabrier et al. 2009). ρ(r) in the polytropic case can be obtained by solving the Lane-Emden equation. This equation can be derived through the usage of (8), the equation for hydrostatic equilibrium and the Poisson equation (Kippenhahn & Weigert 1990). The values for ρ(r) in the polytropic case are taken from Horedt (1986). We determine the appropriate β by simply fitting the power law (7) to the ρ(r) for the polytropic profiles.

It is easily verified that ∇·(ρ₀u₀) = 0 and ∇ ×  { (u₀·∇)u₀ + 2Ω_P × u₀ }  = 0. The curl of the gradient terms in Eq. (5) is trivially zero, and $\mathrm{\nabla }\mathrm{\times }\mathrm{\left\{}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{\nabla }{\mathit{p}}_{\mathrm{0}}\mathrm{\right\}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{\nabla }{\mathit{\rho }}_{\mathrm{0}}\mathrm{\times }\mathrm{\nabla }{\mathit{p}}_{\mathrm{0}}$ which is zero in a polytropic atmosphere in hydrostatic equilibrium. u₀ and ρ₀ as given above are therefore solutions of the Euler equation for some perturbing potential, albeit not necessarily for the perturbing potential of a point mass at a finite distance. One may view ρ₀ either as an approximation to the density profile of the state excited by a tidal perturber idealized as a point mass, or as the exact profile for a perturbing potential which approximates a real tidal potential.

We will now consider the linear stability of the ground state. All quantities are decomposed into their value in the basic state indicated by an index zero and a perturbation which is considered to be small: ϱ = ρ₀ + ρ, v = u₀ + u, P = p₀ + p, and Φ_self = Φ₀ + Φ. The linearized equations are: $\begin{array}{ccc}& & {\mathit{\partial }}_{\mathit{t}}{u}\mathrm{+}\mathrm{(}{{u}}_{\mathrm{0}}\mathrm{·}\mathrm{\nabla }\mathrm{)}{u}\mathrm{+}\mathrm{(}{u}\mathrm{·}\mathrm{\nabla }\mathrm{)}{{u}}_{\mathrm{0}}\mathrm{+}\mathrm{2}{{\Omega }}_{\mathrm{P}}\mathrm{\times }{u}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{\nabla }\mathit{p}\mathrm{+}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{\nabla }{\mathit{p}}_{\mathrm{0}}\mathrm{+}\mathrm{\nabla \Phi } \\ & & {\mathit{\partial }}_{\mathit{t}}\mathit{\rho }\mathrm{+}\mathrm{\nabla }\mathrm{·}\mathrm{(}\mathit{\rho }{{u}}_{\mathrm{0}}\mathrm{+}{\mathit{\rho }}_{\mathrm{0}}{u}\mathrm{)}\mathrm{=}\mathrm{0.}\end{array}$We will restrict ourselves to well mixed ground states ρ₀,p₀ of constant entropy in atmospheres characterized by an adiabatic exponent γ, and only allow perturbations of the ground state which obey the adiabatic equation of state, so that $\mathrm{\nabla }{\mathit{\rho }}_{\mathrm{0}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}}\mathrm{\nabla }{\mathit{p}}_{\mathrm{0}}$ and p/p₀ = γρ/ρ₀ for small deviations ρ,p from the ground state. The right hand side of the linearized Euler equation then simplifies according to $\mathrm{\nabla }\mathrm{\left(}\frac{\mathit{p}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{\nabla }{\mathit{p}}_{\mathrm{0}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{\nabla }\mathit{p}$.

The equation of continuity is next simplified by invoking the anelastic approximation. This approximation is valid if the density deviations from ρ₀ are small and if velocities are small compared with the speed of sound. The latter condition is violated at the surface of a star or a planet where the speed of sound tends to zero. The validity of the anelastic approximation has already been discussed in detail in the more restricted scope of eigenmode calculations as they will be done below. It appears that the region where the anelastic approximation fails is too small to modify global results such as inertial mode frequencies (Ivanov & Papaloizou 2010), so that we shall adopt the anelastic approximation from here on.

The linear stability problem is now reduced to $\begin{array}{ccc}& & \\ & & \mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}}{u}\mathrm{\right)}\mathrm{=}\mathrm{0}\end{array}$with $\mathit{\psi }\mathrm{=}\frac{\mathit{p}}{{\mathit{\rho }}_{\mathrm{0}}}\mathrm{+}\mathrm{\Phi }$. We will solve these equations subject to the boundary condition on the ellipsoidal surface (1). These boundary conditions describe a solid wall in a laboratory experiment more directly than the surface of an atmosphere, but the results below suggest that the detailed choice of boundary conditions does not matter for the main findings of this paper.

We now follow a procedure similar to that used by Gledzer & Ponomarev (1992). We first remove dimensions by rescaling the Cartesian coordinates with their respective semi-major axes (which maps the ellipsoidal surface on a sphere of radius 1) and we rescale time with Ω_F, $\begin{array}{ccc}& & {\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{x}}{\mathit{a}}\mathit{,} {\mathit{y}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{y}}{\mathit{b}}\mathit{,} {\mathit{z}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{z}}{\mathit{c}}\mathit{,} {\mathit{u}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{u}}{{\mathrm{\Omega }}_{\mathrm{F}}\mathit{a}}\mathit{,} {\mathit{v}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{v}}{{\mathrm{\Omega }}_{\mathrm{F}}\mathit{b}}\mathit{,}\\ & & {\mathit{w}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{w}}{{\mathrm{\Omega }}_{\mathrm{F}}\mathit{c}}\mathit{,} {\mathit{t}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{F}}\mathit{t,} \mathit{R}\mathrm{=}\sqrt{\frac{{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{b}}^{\mathrm{2}}}{\mathrm{2}}}\mathit{,} {\mathit{\psi }}^{\mathrm{\prime }}\mathrm{=}\frac{\mathit{\psi }}{{\mathrm{\Omega }}_{\mathrm{F}}^{\mathrm{2}}{\mathit{R}}^{\mathrm{2}}}\mathit{,}\\ & & \end{array}$(13)where u, v, and w are the x, y, and z components of u, respectively. For simplification we consider only the case $\mathit{c}\mathrm{=}\mathit{R}\mathrm{=}\sqrt{\mathrm{(}{\mathit{a}}^{\mathrm{2}}\mathrm{+}{\mathit{b}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{2}}$. It will become plausible below that a change in c will not cause substantial changes. This restriction to certain ellipsoidal shapes will simplify the calculations below. In these new variables, Eq. (11) becomes after omitting primes $\begin{array}{ccc}& & \frac{\mathit{\partial u}}{\mathit{\partial t}}\mathrm{+}\mathit{x}\frac{\mathit{\partial u}}{\mathit{\partial y}}\mathrm{-}\mathit{y}\frac{\mathit{\partial u}}{\mathit{\partial x}}\mathrm{-}\mathit{v}\mathrm{-}\mathrm{2}\mathit{v}\mathrm{\Omega }\frac{\mathit{b}}{\mathit{a}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathit{ϵ}}\frac{\mathit{\partial \psi }}{\mathit{\partial x}}\mathit{,}\\ & & \frac{\mathit{\partial v}}{\mathit{\partial t}}\mathrm{+}\mathit{x}\frac{\mathit{\partial v}}{\mathit{\partial y}}\mathrm{-}\mathit{y}\frac{\mathit{\partial v}}{\mathit{\partial x}}\mathrm{+}\mathit{u}\mathrm{+}\mathrm{2}\mathit{u}\mathrm{\Omega }\frac{\mathit{a}}{\mathit{b}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}\mathit{ϵ}}\frac{\mathit{\partial \psi }}{\mathit{\partial y}}\mathit{,}\\ & & \frac{\mathit{\partial w}}{\mathit{\partial t}}\mathrm{+}\mathit{x}\frac{\mathit{\partial w}}{\mathit{\partial y}}\mathrm{-}\mathit{y}\frac{\mathit{\partial w}}{\mathit{\partial x}}\mathrm{=}\mathrm{-}\frac{\mathit{\partial \psi }}{\mathit{\partial z}}\\ & & \end{array}$

with ϵ = (a² − b²)/(a² + b²) being the ellipticity of the boundaries in the x,y-plane and Ω = Ω_P/Ω_F.

We now switch from Cartesian coordinates (x,y,z) to cylindrical coordinates (s,φ,z). In cylindrical coordinates, the above system of equations can be expressed as ${M}\left(\frac{\mathit{\partial }{u}}{\mathit{\partial t}}\mathrm{+}{H}{u}\right)\mathrm{+}\mathrm{2}\mathrm{\Omega }\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\Lambda }{u}\mathrm{=}\mathrm{-}\mathrm{\nabla }\mathit{\psi ,} \mathrm{\nabla }\mathrm{·}\mathrm{(}{\mathit{\rho }}_{\mathrm{0}}{u}\mathrm{)}\mathrm{=}\mathrm{0}$(16)with

$\begin{array}{ccc}& & {M}\mathrm{=}{I}\mathrm{+}\mathit{ϵ}{T}\mathit{,} {H}\mathrm{=}{I}\frac{\mathit{\partial }}{\mathit{\partial \phi }}\mathrm{+}\mathrm{2}{\Lambda }\mathit{,}\\ & & {T}\mathrm{=}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{,} {\Lambda }\mathrm{=}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathit{,}\end{array}$where I is the identity matrix. It matters at this stage that we restricted c to c = R because otherwise, the equation M = I + ϵT would contain additional terms.

We now use a Galerkin method. We seek solutions to Eq. (16) in the form $\begin{array}{ccc}{u}& \mathrm{=}& \sum _{\mathit{j}}{\mathit{C}}_{\mathit{j}}{{u}}_{\mathit{j}}{\mathrm{e}}^{\mathrm{i}\mathit{ft}}\\ \mathit{\psi }& \mathrm{=}& \sum _{\mathit{j}}{\mathit{C}}_{\mathit{j}}{\mathit{\psi }}_{\mathit{j}}{\mathrm{e}}^{\mathrm{i}\mathit{ft}}\mathit{,}\end{array}$(17)where the u_j are solutions of the unperturbed eigenvalue problem (ϵ = 0) $\begin{array}{ccc}& & {H}{{u}}_{\mathit{j}}\mathrm{+}\mathrm{\nabla }{\mathit{\psi }}_{\mathit{j}}\mathrm{+}\mathrm{2}{\Lambda }{{u}}_{\mathit{j}}\mathrm{\Omega }\mathrm{=}\mathrm{-}\mathrm{i}{\mathit{\omega }}_{\mathit{j}}{{u}}_{\mathit{j}}\mathit{,} \mathrm{\nabla }\mathrm{·}\mathrm{\left(}{\mathit{\rho }}_{\mathrm{0}}{{u}}_{\mathit{j}}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ & & {n̂}\mathrm{·}{{u}}_{\mathit{j}}\mathrm{=}\mathrm{0} \mathrm{ontheunitsphere}\mathrm{.}\end{array}$(18)This is the equation for inertial modes in a sphere. j is a proxy for the indices (n_j,m_j,k_j) which characterize an inertial mode in a sphere (Greenspan 1968). The variables n, m, and k will be indexed by j if it is necessary to distinguish different modes and will appear without an index otherwise. k is the azimuthal wavenumber, and n is the spatial degree or latitudinal wavenumber (on the surface of the sphere, the pressure distribution is given by the spherical harmonic ${\mathit{Y}}_{\mathit{n}}^{\mathit{k}}$). The index m numbers the eigenvalues for any fixed n and k and varies over a finite range. These equations can be solved analytically in the incompressible case (Zhang et al. 2001). They are still separable within the anelastic approximation provided that the density is of the form ${\mathit{\rho }}_{\mathrm{0}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{s}}^{\mathrm{2}}\mathrm{-}{\mathit{z}}^{\mathrm{2}}{\mathrm{)}}^{\mathit{\beta }}\mathit{.}$(19)Repeating the steps of Wu (2005a), one arrives at an eigenvalue problem for ψ which is written as a product in the form ψ(x₁,x₂,φ,t) = ψ₁(x₁)ψ₂(x₂)ζ(φ,t). (x₁,x₂,φ) are the ellipsoidal coordinates used in Wu (2005a) and introduced by Bryan (1889). One must be careful to distinguish the index which characterizes the different inertial modes from the index used in Eqs. (20) and (21) which stands for the two coordinates x₁ and x₂. In this paragraph we omit the index which characterizes the different inertial modes. The range of the three coordinates is x₁ ∈  [μ,1], x₂ ∈  [−μ,μ], with μ = −(ω + k)/(2(1 + Ω)) and φ is the azimuthal angle so that φ ∈  [0,2π]. μ is half the frequency of the inertial mode with respect to the frame rotating with Ω_P + Ω_F and ω is the frequency with respect to the inertial frame. All variables X (pressure and velocity components) depend on time t and on the azimuthal angle φ through the relation X ∝ ζ(φ,t) = exp [i(kφ − 2μ(Ω_P + Ω_F)t)]. We adopt the convention that μ ≥ 0 with k > 0 representing a prograde mode and k < 0 a retrograde mode. Modes with denotation (k,μ) and (−k, −μ) are physically the same modes so that one can restrict either μ or k to positive numbers and avoid redundancy. It is numerically more accurate to solve for the g_i defined by ${\mathit{g}}_{\mathit{i}}\mathrm{=}{\mathit{\psi }}_{\mathit{i}}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{\left|}\mathit{k}\mathrm{\right|}\mathit{/}\mathrm{2}}$, i = 1,2. The eigenvalue problem for the eigenvalue l and the eigenfrequency μ reads $\begin{array}{ccc}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}\frac{{\mathrm{d}}^{\mathrm{2}}{\mathit{g}}_{\mathit{i}}}{{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}{\mathit{x}}_{\mathit{i}}\mathrm{\left(}\mathrm{\right|}\mathit{k}\mathrm{|}& \mathrm{+}& \mathrm{1}\mathrm{)}\frac{\mathrm{d}{\mathit{g}}_{\mathit{i}}}{\mathrm{d}{\mathit{x}}_{\mathit{i}}}\mathrm{+}\frac{\mathrm{2}\mathit{\beta }{\mathit{x}}_{\mathit{i}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}\mathrm{)}}{{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\frac{\mathrm{d}{\mathit{g}}_{\mathit{i}}}{\mathrm{d}{\mathit{x}}_{\mathit{i}}}\\ & \mathrm{+}& \left[{\mathit{\lambda }}^{\mathrm{2}}\mathrm{-}\frac{\mathrm{2}\mathit{\beta }\mathrm{\left|}\mathit{k}\mathrm{\right|}{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}}{{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{2}\mathit{\beta k\mu }}{{\mathit{x}}_{\mathit{i}}^{\mathrm{2}}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\right]{\mathit{g}}_{\mathit{i}}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(20)with the boundary conditions

$\begin{array}{ccc}& & \frac{\mathrm{d}{\mathit{g}}_{\mathrm{1}}}{\mathrm{d}{\mathit{x}}_{\mathrm{1}}}{|}_{{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}\mathrm{=}\frac{{\mathit{\lambda }}^{\mathrm{2}}\mathrm{+}\mathrm{2}\mathit{\beta }\mathrm{\left[}\mathrm{\right(}\mathit{k\mu }\mathrm{-}\mathrm{\left|}\mathit{k}\mathrm{\right|}\mathrm{)}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}}{\mathrm{2}\mathrm{\left(}\mathrm{\right|}\mathit{k}\mathrm{|}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{g}}_{\mathrm{1}}{\mathrm{|}}_{{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}\\ & & \frac{\mathrm{d}{\mathit{g}}_{\mathrm{1}}}{\mathrm{d}{\mathit{x}}_{\mathrm{1}}}{|}_{{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\mu }}\mathrm{=}\mathrm{-}\frac{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{|}\mathit{k}\mathrm{\left|}\mathit{\mu }\mathrm{\right)}}{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{g}}_{\mathrm{1}}{\mathrm{|}}_{{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\mu }}\\ & & \frac{\mathrm{d}{\mathit{g}}_{\mathrm{2}}}{\mathrm{d}{\mathit{x}}_{\mathrm{2}}}{|}_{\mathrm{\left|}{\mathit{x}}_{\mathrm{2}}\mathrm{\right|}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\mu }}\mathrm{=}\mathrm{-}\mathrm{sgn}\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\frac{\mathrm{(}\mathit{k}\mathrm{-}\mathrm{|}\mathit{k}\mathrm{\left|}\mathit{\mu }\mathrm{\right)}}{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{g}}_{\mathrm{2}}{\mathrm{|}}_{\mathrm{\left|}{\mathit{x}}_{\mathrm{2}}\mathrm{\right|}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\mu }}\end{array}$(21)and λ² = l(l + 1) − |k|(| k| + 1). In general, both μ and l have to be determined numerically, but for β = 0, l is simply given by l = n.

Once the eigenfunctions are found, we define the scalar product $\mathrm{\left(}{{u}}_{\mathit{j}}\mathit{,}{{u}}_{\mathit{i}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{V}}{{u}}_{\mathit{j}}{{u}}_{\mathit{i}}^{\mathrm{\ast }}{\mathit{\rho }}_{\mathrm{0}}\mathit{r}\mathrm{d}\mathit{r}\mathrm{d}\mathit{\phi }\mathrm{d}\mathit{z}$(22)which differs from the usual scalar product by the factor ρ₀. The eigenvalues calculated in the anelastic approximation are orthogonal with this scalar product. This can be proven along the same lines as the proof for the orthogonality of inertial modes in the incompressible case (Greenspan 1968).

We substitute (17) in (16) and multiply this equation by ρ₀u_j. We obtain the system

$\begin{array}{ccc}& & \mathrm{(}\mathit{f}\mathrm{-}{\mathit{\omega }}_{\mathit{j}}\mathrm{)}{\mathit{C}}_{\mathit{j}}{\mathit{N}}_{\mathit{j}}^{\mathrm{2}}\mathrm{+}\mathit{ϵ}\mathrm{(}\mathit{f}\mathrm{+}{\mathit{k}}_{\mathit{j}}\mathrm{)}\sum _{\mathit{i}}{\mathit{V}}_{\mathit{j,i}}{\mathit{C}}_{\mathit{i}}\mathrm{=}\mathrm{0}\mathit{,}\\ & & \end{array}$(23)We assume the u_j to be normalized so that ${\mathit{N}}_{\mathit{j}}^{\mathrm{2}}\mathrm{=}\mathrm{1}$. We proceed on the assumption that V_j,i ≠ 0 only for n_j = n_i and k_j = k_i ± 2. The latter condition can be proven in general by noting that the dependence of the variables on the azimuthal angle is given by exp(ikφ). The first condition is less obvious. It was proven analytically by Kerswell (1993) for incompressible fluids. This proof can not simply be adapted to the compressible case because no analytical expressions exist for the modes in this case. However, we calculated the V_j,i numerically for all mode combinations with indices up to n = 20 and | k| = 10. In all cases, the results appeared to converge to zero as the spatial resolution used in the integration was improved.

Based on the above, we can infer that the elliptical instability comes from an interaction of the modes with azimuthal wavenumber k and k + 2, whereas the spatial degree n must be the same for the interacting modes. Modes in resonance will be written as (k,k + 2) and the spatial degree n will be stated separately.

Application of perturbation theory to (23) with the small parameter ϵ leads to the inviscid growth rate σ_inv${\mathit{\sigma }}_{\mathrm{inv}}^{\mathrm{2}}\mathrm{=}\mathrm{-}\frac{\mathrm{\left(}{\mathit{ϵ}}^{\mathrm{2}}{\mathit{V}}_{\mathit{i,j}}^{\mathrm{2}}\mathrm{\right(}{\mathit{k}}_{\mathit{i}}\mathrm{-}{\mathit{k}}_{\mathit{j}}\mathrm{)}\mathrm{+}\mathit{\delta \omega }{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}{\mathit{ϵ}}^{\mathrm{2}}{\mathit{V}}_{\mathit{i,j}}^{\mathrm{2}}{\mathit{q}}_{\mathit{i}}{\mathit{q}}_{\mathit{j}}}{\mathrm{4}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{ϵ}}^{\mathrm{2}}{\mathit{V}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}$(24)with δω = ω_i − ω_j, ${\mathit{V}}_{\mathit{i,j}}^{\mathrm{2}}\mathrm{=}{\mathit{V}}_{\mathit{i,j}}\mathrm{·}{\mathit{V}}_{\mathit{j,i}}$, and q_j = −2μ_j. Eq. (24) is the same as Eq. (3.5) of Gledzer & Ponomarev (1992). It follows from this equation that instability is possible only if | δω| ≤ O(ϵ).

Viscosity is included heuristically in our calculations by adding a damping rate to σ_inv. Two different expressions are used, one for the growth rate with free slip boundaries, σ_fs, and one for no slip boundaries, σ_ns. The two approximations are ${\mathit{\sigma }}_{\mathrm{fs}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{inv}}\mathrm{-}{\mathit{g}}_{\mathrm{1}}\hspace{0.17em}\mathit{n}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{Ek}$(25)and ${\mathit{\sigma }}_{\mathrm{ns}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{inv}}\mathrm{-}{\mathit{g}}_{\mathrm{2}}\mathrm{|}\mathrm{1}\mathrm{+}\mathrm{\Omega }\mathrm{|}\sqrt{\mathit{Ek}}\mathrm{-}{\mathit{g}}_{\mathrm{1}}\mathit{n}\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{Ek,}$(26)where g₁ and g₂ are constants on the order of 1 which in this paper are always chosen as g₁ = g₂ = 1, n is the spatial degree of the two interacting modes and Ek is the Ekman number. Equation (26) contains two dissipative terms. The second one corresponds to dissipation in the bulk, whereas the first one is due to friction inside Ekman layers at the boundaries. The first term is usually computed in the frame of reference in which the boundaries are at rest and expressed in multiples of the rotation rate of that frame, yielding decay rates of inertial modes of the form ${\mathit{g}}_{\mathrm{2}}\sqrt{\mathit{Ek}}$ with g₂ typically between 0.1 and 1 (Greenspan 1968). Transformed to the frame of reference and the unit of time used in Eq. (14), the decay rate is ${\mathit{g}}_{\mathrm{2}}\mathrm{|}\mathrm{1}\mathrm{+}\mathrm{\Omega }\mathrm{|}\sqrt{\mathit{Ek}}$ with Ek = ν/(| Ω_F + Ω_P|R²) and ν the viscosity. In many standard applications of rotating fluid mechanics, n and Ek are both small enough for the term in $\sqrt{\mathit{Ek}}$ to dominate. For free-slip boundary conditions on the other hand, only the bulk dissipation must be considered. We use n(n + 1) in the bulk damping term, because the velocity field of an inertial mode contains on the surface only the spherical harmonic ${\mathit{Y}}_{\mathit{n}}^{\mathit{k}}$. Assuming this is a reasonable approximation to the inertial mode at any radius r, and using ${\mathrm{\nabla }}^{\mathrm{2}}{\mathit{Y}}_{\mathit{n}}^{\mathit{k}}\mathrm{\left(}\mathrm{\Theta }\mathit{,\phi }\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{n}}{{\mathit{r}}^{\mathrm{2}}}{\mathit{Y}}_{\mathit{n}}^{\mathit{k}}\mathit{,}$(27)one justifies the formula for σ_fs (Lorenzani 2001).

3. Numerical implementation

We used numerical methods to calculate the growth rates according to Eq. (24). We first calculate the frequencies for β = 0 of inertial modes in a sphere by solving Eq. (29) from Wu (2005a)$\frac{\mathrm{d}{\mathit{P}}_{\mathit{n}}^{\mathit{k}}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}}{\mathrm{d}{\mathit{x}}_{\mathrm{1}}}{|}_{{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\mu }}\mathrm{=}\mathrm{-}\frac{\mathit{k}}{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{P}}_{\mathit{n}}^{\mathit{k}}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}{\mathrm{|}}_{{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathit{\mu }}$(28)by bisection. We use these eigenfrequencies, together with the corresponding eigenvalues l as starting values for a shooting method to solve Eq. (20). As noted in Wu (2005a), “the number of eigenmodes remains conserved when β varies, with a close one-to-one correspondence between modes in different density profiles”. So we increment β in steps of 0.1 and use the eigenvalues of the foregoing β as an initial point for the shooting method. The integration of the ODE was done by use of a Runge Kutta Cash-Karp method of fifth order with an adaptive stepsize control. Because the endpoints of the ODE for one variable are both singular points, we used the method of shooting to a fitting point (Press et al. 1992).

The next step is to calculate the integrals for Eq. (24). We only have to consider the modes which satisfy the resonance condition. Through a simple calculation the condition for the frequency detuning δω can be made more precise: Instability is possible only if | δω| < 4ϵ + 2ϵ².

Finally we want to solve the integrals V_i,j by Gaussian quadrature. It is not feasible to calculate these integrals in the ellipsoidal coordinate space (x₁,x₂,φ) because in these coordinates the limits of integration depend on the eigenfrequencies, and for the integration of V_i,j we consider two wave functions with in general two different eigenfrequencies. Therefore we perform the integration in (s,z,φ) space. The scalars ψ₁ and ψ₂ were obtained at several points x₁ and x₂ by the numerical integration of the ODE (20). The final integrations were performed using Gaussian quadrature with Gauss-Legendre abscissas and weights with a resolution of 50 points in both s and z directions. The integration in φ direction can be performed analytically. We know x₁(s,z) and x₂(s,z) and obtain ψ₁(x₁(s,z),x₂(s,z))ψ₂(x₁(s,z),x₂(s,z)). The scalars ψ₁ and ψ₂ were obtained at some points x₁ and x₂ which in general do not coincide with the points we need for the Gaussian quadrature. Therefore we need the scalars ψ₁,ψ₂ at arbitrary points x₁,x₂. This was accomplished by a cubic-spline interpolation.

4. Results

We start by investigating the effect of compressibility on individual inertial modes and the growth rate of triad resonances. We will first consider inviscid flows and include viscous effects at the end of this section.

Figure 3 shows the dispersion relation for k = −1 and k = 1 for incompressible fluids and for β = 18. As expected, compressibility modifies the inertial mode frequencies. There is a noteworthy distinction between the cases kq > 0 and kq < 0 in that the modes with the frequency q with the smallest absolute value at any given n always occurs for kq < 0. These low frequency modes will become relevant later on.

The mode frequencies change continuously with β so that individual modes can be tracked as a function of β. Figures 1 and 2 shows a few examples of triads with (k_i,k_j) = (−1,1). It can be seen that a chosen triad may lead to instability or not depending on β. A variation of β modifies both the structure and frequency of the inertial modes and hence the integrals in V_i,j and δω in Eq. (24). As expected from this equation, the positive growth rates in Figs. 1a and 2a are found around the β for which δω = 0. Figures 1b and 2b give the variation of δω with β and demonstrate the effect of compressibility on mode frequencies.

Fig. 1 a) Growth rates of modes (−1,1) for various n as a function of the compressibility parameter β for ϵ = 0.04, Ωp = 0 and Ek = 0. b) Corresponding δω.

Fig. 2 Same as Fig. 1 but for ϵ = 0.16.

Fig. 3 Frequencies of inertial modes in the co-rotating frame, q, for β = 0 (circles and squares) and β = 18 (x and plus) as a function of the spatial degree n. Only the modes with k = −1 (green, respectively plus and circles) and k = 1 (blue, respectively x and squares) are shown.

The case | k| = 1, n = 2 is exceptional because it corresponds to the purely toroidal motion in the spin-over mode. All toroidal modes have no radial velocity component, which implies for a radially dependent background density profile ρ₀ that ρ₀ drops from the continuity equation ∇·(ρ₀u) = 0. The density profile ρ₀ then disappears completely from the eigenvalue and stability problems, so that β does not affect the growth rate or the δω in the triad (−1,1) with n = 2. In addition, these two modes are the only ones in exact resonance with δω = 0 at β = 0.

Fig. 4 Maximum growth rates as a function of Ω for Ek = 0. For this figure we take into account all modes with | k| ≤ 10 and n ≤ 40. The (solid) black line in this figure is the growth rate given by Eq. (29).

It will now be argued that while the properties of inertial modes depend on compressibility, the stability of the flow as a whole does not. For an inviscid fluid, the two control parameters apart from β determining stability are the ellipticity of the deformation, ϵ, and the ratio of the rotation rates Ω defined below Eq. (24). The tidal flow is unstable if any two inertial modes form a triad with positive growth rate. Since there is an infinity of modes, a computational search for instability is necessarily restricted to a finite subset. Figure 4a and b is the result of a scan of all modes with | k| ≤ 10 and n ≤ 40. These two panels display the maximum growth rate found in this subset for the chosen ϵ as a function of Ω for different β. Inertial modes appear only within a limited band of frequencies, so that resonances are impossible in the interval −3/2 < Ω < −1/2 (Craik 1989). Figure 4 shows a part of the complementary interval by way of example. The various curves appear to have a common envelope σ_ud given by ${\mathit{\sigma }}_{\mathit{ud}}\mathrm{=}\frac{\mathrm{(}\mathrm{3}\mathrm{+}\mathrm{2}\mathrm{\Omega }{\mathrm{)}}^{\mathrm{2}}}{\mathrm{16}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{\Omega }{\mathrm{)}}^{\mathrm{2}}}\mathit{ϵ}\mathrm{·}$(29)This expression stems from the calculation of the inviscid growth rate in an incompressible, infinitely extended fluid (Miyazaki et al. 1995). Some deviations remain between Eq. (29) and the computed maximum growth rates in Fig. 4. However, these deviations become smaller and smaller if larger numbers of eigenmodes are included in the search for unstable triads. This can already be seen by comparing panels a and b of Fig. 4: Fig. 4b is for the smaller ϵ. For a given subset of allowed eigenmodes, fewer pairs of modes will meet the criterion | δω| ≤ O(ϵ) required by Eq. (24) for the smaller ϵ, so that larger deviations remain between the envelope and the computed maximum growth rates. It is plausible to assume that the deviations will entirely disappear if all inertial modes are taken into consideration.

We thus arrive at the picture that the stability limit does not depend on β, but β determines which modes are unstable. We will now compute the spatial degree n of the most unstable modes. A resonance can only occur if two modes i and j have frequencies such that ω_i − ω_j < O(ϵ). For a smaller ϵ, we expect for purely statistical reasons that a larger pool of eigenmodes is necessary to find any resonance and that the unstable modes have a higher n for smaller ϵ. Because ω_i = (1 + Ω)q_i − k_i for all i, the frequencies q_i and q_j have to obey ${\mathit{q}}_{\mathit{i}}\mathrm{-}{\mathit{q}}_{\mathit{j}}\mathrm{=}\frac{{\mathit{k}}_{\mathit{i}}\mathrm{-}{\mathit{k}}_{\mathit{j}}}{\mathrm{1}\mathrm{+}\mathrm{\Omega }}\mathrm{+}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{\Omega }}\mathit{O}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{·}$(30)In other words, if there is a mode at frequency q_j, it can only form an unstable triad if there is another mode with a frequency in an interval of size O(ϵ)/(1 + Ω) centered around q_j + (k_i − k_j)/(1 + Ω). Let us first consider the case of small | Ω |, so that q_i − q_j ≈ k_i − k_j + O(ϵ). All eigenmode frequencies obey | q_i| ≤ 2. Assuming that these frequencies are statistically uniformly distributed over the interval available for inertial modes, and that the matrix elements appearing in Eq. (24) do not systematically vary with frequency, one will find on average one resonance with positive growth rate among a number N of modes proportional to 1/ϵ. Only a finite number of modes exist with an n smaller than some limit L: for every n and k, there are n − |k | modes if k ≠ 0 and n − 1 modes if k = 0. For every n, the index k has to obey | k| ≤ n so that there are n² − 1 modes of spatial degree n and N(L) modes exist with n ≤ L, with N(L) given by: $\mathit{N}\mathrm{\left(}\mathit{L}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{L}}\mathrm{(}{\mathit{n}}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{6}}\mathit{L}\mathrm{(}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{-}\mathit{L}\mathit{.}$(31)For large L, one has N ∝ L³, and in order for these modes to contain a resonance, one needs L ∝ ϵ^−1/3, independently of Ω or β. The smaller the deformation ϵ, the higher the typical spatial complexity of the modes involved in resonances.

Fig. 5 Mean value of the spatial degree < n >, for the modes with a growth rate of at least 0.1σud (dashed lines) or 0.8σud (solid lines) as a function of ϵ for β = 0 (green x), β = 2.8 (red circles), and β = 18 (blue squares), together with fits according to Eq. (31).

Figure 5 verifies this scaling. This figure has been obtained as follows: At fixed ϵ, 20 equidistant points have been chosen in the interval of Ω ranging from −0.45 to 0.5. For every point, we determined the triad with the mode with the smallest n for which the growth rate is at least 0.1σ_ud and 0.8σ_ud. This minimal n, averaged over all points, is ⟨ n ⟩ and is shown in Fig. 5 as a function of ϵ. The averaging procedure is intended to remove the occasional outlier which occurs if a mode with a very low n happens to be part of a growing triad. The average ⟨ n ⟩ is rather a measure of the minimum n that is typically necessary to obtain a resonance with a given ϵ. As expected from the argument above, ⟨ n ⟩  ∝ ϵ^−1/3 for small ϵ.

We now turn to the case of large | Ω | in Eq. (30). For | Ω | tending to infinity, q_i − q_j tends to zero. At the same time, we see from Eq. (24) that positive growth rates are only possible for q_iq_j < 0, so that q_i and q_j both have to go to zero as 1/|Ω | if | Ω | goes to infinity. Figure 3 shows that n must be large enough to find any eigenvalue with an absolute value below some prescribed bound. This yields a minimal spatial degree n_min necessary to find a resonance at a given Ω. The corresponding mode frequency needs to lie in the appropriate interval of size O(ϵ)/(1 + Ω) for a positive growth rate. This additional condition is possibly met only for n larger than n_min.

We therefore need a relationship between the smallest frequencies and the spatial degrees of the modes. As already stated the frequencies of inertial modes in a uniform density sphere can be calculated according to Eq. (29) from Wu (2005a) , see Eq. (28).

For large n the asymptotic expansion (Abramowitz & Stegun 1972) $\begin{array}{ccc}{\mathit{P}}_{\mathit{n}}^{\mathit{k}}\mathrm{\left(}\cos \mathit{\theta }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{+}\mathit{k}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{\Gamma }\mathrm{(}\mathit{n}\mathrm{+}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}}{\left(\frac{\mathrm{1}}{\mathrm{2}}\mathit{\pi }\sin \mathit{\theta }\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\\ & & \mathrm{\times }\cos \left[\left(\mathit{n}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right)\mathit{\theta }\mathrm{-}\frac{\mathit{\pi }}{\mathrm{4}}\mathrm{+}\frac{\mathit{k\pi }}{\mathrm{2}}\right]\mathrm{+}\mathit{O}\mathrm{(}{\mathit{n}}^{-1}\mathrm{)}\end{array}$(32)holds. We set μ = cos(θ) and θ = π/2 + ϵ_w with ϵ_w small for μ small. A short calculation shows that for small μ and large n, Eq. (28) is equivalent to $\begin{array}{ccc}\mathrm{-}\mathit{k}\cos [\left(\mathit{n}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right){\mathit{ϵ}}_{\mathit{w}}& \mathrm{+}& \frac{\mathrm{(}\mathit{n}\mathrm{+}\mathit{k}\mathrm{)}\mathit{\pi }}{\mathrm{2}}]\\ & \mathrm{=}& \left(\mathit{n}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right)\sin \left[\left(\mathit{n}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right){\mathit{ϵ}}_{\mathit{w}}\mathrm{+}\frac{\mathrm{(}\mathit{n}\mathrm{+}{\mathit{k}}^{\mathrm{)}}\mathit{\pi }}{\mathrm{2}}\right]\mathit{.}\end{array}$(33)For n + k even this gives $\begin{array}{ccc}{\mathit{ϵ}}_{\mathit{w}}& \mathrm{\approx }\frac{\mathrm{-}\mathit{k}}{{\left(\mathit{n}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right)}^{\mathrm{2}}} \mathrm{for} \mathrm{(}\mathit{n}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}{\mathit{ϵ}}_{\mathit{w}}\mathrm{\to }\mathrm{0}& \end{array}$(34a)$\begin{array}{ccc}{\mathit{ϵ}}_{\mathit{w}}& \mathrm{\approx }\mathrm{\pm }\frac{\mathit{h\pi }}{\mathit{n}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}} \mathrm{for} \mathrm{(}\mathit{n}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}{\mathit{ϵ}}_{\mathit{w}}\mathrm{\to }\mathrm{\pm }\mathit{h\pi }& \end{array}$(34b)and for n + k odd ${\mathit{ϵ}}_{\mathit{w}}\mathrm{\approx }\mathrm{\pm }\frac{\mathit{h\pi }}{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{+}{\frac{\mathrm{1}}{\mathrm{2}}}^{\mathrm{)}}} \mathrm{for} \left(\mathit{n}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\right){\mathit{ϵ}}_{\mathit{w}}\mathrm{\to }\mathrm{\pm }\frac{\mathit{h\pi }}{\mathrm{2}}\mathit{,}$(35)where h is an integer. For large n and small k the μ with the smallest absolute value tends towards zero as k/(n + 1/2)². We can classify these modes as “slow modes”, because if the frequency of modes on any other branch of the dispersion relation in Fig. 3 tends to zero, it does so only in 1/(n + 1/2). We thus have to distinguish two cases when computing n_min: the growing resonance is either between two slow modes, or at least one of the two inertial modes involved is not a slow mode. In the latter case, n_min ∝ |Ω |, whereas in the former case, ${\mathit{n}}_{\min }\mathrm{\propto }\sqrt{\mathrm{\left|}\mathrm{\Omega }\mathrm{\right|}}$. However, since | k₁ − k₂| = 2, and choosing the indices such that k₂ > k₁, a resonance between slow modes can only occur if k₁ = −1 and k₂ = 1, and from (30) we can deduce that Ω must be negative. In summary, we expect n_min ∝ Ω for positive Ω, and for negative Ω, too, except when modes with frequencies of approximately k/(n + 1/2)² can resonate with each other.

Fig. 6 Minimum value of n, nmin as a function of Ω for various ϵ and positive Ω. a) β = 2.8 and β = 18 and b) β = 0.

Fig. 7 Same as in Fig. 6 but for negative Ω.

Fig. 8 Maximum growth rates as a function of Ω for Ek ≠ 0. For this figure we take into account all modes with a maximum | k| ≤ 10 and n ≤ 40. a) Calculated for free slip boundary conditions according to Eq. (25). b) Calculated for no slip boundary conditions according to Eq. (26). The (solid) black line in this figure is the growth rate given by Eq. (29). The (dashed without marker) cyan line in a) is the growth rate given by Eq. (36) with d2 = 3.35.

Figure 6 proves this scaling to be correct for positive Ω. We find n_min = d₁Ω, with a prefactor d₁ which depends on β. We can deduce from that figure that d₁ = 3, 5, and 10 for β = 0, 2.8, and 18, respectively. Some points lie above the line n_min = d₁Ω because some modes with spatial degree n_min, even though their frequency is small enough, do not find another mode to resonate with. For negative Ω (Fig. 7), one also finds points below n_min = d₁|Ω |. These correspond to resonances between two slow modes and always involve the azimuthal wavenumbers k₁ = −1 and k₂ = 1. Since Ω > 0 in most applications, we will use n_min = d₁Ω in the estimation of dissipation effects below.

We have up to here considered ideal fluids only. The remainder of this section will deal with viscous effects. Figure 8 shows the maximum growth rate if dissipation is heuristically taken into account according to Eqs. (25) and (26). The dissipation for no slip boundaries formally depends on n. However, at the Ek and n under consideration, the first term corresponding to friction at boundaries dominates and it is independent of mode structure. Therefore, the maximum growth rates in Fig. 4 are the same as before except for a downward shift, and are given by ${\mathit{\sigma }}_{\mathrm{ud}}\mathrm{-}{\mathit{g}}_{\mathrm{2}}\mathrm{|}\mathrm{1}\mathrm{+}\mathrm{\Omega }\mathrm{|}\sqrt{\mathit{Ek}}$. Instability is expected if ${\mathit{\sigma }}_{\mathrm{ud}}\mathrm{-}{\mathit{g}}_{\mathrm{2}}\mathrm{|}\mathrm{1}\mathrm{+}\mathrm{\Omega }\mathrm{|}\sqrt{\mathit{Ek}}\mathit{>}\mathrm{0}$, independently of β. Values of c depend on the inertial modes and may vary by an order of magnitude, but since the uncertainties on Ek are in practice much larger, there is little incentive to determine the prefactor c more accurately.

Free slip boundary conditions allow only bulk dissipation and are more complicated because dissipation now increases with n. This puts a natural limit to the spatial degrees n it is useful to consider, which restricts altogether the number of modes suitable for instabilities because | k| ≤ n and there is only a finite number of modes for a given k and n. If one includes more and more eigenvalues in the stability analysis, the curve of the maximum growth rate approaches a smooth curve of the same type as σ_ud up to some approximation, until additional modes suffer from such a high damping rate that they can not improve the maximum growth rate any more.

The precise stability limit now depends on the structure of the unstable modes and hence on β and Ω, but as an estimate, we can take the growth rate for free slip boundaries to be the inviscid growth rate reduced by the decay rate of modes with the expected spatial degree, for which we may reasonably take ⟨ n ⟩ for small | Ω | and n_min for large | Ω |, or as a formula interpolating between these two limits, n = d₁|Ω| + d₂ϵ^−1/3. For the parameters of Fig. 8, only the term in d₂ϵ^−1/3 is important and one arrives at ${\mathit{\sigma }}_{\mathrm{fs}}\mathrm{=}{\mathit{\sigma }}_{\mathrm{ud}}\mathrm{-}{\mathit{d}}_{\mathrm{2}}^{\mathrm{2}}{\mathit{ϵ}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\mathit{Ek,}$(36)which is again independent of β. From Fig. 8a, we deduce that d₂ is roughly 3.35.

5. Examples: Io’s tides on Jupiter, the binary system V636 Centauri and the Earth

The results from the previous section are used for three examples. We consider the 0.85  M_⊙ secondary as the perturbing object within the binary system V636 Centauri, the terrestrial core in the Earth Moon system, and the tides on Jupiter raised by Io.

The choice of the correct boundary conditions is generally ambiguous. For the Earth’s core it is clear that no-slip boundary conditions apply. The choice for the two remaining examples is not so clear. One finds in the literature arguments both in favor of no slip (Tassoul 1987, 1995; see also Tassoul & Tassoul 1997) as well as free slip boundary conditions (Rieutord 1992, 2008; see also Rieutord & Zahn 1997).

An even bigger problem is how to deal with additional fluid motion, such as convection, on which the tides are superimposed. The simplest approach is to ignore such motions, which is certainly justified if their amplitudes is much smaller than the amplitude of the tidal flow. Another approach consists in introducing a turbulent viscosity. However, turbulence modeling is always uncertain, and according to at least one well established turbulence model, elliptical instability is enhanced and not suppressed under some circumstances (Fabijonas & Holm 2003). When using a turbulent viscosity for the present problem, one has to deal with the unusual situation that one needs to compute the damping of a motion with a period which is typically shorter than the turn over time of the turbulent eddies. Several authors (Wu 2005b and in a similar form Ogilvie & Lin 2007) have used the following expression for the turbulent viscosity ν_t${\mathit{\nu }}_{\mathrm{t}}\mathrm{~}{\mathit{v}}_{\mathrm{cv}}{\mathit{l}}_{\mathrm{cv}}\frac{\mathrm{1}}{\mathrm{1}\mathrm{+}\mathrm{(}{\mathit{\omega }}_{\mathrm{tide}}{\mathit{\tau }}_{\mathrm{cv}}\mathit{/}\mathrm{(}\mathrm{2}\mathit{\pi }\mathrm{)}{\mathrm{)}}^{{\mathit{s}}_{\mathit{e}}}}\mathit{,}$(37)where v_cv, l_cv, and τ_cv are the characteristic convection velocity, mixing length, and turnover time and ω_tide = 2|Ω_F | is the frequency of the tidal forcing. s_e is a constant, generally s_e = 1 or s_e = 2 is used. Goldreich & Nicholson (1977) and Goodman & Oh (1997) suggested that s_e = 2, while Zahn (1977) argues for s_e = 1. For tides faster than τ_cv, one expects a turbulent viscosity proportional to v_cv and a mixing length given by the distance traveled by particles during one tidal period, l_cvT_f/τ_cv, with T_f = 2π/ω_tide, which corresponds to s_e = 1 in Eq. (37). Numerical simulations by Penev et al. (2009) confirm this expectation, but because of limited spatial resolution they were not able to simulate reliably the case of T_f ≪ τ_cv, the regime for which s_e = 2 has been proposed.

That s_e = 2 should be used at small tidal periods is in accordance with results from Ogilvie & Lesur (2012). Therfore we use s_e = 2 for the case T_f ≪ τ_cv and s_e = 1 for T_f ≳ τ_cv. In order to apply Eq. (37), one needs numbers for v_cv, l_cv, and τ_cv. The mixing length is estimated by l_cv ≈ H, with $\mathit{H}\mathrm{=}\mathrm{-}\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\ln \mathit{\rho }}\mathrm{=}\frac{{\mathit{r}}^{-1}\mathrm{-}\mathit{r}}{\mathrm{2}\mathit{\beta }}$ the density scale height. Therefore an order of magnitude estimation gives ${\mathit{l}}_{\mathrm{cv}}\mathrm{~}\frac{{\mathit{R}}_{\mathrm{c}}}{\mathrm{2}\mathit{\beta }}$. The convective velocity is approximated by v_cv ≈ (F/ρ)^1/3, with the energy flux $\mathit{F}\mathrm{=}\mathit{L}\mathit{/}\mathrm{\left(}\mathrm{4}\mathit{\pi }{\mathit{R}}_{\mathrm{c}}^{\mathrm{2}}\mathrm{\right)}$ and R_c the radius of the central body. An order of magnitude for v_cv is therefore obtained by using ${\mathit{v}}_{\mathrm{cv}}\mathrm{\approx }\mathrm{(}\mathrm{3}{\mathit{R}}_{\mathrm{c}}^{\mathrm{3}}\mathit{F}\mathit{/}\mathrm{(}\mathrm{4}{\mathit{m}}_{\mathrm{c}}\mathrm{)}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}$, with m_c the mass of the central body. Finally, τ_cv = l_cv/v_cv.

Table 1 lists for the three examples the growth rate for an inviscid fluid, and the damping rates computed for molecular and turbulent viscosities. In the no slip case the damping rate D is given by $\mathit{D}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{\Omega }\mathrm{)}\sqrt{\mathit{Ek}}$ according to Eq. (26). We neglect the volume damping term because Ω and Ek are small enough in both cases. In the free slip case, $\mathit{D}\mathrm{=}{\mathit{d}}_{\mathrm{2}}^{\mathrm{2}}{\mathit{ϵ}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}}\mathit{Ek}$ according to Eq. (36). For d₂ we choose 3.35, the value extracted from Fig. 8a. We neglect the terms which are relevant if Ω is large (see the last paragraph in Sect. 4) because Ω is small enough in both cases. The ellipticity is calculated according to $\mathit{ϵ}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{m}}_{\mathrm{p}}}{{\mathit{m}}_{\mathrm{c}}}\left(\frac{{\mathit{R}}_{\mathrm{c}}}{{\mathit{a}}_{\mathrm{cp}}}\right)3$(38)with m_p the mass of the perturbing body and a_cp the distance between the central and the perturbing body. The Ekman number is calculated according to $\mathit{Ek}\mathrm{=}\frac{\mathit{\nu }}{{\mathrm{\Omega }}_{\mathrm{spin}}{\mathit{R}}_{\mathrm{E}}^{\mathrm{2}}}$(39)with R_E = R_co the radius at the edge of the outer conducting region in the case of the Earth and the Jupiter-Io system and R_E = R_c in the case V636 Centauri and Ω_spin = Ω_F + Ω_P.

Instability is expected if the inviscid growth rate exceeds the viscous damping in Table 1. Table 2 summarizes the result of different modeling assumptions. In practice, the turbulence model decides on whether instability is predicted or not. Unfortunately, we have no turbulence model we can safely rely upon.

6. Conclusion and discussion

We computed the linear stability limit of tidal flow within the anelastic approximation through a perturbation calculation in which the small parameter is the deformation of the central body from spherical shape. The instabilities are described as superpositions of two inertial modes of the rotating sphere. The perturbation calculation is tractable if the initially three dimensional problem of the computation of the inertial modes is separable and reduces to the solution of an ordinary differential equation. This is the case if the density profile obeys a power law. For an inviscid fluid, the growth rate of the combination of two chosen modes depends on the density profile, but the growth rate maximized over all possible combinations does not. Furthermore, the maximum growth rate is the same as the one known for elliptical instability in infinitely extended, incompressible fluids. The fact that one finds the same maximum growth rates in so widely different situations justifies some tolerance towards modeling assumptions. For example, we restricted our analysis to density profiles obeying power laws and to the boundary condition that the velocity normal to the boundary be zero at the surface. After the calculations presented here, we expect that different modes will be found to be the most unstable for more realistic assumptions, but that the stability limit of the flow will stay the same. In the same fashion, we do not expect the particular choice of the major axis parallel to the rotation axis in Eq. (14) to affect the stability limit.

We added viscous effects empirically. If friction at a solid boundary dominates dissipation, the stability is easily determined as a function of the rotation rates of the central body and the tidal companion, the tidal deformation, and the Ekman number. If on the other hand bulk dissipation dominates, the dissipation depends on the flow structure of the unstable modes and therefore on all details of the model, in particular the density profile. However, dissipation depends mostly on the spatial degree which on average obeys a simple scaling law as a function of tidal deformation, so that an estimate of the stability limit is still possible.

Viscous damping also depends on boundary conditions, and models of an object’s interior help to decide on which approximation to the viscous damping is most accurate. For instance, the dissipation in Ekman layers usually dominates if a solid core is present. However, if one wants to predict the stability limit of a particular astrophysical object, the biggest uncertainty comes from a possible turbulent viscosity. Progress in the calculation of tidal dissipation thus mostly hinges on advances in the treatment of turbulence.

Acknowledgments

The authors acknowledges research funding by Deutsche Forschungsgemeinschaft (DFG) under grant SFB 963/1, project A5.

References

All Tables

All Figures

	Fig. 1 a) Growth rates of modes (−1,1) for various n as a function of the compressibility parameter β for ϵ = 0.04, Ωp = 0 and Ek = 0. b) Corresponding δω.
In the text

	Fig. 2 Same as Fig. 1 but for ϵ = 0.16.
In the text

	Fig. 3 Frequencies of inertial modes in the co-rotating frame, q, for β = 0 (circles and squares) and β = 18 (x and plus) as a function of the spatial degree n. Only the modes with k = −1 (green, respectively plus and circles) and k = 1 (blue, respectively x and squares) are shown.
In the text

	Fig. 4 Maximum growth rates as a function of Ω for Ek = 0. For this figure we take into account all modes with | k| ≤ 10 and n ≤ 40. The (solid) black line in this figure is the growth rate given by Eq. (29).
In the text

	Fig. 5 Mean value of the spatial degree < n >, for the modes with a growth rate of at least 0.1σud (dashed lines) or 0.8σud (solid lines) as a function of ϵ for β = 0 (green x), β = 2.8 (red circles), and β = 18 (blue squares), together with fits according to Eq. (31).
In the text

	Fig. 6 Minimum value of n, nmin as a function of Ω for various ϵ and positive Ω. a) β = 2.8 and β = 18 and b) β = 0.
In the text

	Fig. 7 Same as in Fig. 6 but for negative Ω.
In the text

Fig. 8 Maximum growth rates as a function of Ω for Ek ≠ 0. For this figure we take into account all modes with a maximum | k| ≤ 10 and n ≤ 40. a) Calculated for free slip boundary conditions according to Eq. (25). b) Calculated for no slip boundary conditions according to Eq. (26). The (solid) black line in this figure is the growth rate given by Eq. (29). The (dashed without marker) cyan line in a) is the growth rate given by Eq. (36) with d2 = 3.35.

In the text