GRAVITY WAVES ON HOT EXTRASOLAR PLANETS. I. PROPAGATION AND INTERACTION WITH THE BACKGROUND

Chris Watkins; J. Y-K. Cho

doi:10.1088/0004-637X/714/1/904

1. INTRODUCTION

A stably stratified atmosphere, characterized by a positive vertical entropy gradient, can support gravity waves, or g-modes. Gravity waves are oscillations that arise from the buoyancy of parcels in the fluid. Such waves are readily excited by flow over thermal and surface topography, convective and shear instabilities, and flow adjustment processes. They propagate through the atmosphere both horizontally and vertically. Gravity waves in the terrestrial atmosphere and ocean are much studied (e.g., Gossard & Hooke 1975; Gill 1982). They have also been observed on other solar system bodies, such as Jupiter (Young et al. 1997) and Venus (Apt et al. 1980).

In the terrestrial atmosphere, a typical gravity wave has an energy flux of approximately 10⁻³ to 10⁻¹ W m⁻². Despite being small compared to the total amount of absorbed solar flux (∼237 W m⁻²), gravity waves are responsible for significantly modifying—even dictating—large-scale flow and temperature structures. Several well-known examples of this are the quasi-biennial oscillation, reversal of mean meridional temperature gradient in the upper middle atmosphere, and generation of turbulence (e.g., Andrews et al. 1987). We expect similar effects to be present in the atmospheres (and oceans) of extrasolar planets. Moreover, due to the greater irradiation and scale heights on them, the acceleration and heating effects of gravity waves can be much stronger on hot extrasolar planets.

There has been much interest in modeling the atmospheric circulation of extrasolar planets (e.g., Joshi et al. 1997; Showman & Guillot 2002; Cho et al. 2003, 2008; Burkert et al. 2005; Cooper & Showman 2005; Dobbs-Dixon & Lin 2008; Koskinen et al. 2007; Langton & Laughlin 2007, 2008; Menou & Raucher 2009; Showman et al. 2008). Accurate simulation of atmospheric circulation is crucial for interpreting observations of extrasolar planets, as well as for improving theoretical understanding in general. For this, the role of eddies and waves in transferring momentum and heat needs to be addressed (Cho 2008). This has long been recognized in solar system planet studies (e.g., Lindzen 1990; Fritts & Alexander 2003).

The plan of the paper is as follows. In Section 2, we derive the governing equation appropriate for linear monochromatic gravity waves on hot extrasolar planets. We also discuss a simple parameterization of the key nonlinear process, saturation. In addition, we present solutions to the equation for simple isothermal atmospheres, with and without shear in the background mean flow. In Section 3, we extend the calculation to a physically more realistic situation, by using background flow and temperature profiles derived from a three-dimensional (3D) hot-Jupiter atmospheric circulation simulation. This is the first such calculation to have been performed for extrasolar planets. Through this, the significant effects of gravity waves on hot extrasolar planet atmospheric mean flows are demonstrated. In this section, we also discuss a way in which gravity waves can transport momentum and heat horizontally—e.g., from the dayside to nightside on tidally locked planets. In Section 4, we discuss the implications of our work for current extrasolar planet atmospheric modeling work. We conclude in Section 5.

2. LINEAR THEORY

2.1. Taylor–Goldstein Equation

The dynamics of a linear gravity wave is described by the Taylor–Goldstein equation (TGE). This equation is derived from the full, 3D hydrodynamics equations (Batchelor 1967). In this work, we restrict the description to two dimensions and neglect rotation:

$\begin{equation} \frac{\mathrm{D}\mathbf {u}}{\mathrm{D}t} = -\frac{1}{\rho }{{\bm \nabla }}p + \mathbf {g} \end{equation} \tag{ 1a }$

$\begin{equation} \frac{\mathrm{D}\rho }{\mathrm{D}t} = -\rho {{\bm \nabla }}\,{\bm \cdot}\, \mathbf {u} \end{equation} \tag{ 1b }$

$\begin{equation} \frac{\mathrm{D}\theta }{\mathrm{D}t} = \frac{\theta }{c_p T}\dot{Q}, \end{equation} \tag{ 1c }$

where D/Dt = ∂/∂t + u · ∇; u = (u, w) is the flow in the horizontal and vertical directions (x, z), respectively; ρ is the density; ∇ = (∂/∂x, ∂/∂z); p is the pressure; g = (0, − g) is the gravity; θ is the potential temperature; c_p is the specific heat at constant pressure; T is the sensible temperature; and $\dot{Q}$ is the net diabatic heating rate. Equations (1) are supplemented with the ideal gas law:

$\begin{equation} p = \rho R T, \end{equation} \tag{ 2 }$

where R is the specific gas constant. Note that

$\begin{equation} \theta \equiv T \left(\frac{p_{\scriptscriptstyle R}}{p}\right)^{\kappa }, \end{equation} \tag{ 3 }$

where $p_{\scriptscriptstyle R}$ is some reference pressure (here taken to be the pressure at the lower boundary of the model) and κ = R/c_p; θ is related to the entropy s by ds = c_p dln θ.

The neglect of the third dimension and rotation requires that we restrict our analysis to waves with horizontal scale L ≲ U/Ω, where U is the characteristic mean flow speed and Ω is the planetary rotation rate. This scale is adequate for all gravity waves, except for large-scale tides (which has been recently considered by Gu & Ogilvie 2009). As an example, U/Ω ∼ 10⁷ m for HD 209458 b, based on U in hot extrasolar planet simulations of Thrastarson & Cho (2010). The resulting value is approximately 1/10 of the planet's radius. In addition, g, R, and c_p are taken to be constant and $\dot{Q}$ is specified. These restrictions do not mitigate the basic application and implications presented in this work. However, for broader applications, relaxation of these and other restrictions will be considered in future work.

The variables in Equations (1) are all expanded as a small perturbation about a mean value, which is a function of height only:

$\begin{equation} \zeta (x,z,t)\ =\ \zeta _0(z) + \zeta _1(x,z,t). \end{equation} \tag{ 4 }$

For the thermodynamic variables, we require that ζ₁/ζ₀ ≪ 1; however, this is not required for the flow variables, u and w. The mean state is assumed to be in hydrostatic balance, dp₀/dz = −ρ₀ g, and contains only horizontal flow so that w₀ = 0. We also assume the anelastic approximation (e.g., Ogura & Phillips 1962), ∇ · (ρ₀u) = 0. This implies the following:

$\begin{equation} \frac{N^2 H^2_{\rho }}{c^2_s} \ll 1 \end{equation} \tag{ 5a }$

$\begin{equation} \frac{\gamma D}{H_p} \lesssim 1 \end{equation} \tag{ 5b }$

$\begin{equation} \left| \frac{u_0}{w_1} \right| \frac{D}{L} \lesssim 1, \end{equation} \tag{ 5c }$

where N(z) = [g(dln θ₀/dz)]^1/2 is the Brunt–Väisälä frequency; H_ρ(z) ≡ |ρ₀(dρ₀/dz)⁻¹| and H_p(z) ≡ |p₀(dp₀/dz)⁻¹| are the density and pressure scale heights, respectively; γ = c_p/c_v is the ratio of specific heats with c_v being the specific heat at constant volume; c_s = (γRT)^1/2 is the speed of sound; and D is the vertical scale of the motion.

With Equations (4) and (5), we obtain

$\begin{equation} \frac{\partial u_1}{\partial t} + u_0 \frac{\partial u_1}{\partial x} + w_1 \frac{\mathrm{{\it d}}u_0}{\mathrm{{\it d}}z} = -\frac{\partial \Phi _1}{\partial x} \end{equation} \tag{ 6a }$

$\begin{equation} \frac{\partial w_1}{\partial t} + u_0 \frac{\partial w_1}{\partial x} = -\frac{\partial \Phi _1}{\partial z} + g \Theta _1 \end{equation} \tag{ 6b }$

$\begin{equation} \rho _0 \frac{\partial u_1}{\partial x} + \rho _0\frac{\partial w_1}{\partial z} + w_1\frac{\mathrm{{\it d}}\rho _0}{\mathrm{{\it d}}z} = 0\end{equation} \tag{ 6c }$

$\begin{equation} \frac{\partial \Theta _1}{\partial t} + u_0\frac{\partial \Theta _1}{\partial x} + w_1\frac{\mathrm{{\it d}}\ln \theta _0}{\mathrm{{\it d}}z} = \frac{\dot{Q}}{c_p T_0}, \end{equation} \tag{ 6d }$

where Φ₁ = p₁/ρ₀ and Θ₁ = θ₁/θ₀. Since the coefficients in Equation (6) are independent of x and t, we assume perturbations of the form,

$\begin{equation} \zeta _1(x,z,t) = \tilde{\zeta }(z)\,\exp \lbrace i(kx - \omega t)\rbrace, \end{equation} \tag{ 7 }$

where it is understood that the real part is to be taken. This leads to the polarization equations

$\begin{equation} {-}ik (c-u_0) \tilde{u} + \frac{\mathrm{{\it d}}u_0}{\mathrm{{\it d}}z} \tilde{w} =\ -ik \tilde{\Phi } \end{equation} \tag{ 8a }$

$\begin{equation} {-}ik (c-u_0) \tilde{w} = {-}\frac{\mathrm{{\it d}}\tilde{\Phi }}{\mathrm{{\it d}}z} + g \tilde{\Theta } \end{equation} \tag{ 8b }$

$\begin{equation} {-}ik \tilde{u} = \frac{\mathrm{{\it d}}\tilde{w}}{\mathrm{{\it d}}z} - \frac{\tilde{w}}{H_{\rho }} \end{equation} \tag{ 8c }$

$\begin{equation} {-}ik (c-u_0) \tilde{\Theta } + \frac{N^2}{g} \tilde{w} = \tilde{F}, \end{equation} \tag{ 8d }$

where $F \equiv \dot{Q}/(c_pT_0)$ is the forcing, c = ω/k is the constant (possibly complex, see Section 2.3) horizontal phase speed, and (c − u₀) is the intrinsic phase speed. Now, transforming into a new variable so that the effect of decreasing density with height is compensated for,

$\begin{equation} \hat{w}(z) = \tilde{w}\, \exp \lbrace -\chi (z) \rbrace \end{equation} \tag{ 9a }$

$\begin{equation} \chi (z) = \int ^z_{z_b} \frac{\mathrm{{\it d}}\xi }{2H_{\rho }(\xi)}, \end{equation} \tag{ 9b }$

where z_b is z at the bottom, we obtain the TGE:

$\begin{equation} \frac{\mathrm{{\it d}}^2 \hat{w}}{\mathrm{{\it d}}z^2} + m^2 \hat{w}\ =\ \frac{\kappa \,\dot{Q}}{H_p \left(c-u_0\right)^2}\, e^{-\chi }. \end{equation} \tag{ 10 }$

Here, m = m(z) is the index of refraction and corresponds to the local vertical wavenumber. It is given by

$\begin{equation} m(z) = \!\left[\!\frac{N^2}{\left(c-u_0\right)^2} + \frac{u_0^{\prime \prime }}{\left(c-u_0\right)} + \frac{u_0^{\prime }}{H_{\rho }\left(c-u_0\right)} - \frac{1}{4H^2_{\rho }} - k^2 \!\right]^{1/2}\!. \end{equation} \tag{ 11 }$

In Equation (10) we have used H_p = RT₀(z)/g and the prime in Equation (11) indicates differentiation with respect to z.

The vertical structure of the perturbations, oscillating in z, signifies that we are dealing with internal waves. In Equation (11), the key terms contributing to m² are the first ("buoyancy") term and the last ("non-hydrostatic") term. Although the other three terms contribute, generally the buoyancy and non-hydrostatic terms control whether the wave propagates since the flow shear and curvature are small and the scale height is large in practice. For large wavelength waves, the waves are hydrostatic and the non-hydrostatic term is small; then, the buoyancy term dominates. In these cases, as long as the atmosphere is stratified (i.e., N²>0) and c ≠ u₀, the wave will propagate vertically. However, for shorter, non-hydrostatic waves, it is possible that k²>[N/(c − u₀)]². In these cases, m is imaginary and the wave does not propagate. This situation is discussed in more detail in Section 3.2.

Equation (10) can be thought of as a driven harmonic oscillator equation. When m is real and constant, its solution is a simple sinusoid. Since the transformation, Equation (9), compensates for the exponential decay of density with height, w grows (decays) rapidly with height (depth). Here, we have dropped the tilde overscript. From hereon we drop all tilde overscripts and "1" subscripts for notational clarity. The growth can be clearly seen in Figure 1, along with the corresponding constant stress (vertical transport of horizontal momentum) when there is no dissipation (see Section 2.3).

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Gravity wave propagating in an isothermal (T₀ = 1350 K) and constant background flow (u₀ = 350 m s⁻¹) atmosphere. The phase speed of the wave c is 100 m s⁻¹ and the horizontal wavelength, 2π/k, is 2500 km. The vertical perturbation velocity w (⋅⋅⋅), horizontal perturbation velocity u (–), and wave stress τ (−·−) are shown. The latter is the vertical transport of the horizontal momentum. Wave amplitudes, u and w, grow exponentially with height, but the wave stress is constant with height, since there is no dissipation. The jump in τ near z/*H_p* = 1 is caused by the forcing.
Download figure:
Standard image High-resolution image

If the temperature or the flow varies with height, m is a function of height z. If m varies slowly, the WKB approximation (Bender & Orszag 1999) can then be used to obtain

$\begin{equation} w(z)\ \approx \ \frac{A\, e^{\, z/2H_{\rho }}}{m^{1/2}} \exp \left\lbrace \pm \, i\! \int _{z_{\rm b}}^z m(\xi)\, \mathrm{{\it d}}\xi \right\rbrace. \end{equation} \tag{ 12 }$

Here, A = w(z_b) [m(z_b)]^1/2. As in the constant m case, the vertical perturbation velocity is wave-like with upwardly and downwardly propagating components; the amplitude of the upward component grows with height and the downward component decays with depth. However, when m varies rapidly, the solution must be obtained numerically.

At the boundaries we use the radiation condition, selecting the upwardly propagating solution at the top boundary z_t and the downwardly propagating solution at the lower boundary z_b. This is achieved using the condition,

$\begin{equation} \frac{\mathrm{{\it d}}\hat{w}}{\mathrm{{\it d}}z} - \left(ism + {\displaystyle {1 \over 2m}} {\displaystyle {\mathrm{{\it d}}m \over \mathrm{{\it d}}z}} \right) \hat{w} \ =\ 0. \end{equation} \tag{ 13 }$

Here, s = ±1, depending on the signs of the horizontal and intrinsic phase speeds and on whether the condition is at the upper or lower boundary. Note that condition (13) requires the WKB approximation to be valid at the boundaries. For example, the boundaries cannot be critical layers, regions where (c − u₀) → 0, since m → ∞ approaching such a layer and the WKB approximation ceases to be valid. Critical layers will be discussed further below.

As alluded to above, we use Gaussian elimination (Canuto et al. 2006) in this work to numerically solve the TGE. In the cases presented here, 3000 equally spaced levels are employed to solve for w. We have checked that this resolution is adequate by performing calculations with 10,000 levels and verifying convergence at the higher resolution. Extensive validations against known analytic solutions, where they exist, have also been performed.

2.2. Polarization Relations and Fluxes

The solution to the TGE is a wave in the vertical velocity perturbation. This can be related to the horizontal and temperature (sensible and potential) perturbations. Understanding these is essential for parameterizing the saturation process in the full nonlinear situation and in general circulation models (GCMs). The perturbation quantities can be obtained from the polarization equations (Hines 1960),

$\begin{equation} u = \frac{i}{k} \left(\frac{\mathrm{{\it d}}w}{\mathrm{{\it d}}z} - \frac{w}{H_{\rho }} \right) \end{equation} \tag{ 14a }$

$\begin{equation} \Phi = \frac{i}{k}\left[ (c - u_0) \left(\frac{\mathrm{{\it d}}w}{\mathrm{{\it d}}z} -\frac{w}{H_{\rho }} \right) + \frac{\mathrm{{\it d}}u_0}{\mathrm{{\it d}}z} w \right] \end{equation} \tag{ 14b }$

$\begin{equation} \theta = {-}\frac{i}{k}\left[\frac{\theta _0 N^2}{g \left(c - u_0\right)}\right] w \end{equation} \tag{ 14c }$

$\begin{equation} T = -\frac{i}{k}\left[\frac{T_0 N^2}{g \left(c - u_0\right)}\right] w, \end{equation} \tag{ 14d }$

where we have introduced the geopotential perturbation function Φ (≡p/ρ₀). Equation (14a) gives the relationship between the vertical and horizontal velocity perturbations. Note that the amplitude of u is larger than the amplitude of w, as it is scaled by kH_ρ. Note also that the geopotential (pressure) perturbation varies with the background flow via the dependence on the intrinsic phase speed, according to Equation (14b). Equations (14c) and (14d) show that the potential and sensible temperature perturbations, respectively, are both π/2 out of phase with the vertical velocity perturbation. But, the phase between u and w varies locally through their dependence on the background.

Gravity waves are an efficient means of transporting both momentum and energy. The (perturbation) momentum and energy fluxes are simply obtained from Equations (14):

$\begin{equation} \tau = \rho _0\,\overline{uw} \end{equation} \tag{ 15a }$

$\begin{equation} F_x = \rho _0\,\overline{\Phi u} \end{equation} \tag{ 15b }$

$\begin{equation} F_z = \rho _0\,\overline{\Phi w}. \end{equation} \tag{ 15c }$

In Equations (15), τ is the vertical flux of horizontal momentum (or, the wave stress); F_x and F_z are the horizontal and vertical fluxes of energy, respectively; and the overline indicates an average over a wavelength,

$\begin{equation} \overline{\alpha \beta } = \case{1}{2}\,{\rm Re}\lbrace \alpha \beta ^* \rbrace, \end{equation} \tag{ 16 }$

where α and β are arbitrary complex functions and "*" denotes the complex conjugate. Note that the energy fluxes depend on the background flow through their dependence on ρ₀ and Φ. However, as can be seen in Figures 1 and 2, the wave stress remains constant, away from the forcing and damping regions—e.g., saturation regions and critical layers. This is in accordance with the second Eliassen–Palm theorem (Eliassen & Palm 1960), which expresses non-interaction of the disturbance in the absence of dissipation and forcing.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Same as in Figure 1, but with negative vertical wind shear: u₀ varies from 350 m s⁻¹ near the bottom to 650 m s⁻¹ near the top, with linear growth in between. Here, c = 500 m s⁻¹; hence, (c − u₀) → 0 at z/*H_p* ≈ 8, where the wave encounters a critical layer and dissipates. In the region just below the critical layer, the wave saturates and transmits momentum into the background flow, as indicated by the drop in τ with height there.
Download figure:
Standard image High-resolution image

In Figure 2, it is important to note that, where the fluxes are changing, the wave is interacting with the background flow in those regions. This should be contrasted with the behavior illustrated in Figure 1, where the stress is not changing. Changes in the wave stress cause accelerations to the mean flow. Correspondingly, changes in the energy fluxes cause the temperature of the region to change. The rates of these changes are given by

$\begin{equation} \frac{\partial u_0}{\partial t} = -\frac{1}{\rho _0} \frac{\partial \tau }{\partial z} \end{equation} \tag{ 17a }$

$\begin{equation} \frac{\partial T_0}{\partial t} = -\frac{1}{\rho _0 c_p} \frac{\partial F_z}{\partial z}. \end{equation} \tag{ 17b }$

In this case, the wave causes the background flow to accelerate from 350 m s⁻¹ toward 500 m s⁻¹.

In the remainder of this paper, when speaking of vertically propagating waves, we consider only waves that propagate upward from the region of excitation. However, it must be remembered that downward propagating waves are also generated. On a giant planet without a solid surface, those waves may not be reflected or absorbed. They can continue to penetrate downward until they encounter a critical level or a convective region. Or, they are simply dissipated since the amplitudes of the downwardly propagating waves decrease exponentially, as already discussed (and as also can be seen in Figure 2). The energy and momentum fluxes are linked by the first Eliassen–Palm theorem (Eliassen & Palm 1960):

$\begin{equation} F_z = \left(c - u_0 \right) \tau, \end{equation} \tag{ 18 }$

which can be derived from Equations (14) and (15). For downwardly propagating waves, we have F_z < 0; and, therefore, via Equations (18) and (17a), we see that for these waves the deposition of momentum also leads to acceleration of the flow toward the phase speed of the wave. Downwardly propagating planetary scale gravity waves (i.e., thermally excited tides) are considered by Gu & Ogilvie (2009).

2.3. Saturation and Critical Layers

The present work concerns inviscid, linear, monochromatic waves. Such waves are infinite in extent and, in principle, can grow without limit when they propagate upward. This is obviously not physical. In reality, such waves become unstable and saturate. The saturation process can be treated by introducing a correction to the solution according to the physical criterion,

$\begin{equation} \frac{\partial }{\partial z} \left(\theta _0 + \theta \right)\ \le \ 0, \end{equation} \tag{ 19 }$

ensuring that a wave does not become convectively unstable and remains at neutral stability. When the saturation region is identified, θ is adjusted so that neutral stability is maintained. This new value for θ is then used in Equation (14c) to obtain the corresponding w in the saturated wave. This value of w is then used in Equations (14a), (14b), (14c), (15), and (17). This saturation condition acquires the simple form,

$\begin{equation} \left| u \right| \ = \ \left| c - u_0 \right|, \end{equation} \tag{ 20 }$

in situations where the WKB approximation is valid (Fritts 1984).

As can be seen from the definitions of N, H_p, and θ, the Brunt–Väisälä frequency can be written as

$\begin{equation} N \left(z\right) = \left[ \frac{g}{T} \left(\frac{\partial T}{\partial z} + \frac{g}{c_p} \right) \right]^{1/2}. \end{equation} \tag{ 21 }$

Since g and c_p are essentially constants in the modeled height range, N depends only on the temperature profile T(z). For isothermal regions, N is a constant. In general, the fractional change of T with height is small compared to g/(Tc_p) throughout the modeled region. Hence, N is nearly constant in the entire domain, with a value that is roughly 2.4 × 10⁻³ s⁻¹. The Brunt–Väisälä frequency for our atmospheric profile is shown in Figure 3. Further discussion of the effects of variation in N is given in Section 3.1.2.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Sample atmospheric mean flow u₀ (–), temperature T₀ (−·−), and Brunt–Väisälä frequency N (⋅⋅⋅) profiles of a typical hot extrasolar planet HD 209458 b, used in this work. The profile is representative of a region at approximately 70°E, 10°N. The planet is a close-in giant planet. The profiles are obtained from a 3D, global circulation model, up to ∼10⁻³ bar level. Above that level, the profiles are simply extended, following loosely the observed profiles of Jupiter (Young et al. 2005; Flasar et al. 2004).
Download figure:
Standard image High-resolution image

If the background flow contains shear, it is possible for the wave to encounter a critical layer at some height. At the critical level, where c = u₀, the TGE becomes singular. However, the equation can be solved using the method of Frobenius (Bender & Orszag 1999), from which it is seen that the wave is drastically attenuated by the critical layer (Booker & Bretherton 1967). The amount of attenuation depends on the Richardson number Ri of the flow,

$\begin{equation} Ri\ =\ \frac{N^2}{\left(\mathrm{{\it d}}u_0/\mathrm{{\it d}}z\right)^2}\;. \end{equation} \tag{ 22 }$

Figure 2 illustrates a critical layer encountered by a wave. If the wave stress has magnitude τ below the critical layer, it then has magnitude, τ · exp{−2π[Ri − (1/4)]^1/2}, after the encounter with the critical layer (Booker & Bretherton 1967). This is a substantial amount of attenuation. For example, in our model atmosphere of Figure 2, Ri ≳ 900. Hence, the wave is essentially completely dissipated at the critical level, with the wave stress falling to practically zero and the momentum deposited into the mean flow. Note that, during its approach to the critical layer, the wave actually saturates and deposits momentum and energy over a finite layer (see Figure 2). It is important to note that while a wave saturates when approaching a critical layer, the presence of a critical layer is not required for saturation. Saturation is a general process referring to wave dissipation by many different mechanisms—e.g., radiation, conduction, breaking, and turbulence.

Numerically, when a critical layer is present, we lift the phase speed from the real axis by adding a small imaginary component: c = c_r + ic_i, where |c_i/c_r| < 10⁻³. This introduces a small amount of linear damping and ensures that the neglected nonlinear terms do not dominate in the regions where waves become steep and eventually break. Of course, adding damping causes the wave stress to decrease with height and the second Eliassen–Palm theorem to be no longer valid. However, the effect is small. This can be seen in Figure 2, where the wave stress falls negligibly over the layer from z/H_p ≈ 1 to z/H_p ≈ 5.

3. EXTRASOLAR PLANET APPLICATION

3.1. Setup

Our aim in this paper is to demonstrate several properties of gravity waves likely to be important for hot extrasolar planets. For this, we choose HD 209458 b as a paradigm planet. This planet has been the focus of many theoretical and observational studies, and it is expected to be generic with respect to the properties discussed here. The physical parameters that characterize the planet's atmosphere are given in Table 1.

Table 1. Parameters Used for HD 209458 b

Parameter	Symbol	Value
Specific gas constant	R	3523 J kg⁻¹ K⁻¹
Specific heat at constant pressure	c_p	12300 J kg⁻¹ K⁻¹
Acceleration due to gravity	g	10 m s⁻²
Rotation rate	Ω	2.08 × 10⁻⁵ s⁻¹

Download table as: ASCII Typeset image

3.1.1. Forcing

In a stratified atmosphere, gravity waves are readily generated by many mechanisms—thermal and mechanical—such as the absorption of stellar radiation, convective release of latent heat, storms, flows over topography and coherent localized "heat islands," and detonation by impacts. In this work, we consider small- and meso-scale thermally excited waves, rather than large-scale waves, as already discussed. The horizontal wavelengths used in this work are 2500 km or less. This is a reasonable range since it is well within the observed range of gravity waves on Jupiter (Young et al. 1997; Hickey et al. 2000). Although not unimportant, we do not dwell on the precise nature of the source of the excitations. The main interest here is in the propagation and deposition of momentum and energy.

The forcing in Equation (10) is simply represented as a Gaussian, modified so that it is zero beyond two half-widths above and below the center. The center is located at z/H_p = 1 above the bottom of the domain. The half-width is 75 km, or ∼0.15H_p. The forcing location and width are chosen so that the vertical scale of the forcing is less than the vertical wavelength of the waves we present here. We have extensively explored the effects of varying the location, width, and strength of forcing and present this case to illustrate several important points. Not surprisingly, the dynamics do depend on the chosen parameter values, but the dependence is broadly predictable. For example, if the forcing scale is much larger than the vertical wavelength, then only a very small amplitude wave is emitted from the forcing region, due to cancellations.

The heating rate, $\dot{Q} / c_p$ , is set to 10⁻³ K s⁻¹. Note that this is a modest value. The forcing corresponds to roughly 300 K per rotation of the planet. This is compared to ∼100 K per rotation at the chosen location in the circulation model of Thrastarson & Cho (2010). A forcing of ∼1100 K per rotation, for a similar latitude–longitude location on the planet (see Section 3.1.2), is used in Showman et al. (2008). The latter value implies that, in the absence of motion, the location on the planet will cool completely in one rotation of the planet. We stress that locally—i.e., scales far below the grid scale of the current circulation models—the forcing could actually be much stronger. The actual value is presently uncertain and likely to be spatially and temporally variable over the planet. To provide a context, for the Earth the heating rate is ∼2 K per day (1 day = 0.29 rotation of HD 209458 b) over large areas; but, locally, at the tops of low clouds on the Earth the rate can be up to ∼50 K per day (Wallace & Hobbs 2006).

3.1.2. Background Structure

Figure 3 shows the mean flow and temperature profiles used to obtain much of the results presented in this section. The lower part of both profiles—approximately the lower six scale heights—is taken from global circulation simulations of the hot extrasolar giant planet HD 209458 b by Thrastarson & Cho (2010), using the NCAR Community Atmosphere Model (Collins et al. 2004). The profile is from a point near the equator, slightly away from the substellar point (70°E, 10°N). This point was chosen as it is within the equatorial jet and the Coriolis parameter f = 2Ωsin ϕ is not large. The temperature profile generally increases with height over the lowest four scale heights and then becomes isothermal. Fortuitously, this provides an opportunity for a loose validation of our model: it is very similar to the temperature structure observed by the Galileo probe in the same region of Jupiter's atmosphere (Young et al. 2005). We extend the profile by keeping the temperature isothermal through the planet's stratosphere and having the thermosphere (beginning of the temperature inversion near the top) start between 12 and 13 scale heights at p₀ ≈ 4 × 10⁻⁶ bar. This profile has a Richardson number of at least 3.4, giving an attenuation of 1.4 × 10⁻⁵. Hence, any critical layer can be considered to fully dissipate the wave.

The chosen flow profile has two local flow velocity maxima. The upper maximum is extended into a jet with a peak at z/H_p ≈ 6. This is similar to the structure of the jet in Jupiter's stratosphere proposed by Flasar et al. (2004), which has a peak velocity between the 10⁻² and 10⁻³ bar levels. We then extend the profile further upward without shear. Note, this is in keeping with the lower boundary of the model of Koskinen et al. (2007). The structure is somewhat different than those of Showman et al. (2008, 2009), where there is just one jet with the peak approximately located at the 10⁻¹ bar level. The peak flow speed is also much greater in those studies at 4 or 5 km s⁻¹. It is important to note, however, that these differences do not change qualitatively the basic points we are making in this paper.

The Brunt–Väisälä frequency profile N(z) is also shown in Figure 3. As already discussed, the profile does not vary much over the whole domain: the maximum Brunt–Väisälä frequency is just 1.2 times the minimum value. Therefore, N does not contribute much to the variation of the index of diffraction m. The main contributor to the variation of m is the variation of the intrinsic phase speed, which is derived from the large variation of flow speeds. This should be compared to the analogous terrestrial situation, where the range of flow speeds is much lower. This allows N to have a larger effect on the variation of m on the Earth.

3.2. Wave–Background Interaction

3.2.1. Critical Layer Encounter

Figure 4 shows a gravity wave encountering a critical level in the upper jet in Figure 3. The wave has a horizontal wavelength, 2π/k, of 2500 km and c = 600 m s⁻¹. Here, since c>u₀ as the wave approaches the critical layer from below, the momentum deposited in the mean flow causes the mean flow to accelerate. This acceleration peaks at over 250 m s⁻¹ per rotation. This is large enough to double the flow speed at this layer in ∼2 rotations—a significant effect. The effect is large enough to require its inclusion in any simulation of the atmospheric circulation (Cho et al. 2008).

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Gravity wave with c = 600 m s⁻¹, propagating in an atmosphere with profiles shown in Figure 3. The horizontal perturbation velocity u (–), intrinsic phase speed c − u₀ (⋅⋅⋅), and the mean flow acceleration du₀/dt (−·−) are shown. The intrinsic phase speed becomes zero at z/*H_p* ≈ 5 and the wave encounters a critical level. In the layers just below the critical level, the wave saturates and sheds momentum into the mean flow, causing it to accelerate, peaking at a rate over 250 m s⁻¹ rotation⁻¹.
Download figure:
Standard image High-resolution image

The waves encountering critical layers are dissipated. Therefore, a flow with a range of flow speeds dissipates all gravity waves with phase speeds within this range. That is, a spectrum of gravity waves is prevented from propagating high into the atmosphere. There are other, secondary effects at the critical layer that may also affect the mean flow, but they are not modeled here. They will form the basis for future work. For example, the deposition of energy into the flow at the critical layer may well lead to the generation of new gravity waves, which then may propagate further, partly mitigating the filtering effect.

3.2.2. Saturation

Figure 5 shows an example of a gravity wave saturating in the atmosphere of Figure 3. In general, it is possible that a wave may not encounter a critical level as it propagates upward. However, such a wave then travels higher into the atmosphere, where it grows large and, if unabated, eventually breaks or suffers dissipation at higher altitudes. Although both momentum and energy are deposited in this case, we focus here on the effects on the temperature field.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** Gravity wave, with c = −40 m s⁻¹, propagating in an atmosphere described by the profiles in Figure 3. The perturbation to the temperature field T (–) and the heating rate dT₀/dt (−·−) are shown. The wave saturates at just above z/*H_p*= 10, where the heating peaks at nearly 80 K rotation⁻¹. The peak energy flux for this wave is approximately 1 W m⁻².
Download figure:
Standard image High-resolution image

The wave launched in Figure 5 has c = −40 m s⁻¹ (i.e., westward). The horizontal wavelength remains at 2500 km, as in the critical layer example of Figure 4. The vertical velocity perturbation grows with height. Therefore, so do the zonal velocity perturbation u and the potential temperature perturbation θ, as expected from Equation (14). In this case, the wave saturates near the top of the jet, where condition (19) is exactly satisfied. Here, condition (20) is approximately satisfied since the zonal perturbation velocity and the intrinsic phase speed (c − u₀) are both approximately 122 m s⁻¹. The saturation deposits energy that causes the atmosphere there to heat up. The heating is significant, peaking at ∼75 K per rotation—a 5% change in one rotation. In the absence of other effects, the ambient temperature can be doubled in about 20 planetary rotations (or orbits, assuming 1:1 spin–orbit synchronization).

A wave that has a phase speed greater than the maximum flow speed will not encounter a critical layer. Those with phase speeds close to, but still above, the maximum flow speed will, in general, saturate in the regions just below the maximum flow, as the intrinsic flow speed will be small in that region. Similarly, waves with phase speeds just less than the minimum flow will saturate as well. In this way, the filtering effect discussed above is extended beyond those waves with phase speeds equal to flow speeds. The main effect of these filtered waves on the flow will be lower in the atmosphere. In the profile given in Figure 3, where the waves emanate from the z/H_p = 1 level, this means that the upper layers of the lower jet will be slowed by gravity waves whereas the lower levels of the upper jet will be accelerated.

Those waves that do not dissipate will be able to propagate into the upper atmosphere depositing their momentum and heat there. Here, the changes to the flow can be very large. For example, the wave shown in Figure 5 causes a deceleration of up to 6.8 km s⁻¹ rotation⁻¹ as it saturates. This clearly dominates the flow at this level.

Moving the location of the wave origin down does not change the basic behavior in the qualitative sense. However, the amplitudes are much larger, compared with the case when the wave originates higher up in altitude. Thus, the possibility exists for stronger effects due to gravity wave interaction with the background.

This is illustrated in Figures 6 and 7, which should be compared with Figures 3 and 5. Here, we have extended the profiles downward. The wave is still launched at z/H_p = 1, but this is now deeper in the atmosphere. The wave again has a phase speed of −40 m s⁻¹, and the horizontal wavelength remains at 2500 km. This wave also saturates near the top of the upper jet, where the energy deposition into the mean flow causes the atmosphere there to heat up. Note that the region of heating is lower than when the wave originates higher up, as in Figure 7. The heating is significant, peaking at ∼3000 K rotation⁻¹. The ambient temperature can be doubled in approximately half of a planetary rotation. In a more realistic scenario, dissipation—which we have not included in our model—will likely reduce the heating rate.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** Same as in Figure 3, but with the bottom of the computational domain extended down to 100 bars. Here, u₀ is extended downward barotropically (independent of height) from the 1 bar level; T₀ is extended downward so that the profile below the 1 bar level is similar to that of Figure 18 in Showman et al. (2009).
Download figure:
Standard image High-resolution image

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** Gravity wave, with c = −40 m s⁻¹, propagating in an atmosphere described by the profiles in Figure 6. The perturbations to the temperature field T (–) and the heating rate dT₀/dt (−·−) are shown. The wave saturates at just above z/*H_p* = 13, where the heating rate peaks at just under 3000 K per rotation. In terms of the pressure level, this location is actually lower than in the case shown in Figure 5, and the magnitude of the peak is nearly 50 times greater. The peak energy flux for this wave is nearly 200 W m⁻¹. Thus, having a source at a lower height can have a much stronger effect.
Download figure:
Standard image High-resolution image

3.2.3. Refraction

So far we have been focusing on the vertical transport of momentum and energy by gravity waves. However, the waves can also transport momentum and energy horizontally (see, Section 3.2.5). Substituting c = ω/k into the index of refraction, Equation (11), and rearranging gives the dispersion relation for gravity waves. Here, we consider the case where u₀ = 0. We will examine cases where u₀ ≠ 0 in the sections following this one.

When there is no background mean flow, the dispersion relation simplifies to

$\begin{equation} \omega = \pm \frac{Nk}{\big[k^2 + m^2 + 1\big/\big(4H_{\rho }^2\big)\big]^{1/2}}. \end{equation} \tag{ 23 }$

We can then use the definitions,

$\begin{equation} u_g = \frac{\partial \omega }{\partial k} \end{equation} \tag{ 24a }$

$\begin{equation} w_g = \frac{\partial \omega }{\partial m}, \end{equation} \tag{ 24b }$

to obtain the group velocities. They are

$\begin{equation} u_g = \pm \frac{N\big[m^2 + 1\big/\big(4H^2_{\rho }\big)\big]}{\big[k^2 + m^2 + 1\big/\big(4H^2_{\rho }\big)\big]^{3/2}} \end{equation} \tag{ 25a }$

$\begin{equation} w_g = \pm \frac{N k m}{\big[k^2 + m^2 + 1\big/\big(4H^2_{\rho }\big)\big]^{3/2}}. \end{equation} \tag{ 25b }$

Thus, for propagating waves (i.e., waves for which m is real), u_g ≠ 0. Therefore, gravity waves always propagate obliquely and cannot strictly propagate vertically when there is no background flow.

From Equations (25) we see that ϑ, the angle of propagation with respect to the horizontal, is given by

$\begin{equation} \tan \vartheta \ =\ \frac{k m}{m^2 + 1/\big(4H^2_{\rho }\big)}. \end{equation} \tag{ 26 }$

Since H_ρ is nearly constant, with a value just under 500 km, tan ϑ varies with 1/m for m ≳ 10⁻⁶ m⁻¹. This gives rise to refraction. As a wave propagates into a region of higher m, it bends to a more horizontal path. Note that, as we are here considering a region with no flow, increasing m is essentially equivalent to increasing N. As shown in Figure 3, N increases in the thermosphere and so the paths of waves in this region will bend toward the horizontal even though the flow is small.

3.2.4. Trapping

From Equation (12) we can see that in regions where m is imaginary, the wave is evanescent: its amplitude decays toward zero and therefore it does not propagate vertically. This can occur when N² < 0 (i.e., when the atmosphere is convectively unstable). But, it can also occur for non-hydrostatic waves when the buoyancy term becomes dominated by the non-hydrostatic term. In addition, when the index of refraction changes between layers, the wave is reflected at the boundary. The amount of reflection is given by the magnitude of the coefficient of reflection |r|, where

$\begin{equation} r\ =\ \frac{m_1 - m_2}{m_1 + m_2}. \end{equation} \tag{ 27 }$

In Equation (27), m₁ and m₂ are the indices of refraction in two adjacent layers. When m₂ is imaginary, total reflection occurs and the wave is evanescent in that region and its amplitude decays to zero. However a region of propagation can exist between two regions of evanescence. This readily occurs for jets, where the intrinsic phase speed varies enough to allow the hydrostatic term to dominate in some regions and not in others. The region can also occur through variations of the Brunt–Väisälä frequency. In Figure 8, we see a wave that is trapped in the relatively quiescent region between z/H_p ≈ 1 and z/H_p ≈ 4. Outside this region, the value of Re(m) is small, or zero. Trapped in this region the wave is able to interact with itself and, under appropriate conditions, resonate.

Figure 8. Refer to the following caption and surrounding text. — **Figure 8.** Gravity wave, with c = 10 m s⁻¹ and horizontal wavelength 2π/k = 955 km, trapped in an atmosphere with the structure presented in Figure 3. The vertical velocity perturbation w (–), mean flow u₀ (−−), and the real part of the index of refraction m (−·−) are shown. The wave is trapped in relatively quiescent region between z/*H_p* ≈ 1.5 and z/*H_p* ≈ 3.5 and does not propagate vertically. The region of trapping corresponds to the region where m is real. The wave is reflected at the boundaries of this region, providing a possibility for resonance.
Download figure:
Standard image High-resolution image

This is another mechanism via which waves may be filtered out by the flow; and so not reach high altitudes at which saturation can occur. However, in this case, a trapped wave does not directly interact with a low level flow that changes its characteristics. Indeed, between the two reflection layers the wave can propagate horizontally—even in the absence of any flow, using the refractive mechanism described above. As long as the layers do not allow much leakage, it is possible for a trapped wave to cover large horizontal distances—transporting momentum and heat zonally (east–west direction).

3.2.5. Ducting

As well as being trapped in relatively quiescent regions, it is possible for waves to be trapped in a jet. As alluded to above, it is possible for such a wave to travel within the region of trapping, which is known as a duct or a waveguide. In Figure 9, a non-hydrostatic wave with c = 700 m s⁻¹ is trapped within the jet (located at ∼5 mbar level) in our model atmosphere. Note the small values of Re(m) outside the jet.

Figure 9. Refer to the following caption and surrounding text. — **Figure 9.** Same as in Figure 8, but with c = 700 m s⁻¹ and horizontal wavelength 2π/k = 1410 km. The horizontal group speed *u_g* (−−), vertical group speed *w_g* (⋅⋅⋅), and the real part of the index of refraction m (−·−) are shown. The wave, not shown, is trapped in the upper jet between z/*H_p* ≈ 3 and z/*H_p* ≈ 9, the region where m is real, and does not propagate vertically above this region; it is however, able to propagate along the jet as the large value of *u_g* within the jet shows. The wave is reflected at the boundaries of this region, providing a possibility for resonance.
Download figure:
Standard image High-resolution image

In this case, the flows are significant. Therefore, we use the full dispersion relation,

$\begin{eqnarray} \omega &=&\big\{ku_0 + 2H_{\rho }k[u_0^{\prime }+H_{\rho }(2(k^2+m^2)u_0+u_0^{\prime \prime })]\nonumber\\ &&\pm 2\big[H^2_{\rho }k^2\big(\big(1+4H^2_{\rho }(k^2+m^2)\big)N^2 +(u_0^{\prime }+H_{\rho }u_0^{\prime \prime })^2\big)\big]^{1/2}\big\}\big/\nonumber\\ &&\big\{1+4H^2_{\rho }(k^2+m^2)\big\}, \end{eqnarray} \tag{ 28 }$

to develop expressions for the group velocities.

However, the expressions obtained are large and rather unilluminating. They can be simplified by assuming that u'₀ and u''₀ are small. This is realistic since the shear is of the order of 10⁻³ s⁻¹ and u''₀ of the order of 10⁻⁸ m⁻¹ s⁻¹. This is small compared with the other terms in the expressions. This then gives

$\begin{equation} u_g = u_0+\frac{N\big[m^2 + 1\big/\big(4H^2_{\rho }\big)\big]}{\big[k^2 + m^2 + 1\big/\big(4H^2_{\rho }\big)\big]^{3/2}} \end{equation} \tag{ 29a }$

$\begin{equation} w_g = \pm \frac{N k m}{\big[k^2 + m^2 + 1\big/\big(4H^2_{\rho }\big)\big]^{3/2}}. \end{equation} \tag{ 29b }$

From these expressions, we can see that u_g follows u₀ as this is the larger term on the right-hand side of Equation (29a) in our model atmosphere. In Figure 9, the values of u_g and w_g are shown. Note that in the center of the duct w_g is very small while u_g is large, so that energy is transported along the flow. At the top and bottom of the duct the vertical group velocity increases, while the horizontal group velocity falls. Therefore, propagation here is nearly vertical. In Figure 9, we show w_g as positive, however this is only for the upward propagation of energy, at the top of the duct the wave is reflected and the vertical component becomes negative. This keeps the wave within the jet. The ray path followed by the wave, before the reflection, is shown in Figure 10.

Figure 10. Refer to the following caption and surrounding text. — **Figure 10.** Path of propagation of the wave in Figure 9 assuming that properties of the duct do not change in the x-direction. The path shown is the first crossing of the duct, that is until the wave encounters the top reflection layer. At this point, the wave will reflect and then propagate downward in a mirror image of this path. Note that the wave travels nearly one planetary radius before reflection. This means with just six reflections the wave will have nearly circumnavigated the planet. Of course in the real situation the properties of the duct will change in the x-direction and the wave will probably either leak out of the duct or dissipate before the circumnavigation is complete.
Download figure:
Standard image High-resolution image

The wave can travel large distances in this duct; but, eventually, the wave will either escape the duct or be dissipated. The range of speeds in the jet may change so that a critical level for the wave is created. The wave will then be reabsorbed into the flow. Alternatively, if the flow or Brunt–Väisälä frequency outside the jet changes so that the buoyancy term is no longer small and is dominated by the non-hydrostatic term, then the reflection is no longer total and the wave will then leak out of the duct. This can, for example, happen when propagating into a colder region, assuming the lapse rate remains constant. As can be seen from Equation (21), a fall in temperature with constant lapse rate will cause an increase in the Brunt–Väisälä frequency and thus an increase in the buoyancy term. This may occur very far from the original region of wave excitation. Indeed, in the example given, it is possible to envisage jets ducting waves and so transporting energy from the dayside of a tidally locked planet to the colder nightside where the waves escape the jet and propagate away from the duct before dissipating.

4. IMPLICATIONS FOR CIRCULATION MODELS

The effects of gravity waves discussed in this work on the larger scale circulation must be parameterized in global models because the spatial resolution—both horizontal and vertical—required to model them is currently prohibitive. The waves important to the large-scale extrasolar planet atmospheric circulation have horizontal length scales ranging from approximately ∼10⁵ m to ∼10⁷ m and vertical wavelengths as small as ∼10⁴ m. Waves with periods of few hours can carry significant momentum and energy fluxes vertically, but the sources of these waves include processes that are not included or resolvable in current circulation models.

The difficulty with representing gravity waves in GCMs exists even for the GCMs of the Earth. For example, the parameterization for convection is not aimed at producing realistic gravity waves (Collins et al. 2004). However, not representing gravity waves can affect the accuracy of the GCMs. The lack of gravity wave drag can lead to the overestimation of wind speeds, resulting in faster and narrower jets than observed (McLandress 1998). Further, gravity waves introduce turbulence with subsequent mixing and thermal transport (Fritts & Alexander 2003). This leads to greater homogenization of the atmosphere with a reduction in, for example, temperature gradients. Gravity waves also interact with planetary waves, playing a role in important transient phenomena (such as sudden stratospheric warming). In the absence of gravity waves, these phenomena are not accurately modeled (Richter et al. 2010).

There are many parameterization schemes currently incorporated or proposed for general circulation modeling (McLandress 1998). In all of the schemes, the basic components are (1) specification of the characteristic of the waves at the source level, (2) wave propagation and evolution as a function of altitude, and (3) effects on and by the atmosphere. All of them are essentially linear and one dimensional, in that waves only propagate vertically and that only vertical variation in the background influences the propagation. As seen in this work, linear theory still requires information such as the wave's phase speed c and wavenumber $\textbf {\it k}$ , for example. A more complete theory would need spatial and temporal spectral information. Intermittency is another crucial feature that would need to be taken into account. The primary differences in various schemes pertain to the treatment of nonlinearity and specificity of wave dissipation mechanisms.

Currently, all global circulation models of hot Jupiters suggest the presence of a low number of zonal jets. However, all the models do not have the resolution required to adequately resolve gravity waves and are subject to all of the limitations described above. This issue has been previously raised by Cho et al. (2008), in which they advocate caution against quantitative interpretation of current model results. For example, without the inclusion of the wave effects discussed in this work, high jet speeds and precise eastward shift of putative "hot spots" can be questioned (e.g., Knutson et al. 2007; Langton & Laughlin 2007)

5. CONCLUSION

Gravity wave propagation and momentum and energy deposition are complicated by the environment in which the wave propagates. For example, spatial variability of the background winds causes the wave to be refracted, reflected, focused, and ducted. Additionally, temporal variability of the background winds causes the wave to alter its phase speed. Still further complications arise due to the wave's ability to generate turbulence, which can modify the source or serve as a secondary source, and the wave's interaction with the vortical (rotational) mode of the atmosphere. Many of these issues are as yet not well understood and are currently areas of active research.

In this work, we have emphasized only some of these issues. We have shown that gravity waves propagate and transport momentum and heat in the atmospheres of hot extrasolar planets and that the waves play an important role in the atmosphere. They modify the circulation through exerting accelerations on the mean flow whenever the wave encounters a critical level or saturates. They also transport heat to the upper stratosphere and thermosphere, causing significant heating in these regions. Moreover, through ducting, they also provide a mechanism for transporting heat from the dayside of tidally synchronized planets.

Before relying on GCMs for quantitative descriptions of hot extrasolar planet atmospheric circulations, further work needs to be performed to ensure that the effects of important sub-scale phenomena, such as the gravity waves discussed here, are accurately parameterized and included in the GCMs.

The authors thank Heidar Thrastarson for generously sharing information from his simulation work on extrasolar planet atmospheric circulation. We also thank Orkan Umurhan for helpful discussions and the anonymous referee for helpful suggestions. C.W. is supported by the Science and Technology Facilities Council (STFC). J.Y-K.C. is supported by NASA NNG04GN82G and STFC PP/E001858/1 grants.