Masahiro Fukuda

Nobuyoshi Koga

Loïc Favergeon

Kinetics and mechanisms of the thermal decomposition of copper(II) hydroxide under different water vapour pressures


B Favier
email: b.favier@damtp.cam.ac.uk

A J Barker

C Baruteau

G I Ogilvie

Nonlinear evolution of tidally forced inertial waves in rotating fluid bodies

Keywords:

We perform one of the first studies into the nonlinear evolution of tidally excited inertial waves in a uniformly rotating fluid body, exploring a simplified model of the fluid envelope of a planet (or the convective envelope of a solar-type star) subject to the gravitational tidal perturbations of an orbiting companion. Our model contains a perfectly rigid spherical core, which is surrounded by an envelope of incompressible uniform density fluid. The corresponding linear problem was studied in previous papers which this work extends into the nonlinear regime, at moderate Ekman numbers (the ratio of viscous to Coriolis accelerations). By performing high-resolution numerical simulations, using a combination of pseudo-spectral and spectral element methods, we investigate the effects of nonlinearities, which lead to time-dependence of the flow and the corresponding dissipation rate. Angular momentum is deposited non-uniformly, leading to the generation of significant differential rotation in the initially uniformly rotating fluid, i.e. the body does not evolve towards synchronism as a simple solid body rotator. This differential rotation modifies the properties of tidally excited inertial waves, changes the dissipative properties of the flow, and eventually becomes unstable to a secondary shear instability provided that the Ekman number is sufficiently small. Our main result is that the inclusion of nonlinearities eventually modifies the flow and the resulting dissipation from what linear calculations would predict, which has important implications for tidal dissipation in fluid bodies. We finally discuss some limitations of our simplified model, and propose avenues for future research to better understand the tidal evolution of rotating planets and stars.

INTRODUCTION

Understanding the gravitational tidal interactions between two orbiting bodies is an important unsolved problem in astrophysics. It is often the case that one or more of the bodies involved in the interaction is wholly or partly composed of fluid layers, such as a star or giant planet, an ice giant with a thick atmosphere, or a terrestrial planet with deep ocean. The processes by which the orbits and spins of these bodies evolve due to tidal interaction are poorly understood. The continuing discovery of short-period extrasolar planets makes a theory of tidal interactions very relevant, since they may have played an important role in producing the observed properties of these systems.

The importance of tidal interactions in a given system, quantified by the relevant timescale for tides to cause the orbit and spin of a body to evolve, is related to the rate at which tidal energy is dissipated. The mechanisms responsible for this dissipation depend strongly on the internal structure of the body. For fluid bodies, the hydrostatic response of a homogeneous fluid body to the varying gravitational potential of the companion is often referred as to the equilibrium tide (Darwin 1880). In the convective envelope of a gas giant planet or low-mass star, the equilibrium tide is thought to be dissipated by the turbulent motions driven by thermal convection, though there are significant uncertainties in the efficiency of this process (Zahn 1966;Goldreich & Nicholson 1977;Penev et al. 2007Penev et al. , 2009;;Ogilvie & Lesur 2012). In stably stratified fluid layers, such as the radiative layers of a planet or star, the resonant excitation of g-mode oscillations (or internal gravity waves) is thought to dominate the tidal dissipation (Zahn 1970;Goodman & Dickson 1998).

Since all astrophysical fluid bodies rotate, it is important to understand the effect that rotation has on the rates of tidal dissipation. This is particularly important because rotating flows support oscillatory motions called inertial waves [START_REF] Greenspan | The Theory of Rotating Fluids[END_REF]). These waves are restored by the Coriolis force and can be excited by low frequency tidal forcing when the absolute value of the forcing frequency (in a frame rotating with the fluid) is less than twice the spin frequency of the body. The importance of these waves at contributing to tidal dissipation has been emphasised in recent years (Ogilvie & Lin 2004;Wu 2005;Ogilvie & Lin 2007;[START_REF] Ivanov | [END_REF]Goodman & Lackner 2009;Ivanov & Papaloizou 2010;Papaloizou & Ivanov 2010). In nearly adiabatically stratified convective regions, internal gravity waves are not supported and the excitation and dissipation of inertial waves might play a dominant role.

In the case of a full sphere of incompressible fluid, a complete set of smooth inviscid inertial modes exist, which can be derived analytically. However, the case of a spherical shell is more relevant for the interiors of giant planets with solid cores, terrestrial planets with deep oceans, and, to some extent, the convective envelope of a solar-type star. In this geometry, regular inviscid inertial modes do not exist, and the viscous eigenfunctions correspond to rays propagating along the characteristics of the Poincaré wave equation. For some frequencies, these wave beams can converge toward limit cycles, known as wave attractors, as they alternately reflect from the inner solid core and the outer boundary (Rieutord & Valdettaro 1997;Rieutord et al. 2001). As a result, the viscous dissipation rate has a very sensitive dependence on the tidal forcing frequency (Ogilvie & Lin 2004;Ogilvie 2009;Riertord & Valdettaro 2010). Moreover, this dissipation is strongly enhanced for particular frequencies, hence it may play an important role in contributing to tidal evolution in rotating fluid bodies.

Most of the previous studies of inertial waves in a spherical shell have considered the linear behaviour of these singular inertial wave solutions. However, as the Ekman number (the ratio of viscous to Coriolis accelerations) is decreased towards the values that it is thought to take in planetary or stellar interiors, and as the amplitude of the tidal forcing is increased, nonlinearities will play an increasingly important role. When a wave attractor is present, it is known from experimental and theoretical studies that the inertial wave beam can become unstable and break, leading to a turbulent flow and a different mechanism of dissipation (Scolan et al. 2013;Jouve & Ogilvie 2013). The dissipation rate may be modified somewhat from that predicted by linear theory.

A related nonlinear mechanism that could contribute to tidal dissipation is the so-called elliptical instability (Pierrehumbert 1986;Bayly 1986;Waleffe 1990;Kerswell 2002). This is an instability of elliptical streamlines, such as those in the equilibrium tidal flow in a rotating fluid body. The elliptical instability leads to the excitation of inertial waves through a parametric resonance. The nonlinear outcome of this instability in isolation has been studied in laboratory experiments (e.g. Lacaze et al. 2004; Le Bars et al. 2007Bars et al. , 2010)), as well as local (Barker & Lithwick 2013a,b) and global numerical simulations (e.g. Cébron et al. 2010;Cébron et al. 2013). These works suggest that the elliptical instability could contribute to tidal dissipation at short orbital periods. However, the elliptical instability has yet to be studied in detail in spherical shell geometry (except briefly in Cébron et al. 2010).

In addition, nonlinearities could generate zonal (net azimuthal) flows, which might also alter the dissipative properties of the flow. An initial study into the effects of cylindrical differential rotation i.e. zonal flows, on inertial waves in a spherical shell has been undertaken by Baruteau & Rieutord (2013). They show that a differential rotation can modify the path of characteristics, and alter the frequency-dependence of the viscous dissipation. In the case of precessing or librating flows, it is well known that zonal flows can be driven by nonlinearities in the Ekman boundary layers (Busse 1968;Noir et al. 2001;Calkins et al. 2010;Sauret & Dizès 2013). The weakly nonlinear model of Tilgner (2007) also illustrated a mechanism by which nonlinearities of inertial waves in a spherical shell can generate zonal flows even in the absence of Ekman boundary layers. Intense axisymmetric flows were observed experimentally in the case of a rotating deformed sphere (Morize et al. 2010). Although the detailed mechanisms that generate them depend on the type of forcing and on the boundary conditions, the existence of zonal flows is a generic feature of forced rotating flows, hence it is important to study their effect on tidal dissipation in fluid bodies.

In this paper, we consider a simple model of a rotating fluid body subject to tidal forcing first studied in Ogilvie (2009). We explore the nonlinear behaviour of this model using high-resolution three-dimensional numerical simulations. To our knowledge, this is one of the first attempts to numerically compute forced inertial waves in a spherical shell as an initial value problem (see also Papaloizou & Ivanov 2010), and the first to study the nonlinear problem. Our model and the numerical methods adopted are presented in Section 2. We first illustrate the model and compare with previous results in the linear regime in Section 3. The main results concerning the nonlinear regime are discussed in Section 4. Finally, our conclusions and future directions are given in Section 5.

PHYSICAL MODEL

Model and governing equations

Our model is based on Ogilvie (2009). We suppose that the body is in uniform rotation with angular velocity Ω = Ωez where ez is the unit vector in the vertical direction. The fluid is incompressible with uniform density ρ, henceforth taken to be unity, without loss of generality. This assumption is adopted at this stage for simplicity. We neglect centrifugal distortion, effectively limiting us to studying bodies that rotate much slower than their dynamical frequency.

We define a spherical polar coordinate system centred on the body with coordinates (r, θ, φ) in the rotating frame. We aim to study the response of the fluid contained within a spherical shell with ri < r < re to a tidal gravitational

potential Ψ ∝ Aℜ Y m l (θ, φ)e -iωt . Y m l
is the spherical harmonic of degree l and order m with l m, ω is the tidal frequency in the rotating frame, related to that in the inertial frame ωi by ω = ωi -mΩ. A is an arbitrary real amplitude. Our system of units is defined such that ρ = 1, Ω = 1 and re = 1. The equations of motion are then

∂u ∂t + u • ∇u + 2ez × u = -∇W + E∇ 2 u (1) ∇ • u = 0 , (2)
where

E = ν Ωr 2 e (3)
is the Ekman number (the ratio of viscous to Coriolis accelerations), ν is kinematic viscosity of the fluid, assumed to be uniform. W = p+Ψ+Φ ′ , where p is the pressure, Ψ is the external tidal potential and Φ ′ is the internal self-gravitational potential perturbation. Φ ′ is an order-unity multiple of Ψ, it is in fact (3/2)Ψ for a homogeneous sphere, so it just amplifies the tidal potential and increases the amplitude of the tide.

The linear properties of inertial waves depend strongly on the size of the core (e.g. the dissipation rate appears to scale as the fifth power of the core size - Goodman & Lackner 2009;Ogilvie 2009). Throughout this paper, the core size is however fixed to ri = 0.5 in order to reduce the number of parameters. This core size is likely to limit the appearance of the elliptical instability, which is not the main focus of this work. This is because the nonlinear coupling between the equilibrium tide and the singular inertial waves in a spherical shell may be weaker than with the regular inertial modes of a core-free body. In addition, we will only consider moderate amplitude tidal forcing in this paper, with A ≪ 1. At these amplitudes, it is likely that the elliptical instability will be weak, since the maximum (inviscid) growth rate of the elliptical instability is proportional to the amplitude of the forcing (∼ A/Ω). Choosing such a small A is therefore likely to eliminate the elliptical instability from playing a major role in our calculations. It may also mean that it is difficult to observe any parametric subharmonic instabilities of the inertial wave beams, which are likely to limit the wave amplitudes when E → 0. We return to this point later in Section 4.

Boundary conditions

The outer boundary condition is used to drive the flow. Here we choose to impose a radial velocity on a rigid spherical boundary at r = re, which represents radial motions associated with rising and falling of the equilibrium tide at the surface of the body. This model is valid in the limit of frequencies low compared to surface gravity wave frequencies and small Ekman number, as demonstrated for the linear problem in Ogilvie (2009). In the case of a tidal potential dominated by the l = m = 2 spherical harmonic, the radial velocity at the outer boundary can be written as

ur(r = re) = 32π 15 Aℜ Y 2 2 (θ, φ)e -iωt , (4)
where A is an arbitrary (real) amplitude. Using this definition, the maximum radial velocity at the outer boundary is equal to unity when A = 1. Note that this choice of boundary condition is compatible with the incompressibility condition, as there is no net mass flux through the boundary.

In reality, as the angular momentum of the system evolves, the amplitude of the tidal forcing will also vary with time (as the fluid synchronises its spin with the tidal deformation, the radial displacement decreases). Our model does not capture this effect. However, since we do not intend to evolve the system until synchronism is reached (the spin typically evolves by only a few per-cent in our simulations), the differences between our results and those for which this effect is taken into account should be minor (we have verified that this is the case for a couple of example simulations in which this effect was taken into account.).

If the companion has an eccentric or inclined orbit, other spherical harmonics than the one present in equation ( 4) can be excited at various frequencies (e.g. Ogilvie & Lin 2004;Barker & Ogilvie 2009). However, here we focus on the particular component of the tidal potential as defined by equation ( 4), since it is usually the dominant one. In particular, this is the sole relevant component for the synchronisation problem, in which a companion on a circular orbit in the equatorial plane of the primary orbits with an orbital frequency that is not synchronous with the spin frequency of the fluid, so that ω = 0.

At the inner boundary, we neglect the elasticity of the core, and treat it as a rigid spherical boundary on which an impenetrability condition is imposed, so that ur(r = ri) = 0. The radius of the inner core strongly influences the excitation of inertial waves (e.g. Goodman & Lackner 2009;Ogilvie 2013). We did not change the geometry of the core in this paper, in order to focus on studying the dominant nonlinear effects. For numerical simplicity, we assume that both the inner core and outer boundary are stress-free. We therefore have er × σer = 0 at r = ri and r = re, where σ = ν ∇u + (∇u) T is the stress tensor and er is the unit vector in the radial direction. This is a good approximation for the outer boundary, but is only a simplification for the inner boundary, which is adopted mainly to avoid numerical constraints due to the generation of Ekman boundary layers. We later discuss the case of a no-slip inner core in section 4.6, since this is probably a better approximation to the boundary of the solid core of a giant planet.

In the case of a non-rotating body, the response of the fluid to this boundary forcing can be seen in Fig. 1. In a uniformly rotating body, inertial waves are generated at the critical latitude on the inner core (Goodman & Lackner 2009;Sauret 2013) whenever -2 < ω/Ω < 2.

Energy and angular momentum

To help with the analysis of our simulations, we define the volume-integrated kinetic energy, work done by the tidal forcing, and dissipation rate, by

K = 1 2 V |u| 2 dV, (5)
I = - V u • ∇W dV, (6)
D = -ν V u • ∇ 2 u dV. ( 7
)
The kinetic energy evolves in time according to

∂tK = I -D. (8)
We have checked that this is satisfied in the numerical simulations, indicating that our solutions are numerically converged. The total angular momentum, defined as

L = V x × u dV, (9)
is another important quantity in calculations with stress-free boundaries (e.g. Jones et al. (2011)). It evolves according to

∂L ∂t = -2 V x × (Ω × u) dV - V x × (u • ∇u) dV . ( 10
)
Note that with spherical boundaries, the pressure torque vanishes, and with stress-free conditions, so does the viscous torque. Note also that, if we write ũ = u + Ω × x, then the only contribution to the right hand side of equation ( 10) is

V x × ( ũ • ∇ ũ) dV = - ∂V (x × ũ) ũ • n dS = 0 , ( 11
)
where n is the normal to the outer boundary. This term is nonzero, in general, because the horizontal components of the velocity are nonzero at the boundaries, and a nonzero ur is imposed, i.e., this leads to an angular momentum flux through the outer boundary. Hence, this term is responsible for causing secular angular momentum evolution. This term is nonzero even for a linear problem, as long as Ω = 0. However, for a linear problem, the period-averaged contribution from this term vanishes and there is no net change in the angular momentum. Note however that linear theories can still be used to calculate the tidal torque, which is a quantity of second order in the tidal amplitude.

For linear calculations with l = m = 2 forcing, we correctly observe only small amplitude oscillatory behaviour in the components of L, with zero mean. However, in the nonlinear simulations, we observe net growth in vertical angular momentum |Lz|, as is consistent with our outer boundary condition, and represents the process of tidal synchronisation of the spin and orbit. We do not observe |Lx| and |Ly| to grow appreciably, which shows that numerical errors are playing a negligible role in the simulations (simulations with no-slip boundary conditions have an oscillatory Ekman layer, which exhibits oscillatory behaviour in Lx, Ly with zero mean, which we present in section 4.6).

Numerical methods

In this paper we use two different numerical approaches to solving equations (1)-( 2), which we will now briefly describe. Since the two numerical methods are very different, a careful comparison is required, in particular for the nonlinear solutions where previously published results are unavailable.

PARODY

In the following, we adapt and use the code PARODY1 to solve Eqs. (1)-(2). This code was originally written by E. Dormy [START_REF] Dormy | [END_REF]) and later improved by J. Aubert (Aubert et al. 2008). PARODY has been benchmarked against other numerical codes in the context of a convectively-driven dynamo problem (Christensen et al. 2001). The code is parallelised using both OpenMP and MPI. The time stepping is achieved using a mix of semiimplicit Crank-Nicholson scheme for the linear terms and a second-order Adams-Bashforth scheme for the nonlinear terms.

The velocity field is written using a poloidal-toroidal decomposition, thus ensuring incompressibility, with

u = ∇ × ∇ × (Ser) + ∇ × (T er) , ( 12
)
where T is the toroidal component and S is the poloidal component, and er is the unit vector in the radial direction.

Each of these scalars is decomposed onto spherical harmonics,

S = S m l (r)Y m l (θ, φ) and T = T m l (r)Y m l (θ, φ) , (13)
where the sum is carried over integers such that 0 m l. The radial functions T m l (r) and S m l (r) are represented by their discretized values on a non-uniform radial grid between the inner core located at r = ri and the outer core located at r = re. The radial derivatives are computed using second order finite-differences. The grid is denser close to the inner and outer boundaries in order to appropriately resolve flows near the boundaries. This mesh refinement at the boundaries is particularly important when using no-slip boundary conditions in order to accurately resolve thin Ekman boundary layers.

Due to the poloidal-toroidal decomposition ( 12), the implementation of a non-vanishing radial velocity at one of the boundaries is non-trivial. The stress-free condition imposes

∂ (T m l /r) ∂r = 0 (instead of T m l = 0 in the no-slip case) (14)
for all l, m at both r = ri and r = re. We focus here on a radial forcing corresponding to l = 2 and m = 2 as defined by equation ( 4). The boundary condition for the poloidal component at the outer radius re is in this case

∂ 2 S m l ∂r 2 = -A 16 r 2 e 2π 15 S m l if l = m = 2, ( 15
)
or zero otherwise, and

S m l = A 2 3 2π 15 ree iωt if l = m = 2, ( 16
)
or zero otherwise. The inner core is stress-free and impenetrable, with

∂ 2 S m l ∂r 2 = S m l = 0 (17)
for all l, m at r = ri.

The main advantage of PARODY is the spectral decomposition in the azimuthal direction, which allows us to focus on a limited range of azimuthal modes. In the linear case, one can take advantage of the fact that only modes with m = 2 are excited, effectively reducing the problem to two dimensions. The numerical convergence of the code has been checked for various Ekman numbers and forcing frequencies, and the typical resolution in the meridional plane used for production runs is shown in Table 1. In the nonlinear case, we also checked the numerical convergence with the number of azimuthal modes and we typically use mmax = 64 (which is sufficient for amplitudes A 10 -2 ).

Nek5000

For comparison with the nonlinear solutions obtained using PARODY, we also use the efficiently parallelised spectral element code Nek5000, written by Paul Fischer and collaborators (Fischer et al. 2007;Fischer et al. 2008). Spectral element methods combine the geometric flexibility of finite element methods with the accuracy of spectral methods, and solve the weak variational form of the equations of motion, similarly to finite element methods [START_REF] Deville | High-Order Methods for Incompressible Fluid Flow[END_REF]. This method partitions the domain into a set of E non-overlapping sub-domains, called elements, whose union is the entire domain. Within each element the velocity components and the pressure are represented as tensor product Lagrange interpolation polynomials of order N and N -2, respectively, defined at the Gauss-Lobatto-Legendre and Gauss-Legendre points. Such a method has algebraic convergence with increasing E , but spectral (exponential) convergence with increasing N (for smooth solutions), with the total number of grid points in 3D being E N 3 .

Temporal discretisation is based on a semi-implicit formulation, where the nonlinear and Coriolis terms are treated explicitly and the viscous terms are treated implicitly. In particular, a 3rd order backward-difference formula is used for the viscous & pressure terms, and a 3rd order extrapolation is used for the explicit terms. Dealiasing is used, with the polynomial order N increased by a factor of 3/2 for the evaluation of non-linear (advective) terms. In our spherical shell computations, the points lie on spherical shells, to double precision. Typical resolutions adopted are E = 576 and N = 20 (30 for the nonlinear terms).

For the computations reported in this paper, PARODY tends to be somewhat more efficient, primarily due its ability to simulate a restricted range of azimuthal mode numbers m, which is found to be sufficient when A ≪ 1. However, Nek5000 is likely to be more efficient at very high resolution, since the Legendre transform requires global communication of all spectral coefficients, which is not required in spectral element methods. Both codes have been compared for both linear and nonlinear problems, as we will describe below.

LINEAR REGIME

In this section, we illustrate the model described in section 2 and compare some of our results with previous works in the linear regime. This is to illustrate both the transient phase before our solution converges to a steady state, and to check that our predictions for the steady state in the linear regime are consistent with previously published results. Numerically, the linear regime is recovered by effectively switching off the second term on the left-hand side of equation ( 1) and fixing the amplitude A to an arbitrary value equal to unity. In this section, we primarily use PAR-ODY, since this allows us to focus on the linear response of the m = 2 modes. PARODY and Nek5000 have been found to agree accurately for the linear problem for all frequencies compared when E = 10 -5 . Nek5000 was not used to explore smaller Ekman numbers because the mesh adopted is fully three-dimensional, whereas PARODY can exploit the m = 2 symmetry, which reduces the dimensionality of the linear problem. The typical resolution required for such linear simulations can be found in Table 1.

Illustration of the model

We first illustrate how this model behaves for a particular frequency ω/Ω = √ 2. This frequency was considered by Riertord & Valdettaro (2010), because the path of characteristics generated at the critical latitude converges towards a simple rectangular shape (see the right panel in Fig. 3). The group velocity of an inertial wave is proportional to its

E = 10 -5 E = 10 -6 E = 10 -7
Figure 3. Azimuthally-averaged velocity magnitude |u| is the meridional plane. The solution is symmetric with respect to the equator so that only the upper quadrant is shown. In all cases, black corresponds to zero whereas red corresponds to the largest value. The frequency is ω/Ω = √ 2 and the Ekman number is E = 10 -5 , 10 -6 and 10 -7 moving from left to right in the panels. The black straight line shows the location of the critical latitude defined by the polar angle θc = arccos(ω/2) at the inner core. The dotted white line is tangent to the inner core and corresponds to the direction of propagation of the inertial waves emitted at the critical latitude on the inner core. On the right, we show the path of characteristics initialised at the inner critical latitude. wavelength and is inclined at an angle arcsin(ω/2Ω) to the rotation axis, which is the angle at which the (inviscid) rays propagate. The critical latitude is the location where the inertial waves propagate tangentially to the boundary. We compute the response of the fluid for four different Ekman numbers from E = 10 -4 down to E = 10 -7 .

We plot in Fig. 2 the total kinetic energy and the total dissipation rate versus time. In all cases, the system reaches a steady-state after a significant fraction of a global viscous timescale. Note that while we reach a steady state for E = 10 -4 , E = 10 -5 and E = 10 -6 , transients are still present at the end of our E = 10 -7 simulation. For this particular frequency, the dissipation rate decreases (not linearly) as the Ekman number is decreased. As we will see in the next section, this is not the case for all frequencies. As already discussed in section 2.3, there is no net evolution in the vertical component of the angular momentum.

We show the azimuthally-averaged velocity magnitude |u| in the meridional plane in Fig. 3. This flow pattern corresponds to the steady state in each case.

As the Ekman number decreases, it becomes particularly clear that an inertial wave beam is emitted at the critical circle at the inner boundary defined by θ = arccos(ω/2Ω) and r = ri. The emergence of waves from the critical circle has been noticed before (Tilgner 1999;Ogilvie & Lin 2004;Ogilvie 2009). In the inviscid limit, the flow is singular there (Stewartson & Rickard 1969) and the solution is regularised by viscosity. The path of characteristics emerging from the inner critical latitude in the inviscid limit is shown in the rightmost panel in Fig. 3.

Comparison with previous works

In this section, we compare our results with Ogilvie (2009). Three close frequencies are compared: ω/Ω = 1.05, ω/Ω = 1.1 and ω/Ω = 1.15. While these frequencies are very similar, the dissipation rate was shown to crucially depend on the Ekman number in a very different way for each of these frequencies. More specifically, the dissipation was roughly independent of the Ekman number for ω/Ω = 1.1, whereas it increases as the Ekman number is decreased for ω/Ω = 1.05, with the opposite behaviour for ω/Ω = 1.15. We reproduce this result here using our initial value approach. The Ekman number is varied from E = 10 -2 down to E = 10 -8 and the aspect ratio of the spherical shell is again ri/re = 0.5. The simulations are run until a steady state is reached, and we then measure the corresponding dissipation rates. We show in Fig. 4 the results reproduced from Ogilvie (2009) superimposed with the results from our initial value approach. The agreement with the results from PARODY is excellent, with slight differences at the smallest Ekman number of 10 -8 , for which the high resolution used (see table 1) implies millions of iterations before reaching a steady-state. This comparison is to be considered as a consistency check only, as our initial value problem is much more numerically demanding than the associated direct linear calculation of a steady-state response.

Scalings with the Ekman number:

nonlinearities in the astrophysical regime Stewartson (1972) showed that the singularities associated with the solutions to the Poincaré equation should turn into oscillating shear layers through viscosity. The width of these regions depends on the mechanism that generates them and must scale with the Ekman number. Finding scalings of this kind is important because it allows us to extrapolate our results to the astrophysically relevant regime, in which the Ekman number takes extremely small values Most of the previous analytical studies were derived in the context of no-slip boundaries so that Ekman boundary layers scaling as E 1/2 are present, which are not expected in our problem with stress-free boundaries. However, we can numerically determine the scalings from the results of the linear calculations presented in Section 3.1 for ω/Ω = √ 2. From Fig. 3, it is clear that the width of the shear layers scales as a positive power of the Ekman number. In order to quantify this scaling, we plot in Fig. 5 the azimuthally averaged velocity amplitude along two lines inclined with an angle arcsin(ω/2Ω) = π/4 with respect to the vertical axis. One of the lines intersects the origin r = 0 and passes through the inner critical latitude whereas the second one intersects the z-axis at z = 0.4 (see embedded plot in Fig. 5). In each of these plots, the spatial coordinate and the amplitude are compensated by some power of the Ekman number.

As expected from previous analysis (Kerswell 1995), the radial width of the shear layer at the critical latitude on the -2 0 2 4

E 1/5 |u(r, θ = π/4)| φ E -2/5 (r -r i ) E = 3 × 10 -5 E = 10 -6 E = 3 × 10 -6 E = 10 -7 E = 3 ×
E 1/6 |u(s)| φ E -1/3 (s -s c )
Figure 5. Scalings with the Ekman number E. Top: Azimuthallyaveraged velocity amplitude along radius at fixed polar angle θ = π/4 which corresponds to the critical latitude. We focus on the region close to r = r i which corresponds to the critical latitude for ω/Ω = √ 2. The radius is normalised by E 2/5 whereas the amplitude is normalised by E -1/5 . Bottom: Azimuthally-averaged velocity amplitude along the blue line shown in the embedded panel. The spatial coordinate along this line s is centred around the second intersection with the path of characteristics shown in the embedded figure and normalised by E 1/3 . The amplitude is normalised by E -1/6 . inner boundary scales as E 2/5 while the width of the internal shear layers in the direction normal to the path of characteristics scales as E 1/3 . Note that E 1/4 is another expected scaling but our results do not allow us to distinguish between E 1/3 and E 1/4 (Rieutord et al. 2001). The amplitude of the flow at the critical latitude and in the shear layers is another important issue. The amplitude of the velocity at the critical latitude scales as E -1/5 whereas the amplitude of the velocity inside the shear layers scales as E -1/6 (see the scalings used in Fig. 5). Note that we could not check the dependence of these scalings on the Doppler-shifted frequency ω. We only manage to unambiguously determine the scalings for ω/Ω = √ 2 because the path of characteristics is spatially very simple in that particular case (see Fig. 3). For other frequencies, the proximity and intersection between different wave beams make the scaling analysis difficult, especially at the Ekman numbers available numerically.

We may use the scalings determined above to predict at what values of the input parameters, and at what location in the flow, nonlinearities are likely to become important first.

To do this, we can compare the wave velocity amplitude u with the phase velocity of the local wave packet ω/k, where k ∼ l -1 is the perpendicular wave number, and l is the width of the wave beam, and suppose that nonlinear terms become important when uk ω 1. This can be thought of as a dimensionless measure of the nonlinearity in the wave beams, and it is very likely that an inertial wave beam with an amplitude larger than this will become unstable and break (see Clark & Sutherland 2010;Scolan et al. 2013 for inertial gravity waves and Bordes et al. 2012;Jouve & Ogilvie 2013 for inertial waves). Similar estimates are also appropriate for internal gravity waves (e.g. Barker & Ogilvie 2010;Bourget et al. 2013).

Using the scalings determined above, the nonlinearity scales as E -3/5 at the critical latitude and E -1/2 , at most, in the internal shear layers/wave beams (taking the thinnest E 1/3 scaling). This suggests that in the astrophysical regime, in which E → 0, the nonlinearity in the vicinity of the critical latitude is probably the most important. We therefore expect this location to be the one with the dominant nonlinear interactions. Inertial waves launched from this location might undergo instabilities and break before they can reflect from the boundaries. If this occurs, this will very likely modify the frequency dependence of the dissipation rate (Goodman & Lackner 2009). However, this might be difficult to capture numerically, since these instabilities probably have much smaller scales than the primary wave beams, therefore these instabilities might be difficult to capture in our simulations. Nevertheless, the aforementioned scalings give some insight into when non-linearities should become prevalent. From now on, we simulate the effects of these nonlinearities directly.

NONLINEAR REGIME

We now move on to the main focus of this work, which is to study the effects of nonlinearities. The nonlinear term in equation ( 1) is now taken into account. This leads to additional numerical constraints. In particular, it is necessary to include a range of azimuthal wave numbers in PARODY, since the m = 2 symmetry is no longer preserved in the presence of nonlinear couplings. The typical resolution in the meridional plane must be also increased for the same reason. In addition, the Courant-Friedrichs-Lewy stability constraint requires a smaller time step. Consequently, it is not possible to numerically reach the very low Ekman numbers that were obtainable in the purely linear regime. Most of the nonlinear simulations presented in this section correspond to E 10 -5 . The next subsections are devoted to the effect of nonlinearities while varying some of the relevant parameters. Our aim is to determine how nonlinear effects modify the solution from the predictions of linear theory.

Illustration for three frequencies and comparison of the two codes

To validate the results of our numerical simulations in the nonlinear regime, we compare the results using both codes, for the following simulations. The Ekman number is fixed to be E = 10 -5 whereas the amplitude of the forcing is A = 3.86 × 10 -3 . We compare the three frequencies already discussed in Section 3.2, namely ω/Ω = 1.05, ω/Ω = 1.1 and ω/Ω = 1.15. For PARODY, the resolution used is 32 Fourier modes in the azimuthal direction, 194 Legendre polynomials and 480 grid-points in the radial direction. For Nek5000, the resolution is E = 576 and N = 20 (30 for the nonlinear terms). Note the advantage of using a spectral method in the azimuthal direction, nonlinear couplings are weak for this particular amplitude, so that a relatively small number of azimuthal modes are required in order to reach numerical convergence. The total kinetic energy and viscous dissipation rate are plotted versus time in Fig. 6. The solid lines correspond to the results obtained with PARODY whereas the symbols correspond to the results obtained with Nek5000.

The thin dotted lines correspond to the previous linear results. Note that in the case ω/Ω = 1.15, the dissipation is This transition occurs for amplitudes A = 10 -2 and A = 5×10 -2 at times Ωt ≈ 300 and Ωt ≈ 3000 respectively. Note that during the early evolution of the solution, δΩ scales like A 2 . much lower than for the two other frequencies and there are barely any differences between the linear and nonlinear solutions for this amplitude. The agreement between the two codes is excellent, even after 500 periods (O(10 4 ) time steps), bearing in mind that the two numerical schemes are based on very different approaches 2 . We can therefore explore the nonlinear behaviour of the model with confidence.

Increasing the amplitude

In this section, we fix ω/Ω = √ 2. This frequency was already considered in the linear regime in Section 3. We now consider the nonlinear regime by progressively increasing the amplitude of the forcing from A = 10 -4 to A = 10 -2 . The dissipation rate versus time is shown in Fig. 7 for the various amplitudes along with the purely linear result. Con- trary to the linear regime, no steady state is reached and the dissipation rate is now a time-dependent quantity, with variations of more than one order of magnitude. This is primarily because angular momentum is continuously injected through the outer boundary, as we discussed in Section 2.3, which is an effect not present in the linear calculations. This can be seen in Fig. 7, where we plot the increase in the mean rotation rate of the fluid versus time. In the frame initially rotating at a frequency Ω, we define the increase in the volume-averaged rotation rate of the fluid as

δΩ = 1 V u φ r sin θ dV , (18)
where u φ is the azimuthal velocity and V is the volume of the spherical shell. δΩ is initially zero, and evolves with time so that the total rotation rate of the fluid at a given time is Ω * (t) = Ω + δΩ(t). Contrary to the linear regime, the contributions from the forcing do not average out over one period so that angular momentum is injected into the fluid (if ω > 0, extracted otherwise, see Section 4.7). This process will continue until the fluid is spinning synchronously with the forcing, i.e. when ω = ωi -2(Ω + δΩ) = 0. Even if both the amplitude and the Ekman number in our simulations are typically larger than we would expect in stellar or planetary interiors, it is not possible computationally to run the simulation until complete synchronisation is achieved. By changing the spin frequency Ω + δΩ of the fluid, we continuously sweep across different values for the Doppler- Figure 9. Azimuthally averaged azimuthal velocity in a meridional plane at various times for various forcing frequencies. In all cases, the amplitude of the forcing is A = 10 -2 whereas the Ekman number is E = 10 -5 . Black corresponds to the minimum of the azimuthal velocity whereas red corresponds to the maximum. The black lines perpendicular to the inner core shows the location of the critical latitude at t = 0 whereas the white lines shows some of the ray paths emitted from that particular latitude assuming uniform rotation. This shows that tidal forcing generates cylindrical differential rotation in the body in the synchronisation process, and the fluid certainly does not spin up as a solid body.

shifted forcing frequency ω, so that the properties of the linear excitation of inertial waves varies with time. It is known from previous linear calculations that the dissipation rate at low Ekman number is a complicated function of the forcing frequency (Ogilvie 2009). By spinning up the fluid, it is therefore unsurprising to observe strong variation in the dissipation rate as time evolves. In addition, the generation of differential rotation in the bulk i.e. zonal flows, also changes the properties of the inertial wave response (Baruteau & Rieutord 2013). This is discussed further in section 4.3.

As can be seen in Fig. 8, the wave structure observed in Fig. 3 in the linear regime is still visible in the nonlinear regime, even after thousands of periods. In fact, there are barely any differences between the linear and nonlinear regimes when looking at the radial or tangent components of the velocity at the same time. As will be discussed in Section 4.3, the main difference between the linear and nonlinear regimes is related to mean zonal flows. This is partly because the total rotation rate Ω + δΩ does not vary significantly during this simulation (see the amplitudes in the bottom panel of Fig. 7). In addition, the internal shear layers or the flow near the critical latitude do not appear to be unstable to small-scale instabilities. Instead, the dominant nonlinearities in the regime probed by these simulations, appear to be associated with the generation of zonal flows, as we will now describe.

Generation of zonal flows

In all of our nonlinear simulations, we observe a non-uniform deposition of angular momentum in the fluid. Since this angular momentum injection (or extraction) by the forcing is not distributed homogeneously in the fluid, this leads to the generation of differential rotation in the form of zonal flows, which is a generic feature of all nonlinear simulations reported in this paper. We show in Fig. 9 an example of such zonal flows for various positive and negative Doppler-shifted frequencies ω. The azimuthal component of the velocity is azimuthally-averaged and plotted in the meridional plane. The results shown in Fig. 9 correspond to A = 10 -2 and E = 10 -5 in all cases, and are plotted at an arbitrary time in the range 500 < Ωt < 1300. We stress, however, that these results are time-dependent and the azimuthal flows evolve significantly during the course of a given simulation.

It is clear that the azimuthal flow significantly departs from a purely solid body rotation. The white lines in Fig. 9 correspond to the path of characteristics emitted at the inner critical latitude. Angular momentum seems to be preferentially deposited where the inertial waves reflect on the outer boundary (this is particularly visible for ω/Ω = -0.2, ω/Ω = -0.8 and ω/Ω = 1.87). This differential rotation affects the local properties of inertial waves, since the frequency ω, and therefore the local direction of propagation of inertial waves, now depends on the cylindrical radius. The azimuthal component of the velocity is indeed nearly vertically invariant in most of our simulations, corresponding to geostrophic flows with rotation constant on cylinders. The properties of small-amplitude inertial waves propagating in a differentially rotating incompressible fluid con-

-10 -4 0 10 -4 2 × 10 -4 0 0.2 0.4 0.6 0.8 1 u φ /A 2 Cylindrical radius s ω/Ω = 1 A = 10 -3 A = 3 × 10 -3 A = 10 -2 A = 3 × 10 -2 -2 × 10 -4 0 2 × 10 -4 4 × 10 -4 0 0.2 0.4 0.6 0.8 1 u φ /A 2 Cylindrical radius s ω/Ω = √ 2
Figure 10. Averaged zonal velocity for ω/Ω = 1 and ω/Ω = √ 2 and various amplitudes of the forcing. . corresponds to a vertical and azimuthal average. Both plots correspond to Ωt = 50. The zonal velocities are scaled with the square of the amplitude of the boundary forcing. The two vertical dotted lines correspond to the location of the reflections (see Fig. 3) on the outer boundary whereas the solid line corresponds to the critical latitude on the outer boundary.

tained in a spherical shell have been recently investigated by Baruteau & Rieutord (2013), who note that inertial waves can now propagate along curved paths since the Dopplershifted frequency is now a function of space.

It has been suggested that the interactions of m = 2 inertial modes in a spherical shell can generate a zonal flow (Tilgner 2007). As expected from the weakly nonlinear origin of these zonal flows, their amplitudes scale as the square of the amplitude of the forcing. This is confirmed by the results presented in Fig. 10. The azimuthal component of the velocity is averaged along the vertical and azimuthal directions and plotted against the cylindrical radius s = r sin θ. We plot the results for two different frequencies, ω/Ω = 1 and √ 2, at E = 10 -5 and various amplitudes. All results corresponds to the early time Ωt = 50. As the zonal velocities are scaled with A 2 , they all collapse onto approximately the same curve. Note that this scaling was already observed in Fig. 7 at early times. At later times, the complicated frequency dependence of the inertial wave excitation, as ob- served in the linear problem, makes these solutions depart from each other more strongly.

These results also give an indication as to where the nonlinearities are dominant. The case ω/Ω = 1 is peculiar, since the location where the waves reflect on the outer boundary is the critical latitude of the outer sphere, i.e. s = 0.87. In any case, the angular momentum seems to be dominantly deposited near to that radius. The case ω/Ω = √ 2 is more representative, since the locations of the wave reflection (see Fig. 3) are distinct from the outer critical latitude. The former are shown in Fig. 10 as dotted lines whereas the latter is shown as a solid line. Clearly, the wave reflections on the outer boundary close to the equator play a dominant role compared to the outer critical latitude.

The mechanism responsible for the zonal flows in our case is at odds with the no-slip librating case where it has been shown that the torque in the Ekman boundary layers and its eruption at the outer critical latitude are dominant and drive the zonal flow (Calkins et al. 2010). In that case, the amplitude of the zonal flow scales as the square of the amplitude of the forcing but is independent of the Ekman number. In our case however, there are no Ekman boundary layers since we adopt stress-free boundaries. We therefore do no expect the same mechanism to drive zonal winds. As already mentioned, the reflection of the waves generated at the inner critical latitude on the outer boundary seems to be responsible for the generation of the zonal flow. While our result depends strongly on the nature of wave reflection at the boundaries, it is known that these reflections are very similar between a solid wall and a free surface as long as the frequency is much smaller than the one of surface gravity waves (Phillips 1963). However, since we allow for a radial flow across the outer boundary, it is not clear to what extent our choice of boundary conditions influences the nonlinear behaviour of the waves. This can only be resolved by considering a more realistic ellipsoidal geometry or by changing the way the waves are forced. This is left for future work.

We now attempt to quantify the variation of the amplitude of the differential rotation as a function of the Ekman number. It is difficult to define the amplitude of the differen-tial rotation at low Ekman numbers, since there is generally no quasi-steady state. Nevertheless, it is still helpful to define the energy in the differential rotation as (Tilgner 2007),

E dr = 1 2 V u φ φ -δΩ r sin θ 2 dV , (19)
where u φ φ is the azimuthal average of the zonal velocity and δΩ is the rotation rate of the fluid at a particular time, as seen in the frame rotating at the rate Ω. In Fig. 11, we plot this quantity versus the Ekman number for A = 10 -2 and three different frequencies ω/Ω = 0.88, ω/Ω = 1 and ω/Ω = 1.1. At large enough Ekman numbers (i.e. E > 5 × 10 -4 ), the system can reach a nonlinear steady state allowing for an unambiguous measure of the amplitude of the zonal flows. As one decreases the Ekman number, the system becomes periodic or chaotic which explains the presence of error bars in Fig. 11. In any case, the energy in the differential rotation increases as a negative power of the Ekman number, as is also found in the numerical model of Tilgner (2007) and in the experiments of a coreless deformed sphere by Sauret (2013) and Sauret et al. (2013). The two straight lines in Fig. 11 correspond to the arbitrary scalings E -3/2 and E -1/2 , and are shown for illustration. Note that the amplitude of forcing of A = 0.01 that we are considering in Fig. 11 can be considered to roughly represent the amplitude of the synchronisation tide inside a hot Jupiter orbiting a solar-type star in a one day orbit (if the tidal frequency is comparable with the spin frequency). A naive extrapolation of our results (for any scaling exponent in the given range) inside a hot Jupiter (taking an approximate value of E ≈ 10 -18 ) would suggest that the energy in the differential rotation would far exceed that associated with the solid body rotational kinetic energy, i.e. very large shears would be predicted. However, it is very likely that shear instabilities would become important as the Ekman number is decreased (see section 4.5 below), thereby modifying the scaling behaviours and preventing such large shears from developing. In addition, angular momentum redistribution by turbulent convection or magnetic fields, which we have omitted, could limit the amplitude of these zonal flows. Finally, in order to show that the differential rotation is indeed the dominant nonlinear feature observed in our simulations, we further compare our nonlinear results with predictions obtained from direct linear calculations. In particular, we focus on the total viscous dissipation rate calculated as a function of the Doppler-shifted frequency of the forcing. During the course of our nonlinear simulations, we monitor the total rotation rate of the fluid Ω * = Ω + δΩ, the viscous dissipation rate and the energy in the differential rotation as defined by equation ( 19). In Fig. 12, we plot the dissipation rate as a function of the normalised frequency ω/Ω * . In the linear case, the rotation rate of the fluid is fixed and the curved is obtained by direct linear calculation for each value of ω/Ω * (using the same method as in section 3.2 and in Ogilvie 2009). In the nonlinear case, the total rotation rate and the viscous dissipation rate are changing with time. We show the results corresponding to three different initial frequencies: ω/Ω = -1.8, ω/Ω = 0.6 and ω/Ω = 1.275. In all cases, we observe an initial transient phase where the dissipation rapidly increases (as observed in Fig. 7 and Fig. 13). The dissipation then saturates close to the value predicted by linear theory. As time increases, the differential rotation builds up and the dissipation rate predicted by our nonlinear simulations departs from that predicted by linear theory at the same Doppler-shifted frequency. This indicates that the linear predictions are accurate in the early stage of the synchronisation process, but that the nonlinear path to synchronisation will be very different, for example due to the differential rotation driven by nonlinearities. Note however that for some frequencies, departure or similarities between linear and nonlinear predictions are observed irrespective of the differential rotation (see ω/Ω = 0.6 in Fig. 12 for example, where a relatively strong differential rotation does not result in a dissipation rate very different from linear predictions).

Varying the frequency

While it is not possible to consider every frequency within the range -2 < ω/Ω < 2 in the nonlinear case, we consider several frequencies while fixing the forcing amplitude to be A = 10 -2 . The viscous dissipation rates versus time are shown in Fig. 13. Note that for all frequencies considered here, no steady state is obtained for E = 10 -5 even after 5000 periods. For ω/Ω = 0.4, 0.6 and 0.8, we also plot the linear results as dotted lines for comparison. The dissipation in the nonlinear regime can be larger than in the linear regime (as it is the case for ω/Ω = 0.4), smaller (as it is the case for ω/Ω = 0.8) or both depending on time. The Figure 13. Dissipation rate (in units of ρr 3 e A 2 Ω) versus time for various positive frequencies, E = 10 -5 and A = 10 -2 . For ω/Ω = 0.4, 0.6 and 0.8, we also plot the linear results as dotted lines for comparison. complicated dependence of the viscous dissipation on frequency was already observed in the linear regime (Ogilvie 2009) and is also a property of the nonlinear regime. This is related to the fact that the dynamics is still strongly dominated by wave beams generated at the inner critical latitude and reflecting on the boundaries.

Secondary shear instabilities

We have observed the amplitude of the zonal flows to increase as we decrease the Ekman number in Fig. 11. This suggests that strong localised shear could be generated in the regime of very small Ekman numbers, which might be unstable to hydrodynamical shear instabilities. In this section we present the results from a particular simulation with A = 0.02, ω/Ω = √ 2 and E = 10 -5 . As already observed, we find an increase in the vertical component of the angular momentum, and a zonal flow is driven, primarily close to the location of the reflection on the outer boundary, which in this case in close to the equator. We show on the left four panels of Fig. 14 the amplitude of the velocity |u| is a meridional slice and in the equatorial plane at Ωt ≈ 118 and Ωt ≈ 137. We also show the time evolution of the total viscous dissipation on the bottom right panel of Fig. 14. At time Ωt ≈ 130, we observe a sudden increase in the total dissipation. This corresponds to the onset of a hydrodynamical shear instability, which preferentially excites modes with large azimuthal wave numbers m. The kinetic energy spectrum in the azimuthal direction is shown on the top right panel of Fig. 14. We observe the growth of large m components in the flow in both the spectrum and the equatorial slice. Large m components were not present during the early nonlinear evolution, where the slope is much steeper in spectral space, indicating only weak nonlinear transfers to modes with smaller azimuthal scales than m = 2. The axisymmetric m = 0 component of the flow continuously grows until energy is transferred into these non-axisymmetric modes, after which is appears to saturate. After the instability saturates, the energy at small azimuthal wavelength is dissipated, lead-ing to a steeper kinetic energy spectrum (see the spectrum at Ωt ≈ 300 in Fig. 14). Although we didn't pursue the simulation further, it is possible that once the turbulence is dissipated, another shear instability kicks in, leading to a cyclic behaviour between the laminar and turbulent states.

To our knowledge, this is the first numerical evidence of unstable zonal flows driven by inertial waves in spherical shells. Similar secondary shear instabilities have been recently observed experimentally in a rotating deformed sphere with or without an inner core (Sauret 2013;Sauret et al. 2013). In an astrophysical object, it is possible that magnetic stresses, magnetohydrodynamical shear instabilities, or convection could suppress this differential rotation before it could become unstable to such a hydrodynamical shear instability.

No-slip inner core

In this section we briefly present our results varying the inner boundary condition on the core. In particular, we choose to adopt no-slip, rather than stress-free conditions. This might be more relevant to terrestrial planets and the solid cores of a giant planet, as well as to laboratory experiments. In Fig. 15, we plot the normalised kinetic energy K and viscous dissipation rate D for a set of nonlinear simulations with A = 3.86 × 10 -3 , for the three frequencies ω/Ω = 1.05, 1.10 and 1.15. Again, we compare the results from both numerical methods described in Section 2.4. The kinetic energy evolves very similarly to the stress-free core case. There are some differences in D, primarily due to the presence of oscillatory Ekman boundary layers that cause periodic oscillations in Lx and Ly. The similarity between Fig. 6 and Fig. 15 however indicates that the generic nonlinear properties that we have observed in the rest of this paper seem to be robust to this change in the inner boundary condition.

We show the azimuthal velocity in a meridional plane for ω/Ω = 1.1 and 0.6 and both boundary conditions in Fig. 16. The time is the same in both cases, Ωt = 1225. The zonal flow is initially generated locally where the waves reflect from the outer boundary and then spreads in cylindrical radius due to viscosity. In the no-slip case, the zonal flow only exists outside the cylinder tangent to the inner core whereas it continues to spread in the stress-free case. Note however that for ω/Ω = 0.6, a weak zonal flow is generated within the tangent cylinder in both cases. This shows that a no-slip inner core has only a weak effect on the dissipation mechanism and generation of zonal flows for the parameters considered in this paper.

Limitation of the model

When the Doppler-shifted frequency ω is positive, one would expect the fluid to spin up as the forcing injects angular momentum into the system. This can be thought to represent the synchronisation of the spin of the body with the orbit of its companion when the companion is orbiting at a faster rate than the primary body is spinning, so that the tidal torque causes the angular momentum of the body to increase. The opposite is supposed to happen for ω < 0, which corresponds to the case of a companion orbiting at a slower rate than the spin of the body, so that the tidal torque acts to spin down the body. While this is true for most of the frequencies and amplitudes that we have investigated, some striking exceptions were observed. This is the case when ω/Ω = 0.75, for example. In this case, the angular momentum should increase, whereas the opposite is observed in our simulations, leading to a desynchronisation of the fluid, which is unexpected. This anomalous evolution of the angular momentum is also observed for ω/Ω = 1.87, see Fig. 9, where the vertical angular momentum is mostly negative. We compare the time evolution of the vertical component of the angular momentum Lz in Fig. 17. The results from both codes are shown for E = 10 -5 and A = 3.68 × 10 -3 . Note that this unexpected result persists for other Ekman numbers and forcing amplitudes. Although there is a slight disagreement between the codes at large times, the fact that the vertical component of the angular momentum is negative is a robust feature (the discrepancy between the codes appears to be caused by differences in the excitation of transients at t = 0, as the forcing is instantaneously switched on). This behaviour is observed for particular frequencies only. It is not a result of numerical errors, since both codes exhibit the same behaviour, and the spatial and temporal resolution of these results has been carefully checked for convergence. In addition, this does not appear to be only a transient phase in the simulations, and is observed to persist after several thousand rotation periods. It is therefore a surprising property of the model for our nonlinear simulations at certain specific frequencies.

This behaviour is ultimately due to the fact that we consider an open system with a non-vanishing radial velocity at the outer boundary. While we constrain the velocity of the fluid going in and out of the system to mimic the elliptical tidal deformation of the body, we do not constrain the ingoing (vertical) angular momentum flux. The angular momentum injected in the system is related to Reynolds stresses uru φ φ at the outer surface as described in Section 2.3. While the radial velocity is imposed at r = re, the azimuthal component is not constrained, since we adopt stressfree boundary conditions on a sphere instead of an ellipsoid. This seems to be insufficient to realistically constrain the angular momentum evolution in some cases. In addition, while we clearly observe a zonal flow driven by the wave reflection on the outer boundary, it is not clear to what extent this conclusion depends on our choice of boundary conditions. Recent experimental results seem to be qualitatively consistent with our conclusions (Sauret 2013;Sauret et al. 2013), but the mechanism responsible for the zonal flows ought to be clarified in a more realistic model, either using the realistic ellipsoidal geometry, or introducing a body force in order to avoid the unrealistic angular momentum source term at the outer boundary. These examples clearly illustrate the limitations of the current model. Nevertheless, the qualitative behaviour of the nonlinear system described by our simulations is very likely to be robust, as the similitude between our results and the ones reported by Tilgner (2007) and Sauret (2013); Sauret et al. (2013) suggests. This includes the generation of zonal flows, complicated time evolution of the dissipation rates with time and secondary shear instability, which are expected to all persist in a more realistic model.

DISCUSSION

In this paper, we have investigated numerically, for the first time, the nonlinear behaviour of tidally forced inertial waves in a spherical shell, as an initial value problem. The forcing corresponds to an imposed radial velocity at the outer boundary, mimicking the radial displacement the fluid would experience in response to a tidal gravitational potential. We have focused here on the dominant l = m = 2 spherical harmonic, which is usually the dominant one for astrophysical applications, such as tidal synchronisation. This intentionally simplified model is designed to be an idealised representation of the convective regions of a giant planet, the fluid regions of a lower mass planet such as a Neptune mass planet, or to the convective envelope of a solar-type star. This allows us to perform a detailed study of the nonlinear evolution of tidally forced inertial waves.

First, the purely linear regime is studied, and we recover previous results obtained using a direct solution of the steady-state response (Ogilvie 2009). In the nonlinear regime, the forcing injects or extracts angular momentum (depending on the sign of the Doppler-shifted frequency ω) in a non-uniform manner, leading to the generation of significant differential rotation in the interior of the body. As the fluid spins up (or down), the fluid experiences highly time-dependent dissipation rates, which are a complicated function of the forcing frequency. The amplitude of these zonal flows scales as A 2 E -γ where γ is some positive number depending on the frequency of the forcing. This suggests possible shear instabilities in the low Ekman number or large amplitude regime that may be expected in the interiors of short-period extrasolar planets. A hydrodynamical shear instability was indeed observed in some of our simulations, which is similar to an instability recently observed in laboratory experiments (Sauret 2013;Sauret et al. 2013). Note however that we cannot at this stage perform a quantitative comparison between experiments and numerical simulations as there are significant differences between the two, the most important being the different type of boundary conditions used (stress-free plus inflow on a sphere in our simulations versus no-slip on a deformed sphere in experiments).

We have so far considered the regime of moderate Ekman numbers, which is the regime that is possible to access numerically, as well as relatively low amplitudes, to try to understand the dominant nonlinearities. It is possible that this restriction constrains the nonlinearities to promote the importance of zonal flows. This restriction might also eliminate the possibility of the internal wave beams becoming unstable to small-scale parametric instabilities. At larger amplitudes, and more importantly, at lower Ekman numbers, parametric subharmonic instabilities, as recently observed by Bordes et al. (2012) in the laboratory, or nonlinearities at the critical latitude itself (Goodman & Lackner 2009), might become more important, and might wash out some of the frequency dependence of the dissipation rate. These possibilities will be explored in future studies. Nevertheless, the ultimate nonlinear outcome of these instabilities is likely to also result in the formation of zonal flows such as those that we have observed.

Due to the complexity of this problem and the limitations of our intentionally simplified model, we cannot currently compare our results directly to astrophysical observations. However, it is worth mentioning those aspects that are likely to be relevant for the astrophysical problem. When nonlinearities are considered, the flow can depart significantly from what would be predicted from linear theory, primarily due to the generation of differential rotation in the initially uniformly rotating body, in the form of zonal flows. These zonal flows can play a role in modifying the subsequent dissipation rate, which is what is important in relating our work to astrophysical observations. This is one example to illustrate that the simple picture assuming the body approaches synchronism as a uniformly rotating body is likely to be incorrect. This is also true when tidal forcing is able to excite internal gravity waves (Goldreich & Nicholson 1989;Barker & Ogilvie 2010). In all cases, the inclusion of nonlinearities and time evolution of the angular momentum leads to complicated time-dependence of the flow and the dissipation rate. There are large differences between the dissipation rate as the amplitude and frequency of the tidal forcing is varied, as well as the Ekman number of the body. There does not appear to be a general trend for the normalised dissipation rate: sometimes nonlinear terms increase the dissipation rate over the linear prediction, and they sometimes decrease it.

Tidal forcing is likely to lead to complicated (cylindrical) differential rotation in the interior of the body. Such differential rotation could be acted on by hydrodynamical or magnetohydrodynamical shear instabilities, or suppressed by magnetic stresses and convection. Those secondary processes could play an important role by contributing to tidal dissipation in fluid bodies. Future work will include studying more realistic mechanisms of tidal forcing, which do not exhibit the limitations of the current model discussed in Section 4.7, in order to more accurately capture the tidal forcing of inertial waves in fluid stars and planets. In addition, topics of further study include studying the influence of stratification and convection, the latter of which could play a role in dissipating inertial waves, and the internal structure of the body and the imperfect rigidity of the core.

Figure 1 .

 1 Figure 1. Response of the fluid due to the boundary forcing with A = 10 -2 for a non-rotating body. Top: Radial velocity in the equatorial plane. Bottom: Radial velocity in the meridional plane corresponding the dotted white line in the top figure. With our scaling, max(ur) = 1 when A = 1. The dominant azimuthal mode m = 2 and meridional mode l = 2 are clearly visible. As time evolves, this pattern rotates around the vertical axis at the frequency ω, following the motion of the companion.

Figure 2 .

 2 Figure 2. Kinetic energy (in units of ρr 3 e A 2 ) and dissipation rate (in units of ρr 3 e A 2 Ω) versus time for ω/Ω = √ 2 and various Ekman numbers. The time is expressed in units of the spin frequency Ω.

ForcedFigure 4 .

 4 Figure 4. Dissipation rate (in units of ρr 3 e A 2 Ω) for ω/Ω = 1.05, 1.1 and 1.15 and various Ekman numbers. Lines are reproduced from the Fig.7 of Ogilvie (2009) whereas symbols correspond the steady-state of our initial value problem solved with PARODY.

Figure 6 .

 6 Figure6. Comparison between the two codes for a simulation with A = 0.00386 and E = 10 -5 for the three frequencies 1.05, 1.10 and 1.15. The solid lines correspond to the results obtained with PARODY whereas the symbols correspond to the results obtained with Nek5000. The viscous dissipation rate D is shown in units of r 3 e A 2 Ω whereas the kinetic energy K is shown in units of r 3 e A 2 . The thin dotted lines correspond to the linear results already presented in Fig.4

Figure 7 .

 7 Figure 7. Dissipation rate (in units of ρr 3 e A 2 Ω) and spin frequency of the fluid δΩ in the rotating frame for ω/Ω = √ 2 and various amplitudes A of the forcing. The horizontal line shows the transition between two different regimes as soon as δΩ > 10 -2 .This transition occurs for amplitudes A = 10 -2 and A = 5×10 -2 at times Ωt ≈ 300 and Ωt ≈ 3000 respectively. Note that during the early evolution of the solution, δΩ scales like A 2 .

Figure 8 .

 8 Figure 8. Normalised tangential velocity u θ in a meridional slice for ω/Ω = √ 2 and E = 10 -5 . We compare the linear results (top) during the steady state and the nonlinear regime (bottom) with A = 5 × 10 -3 at the arbitrary time Ωt ≈ 3000, well into the nonlinear regime.

Figure 11 .

 11 Figure11. Energy in the differential rotation as defined by equation (19). The error bars correspond to the maximum and minimum values between Ωt = 10 3 and Ωt = 10 4 . A large Ekman numbers (i.e. E > 10 5 ), a nonlinear steady state is reached for Ωt > 10 3 which is why no error bars are present in that case.

Figure 12 .

 12 Figure12. Dissipation rate as a function of ω/Ω * , where Ω * = Ω + δΩ is the total rotation rate of the fluid. The linear prediction is shown as a black line whereas the time evolution obtained from nonlinear simulations is shown as a colored line, where the color corresponds to the amplitude of the differential rotation, as defined by equation (19). From top to bottom, the initial Dopplershifted frequencies are ω/Ω = -1.8, ω/Ω = 0.6 and ω/Ω = 1.275

Figure 14 .

 14 Figure 14. Secondary shear instability. We show the velocity amplitude in meridional and equatorial planes just before (Ωt ≈ 118) and after (Ωt ≈ 137) the instability. The kinetic energy spectrum as a function of the shifted azimuthal wave number m + 1 is shown on the top right. The thin dotted line corresponds to the slope m -5/3 . The time evolution of the dissipation rate is shown on the bottom right. The two vertical dotted lines correspond to the two times at which the flow is visualised on the left.

Figure 15 .

 15 Figure 15. Same as in Fig. 6 but with no-slip boundary condition at the inner core.

Figure 16 .

 16 Figure 16. Azimuthally-averaged zonal velocity for ω/Ω = 1.1 (top) and ω/Ω = 0.6 (bottom). The amplitude of the forcing is A = 10 -2 and E = 10 -5 . The boundary condition on the inner core is stress-free on the left and no-slip on the right. The poloidal component of the flow is very similar in both cases.

Figure 17 .

 17 Figure 17. Vertical component of the angular momentum for E = 10 -5 , A = 10 -2 and ω/Ω = 0.75. Although the forcing frequency is positive, we observe a decrease in the spin frequency of the fluid. We show the results from both numerical methods.

Table 1 .

 1 Typical resolutions used within PARODY in the linear regime. N R is the radial resolution whereas N L is the maximum degree of the Legendre polynomials. For the linear calculations we only include one azimuthal mode number, m = 2.

	Ekman number	N R	N L
	10 -4	240	32
	10 -6	960	256
	10 -8	3840	1024

http://www.ipgp.fr/ ~aubert/DMFI.html

The most likely source of discrepancy between the two codes is that the incompressibility condition is exactly satisfied in PAR-ODY, as a result of the poloidal-toroidal decomposition of the velocity field. On the other hand, Nek5000 solves a discrete Poisson equation for pressure, which it solves by preconditioned conjugate gradient iteration with a given tolerance, which usually results in errors 10 -7 per time step. However, these errors can accumulate over long duration simulations with ∼ 10 5 time steps, which could explain this minor discrepancy.

Acknowledgements This work has been financially supported by STFC. CPU time was provided by the UKMHD supercomputing facility located in Warwick. BF thanks the Cambridge Newton Trust for financial support. CB acknowledges support from a Herchel Smith Postdoctoral Fellowship of the University of Cambridge.