Three-dimensional radiative-hydrodynamical simulations of the highly irradiated short-period exoplanet HD 189733b (original) (raw)

Journal Article

,

1Department of Astronomy, Box 351580, University of Washington, Seattle, WA 98195, USA

2NASA Astrobiology Institute, Box 351580, University of Washington, Seattle, WA 98195, USA

Search for other works by this author on:

1Department of Astronomy, Box 351580, University of Washington, Seattle, WA 98195, USA

Search for other works by this author on:

Received:

05 November 2012

Revision received:

08 August 2013

Published:

06 September 2013

• PDF
• Split View
• • Article contents
  • Figures & tables
  • Video
  • Audio
  • Supplementary Data
• Cite

Cite

Ian Dobbs-Dixon, Eric Agol, Three-dimensional radiative-hydrodynamical simulations of the highly irradiated short-period exoplanet HD 189733b, Monthly Notices of the Royal Astronomical Society, Volume 435, Issue 4, 11 November 2013, Pages 3159–3168, https://doi.org/10.1093/mnras/stt1509
Close

Navbar Search Filter Mobile Enter search term Search

Abstract

We present a detailed three-dimensional radiative-hydrodynamical simulation of the well-known irradiated exoplanet HD 189733b. Our model solves the fully compressible Navier–Stokes equations coupled to wavelength-dependent radiative transfer throughout the entire planetary envelope. We provide detailed comparisons between the extensive observations of this system and predictions calculated directly from the numerical models. The atmospheric dynamics is characterized by supersonic winds that fairly efficiently advect energy from the dayside to the nightside. Super-rotating equatorial jets form for a wide range of pressures from 10−5 to ∼10 bars, while counter-rotating jets form at higher latitudes. The calculated transit spectrum agrees well with the data from the infrared to the ultraviolet (UV) including the strong Rayleigh scattering seen at short wavelength, though we slightly underpredict the observations at wavelengths shorter than ∼ 0.6 μm. Our predicted emission spectrum agrees remarkably well at 5.8 and 8 μm, but slightly overpredicts the emission at 3.6 and 4.5 μm when compared to the latest analysis by Knutson et al. Our simulated Infrared Array Camera (IRAC) phase curves agree fairly well with the amplitudes of variations, shape, and phases of minimum and maximum flux. However, we overpredict the peak amplitude at 3.6 and 4.5 μm, and slightly underpredict the location of the phase curve maximum and minimum. These simulations include, for the first time in a multi-dimensional simulation, a strong Rayleigh scattering component to the absorption opacity, necessary to explain observations in the optical and UV. The agreement between our models and observations suggests that including the effects of condensates in simulations as the dominant form of opacity will be very important in future models.

1 INTRODUCTION

Among the ∼630 known exoplanets, HD 189733b currently boasts the most numerous and comprehensive set of observations to date. A member of the well-known class of short-period, highly irradiated gaseous planets (hot Jupiters), HD 189733b's relative proximity to earth and the favourable star to planet radius ratio have made it the preferred target for many different groups studying radial velocity, primary and secondary transits, transit spectra, emission spectra, variability, and phase curves across a wide range of wavelengths. Given this wealth of knowledge, it is the best candidate for detailed comparisons between model predictions and observations. In this paper, we compare the results from our three-dimensional (3D) radiative-hydrodynamical simulations of HD 189733b with this wide array of observations.

Hot Jupiters have been found around approximately 10 per cent of solar-type stars. Their short orbital periods (∼1–4d) suggest that these planets are subject to strong tidal forces that synchronize and circularize their orbits. As a result one side perpetually faces its host star, receiving intense irradiation, while the other side receives no stellar irradiation. In the absence of any hydrodynamics, the dayside reaches temperatures of several thousand degrees, while the nightside temperature, maintained by the internal energy of the planet, would remain at several hundred degrees. However, this large temperature differential across the planet acts as an enormous driving force for the atmospheric gas generating supersonic winds that advect energy across the planet. The resulting temperature distribution depends crucially on the non-linear coupling between the radiative transfer of energy and hydrodynamic transfer of energy. The spatially varying temperature across the planet leads to varying chemical composition, and emission and absorption efficiencies and is thus directly tied to observable properties. HD 189733b is a fairly typical example of such a planet.

There are a wide range of numerical approaches that have been taken to address the coupled radiation hydrodynamics of irradiated gas giants. Though coupled, it is convenient to separate a given studies approach for the radiative transfer from that for the hydrodynamics. Though not an exhaustive list, previous approaches to radiative transfer include relaxation methods (i.e. Newtonian heating) (Showman & Guillot 2002; Cooper & Showman 2005, 2006; Langton & Laughlin 2007, 2008; Showman et al. 2008; Menou & Rauscher 2009; Perna, Menou & Rauscher 2010; Rauscher & Menou 2010), kinematic constraints designed to represent incident flux (Cho et al. 2003, 2008; Rauscher et al. 2008), 3D flux-limited diffusion (FLD) with (Dobbs-Dixon, Cumming & Lin 2010) and without (Burkert et al. 2005; Dobbs-Dixon & Lin 2008) a separate radiation component, dual grey 1D two-stream approximation (Heng, Frierson & Phillipps 2011; Rauscher & Menou 2012, 2013), 3D FLD coupled to wavelength-dependent stellar irradiation (Dobbs-Dixon, Agol & Burrows 2012) and 1D frequency-dependent radial radiative transfer (Showman et al. 2009). Previous methods used to solve the hydrodynamical portion include solving the equivalent barotropic equations (Cho et al. 2003, 2008; Rauscher et al. 2008; Rauscher & Menou 2012), the shallow water equations (Langton & Laughlin 2007, 2008), the primitive equations (Showman & Guillot 2002; Cooper & Showman 2005, 2006; Showman et al. 2008; Menou & Rauscher 2009; Showman et al. 2009; Perna et al. 2010; Rauscher & Menou 2010; Thrastarson & Cho 2010; Heng et al. 2011; Thrastarson & Cho 2011; Polichtchouk & Cho 2012; Rauscher & Menou 2013), Euler's equations (Burkert et al. 2005; Dobbs-Dixon & Lin 2008) and the Navier–Stokes equations (Dobbs-Dixon et al. 2010, 2012). There are a wide range of assumptions in all these approaches. Currently, the most sophisticated approach to hydrodynamics utilizes the Navier–Stokes equations, while that for radiative transfer are 1D frequency-dependent models. However, it should be noted that the inclusion of additional physics (varying chemistry, clouds, magnetic fields, etc.) may prove equally if not more important.

Given the wide range of approaches utilizing a wide range of numerical codes, HD 189733b provides an excellent target for not only testing how well these models perform but also exploring the physical processes included. Therefore, in this paper we present 3D models specifically tuned for HD 189733b utilizing a model coupling wavelength-dependent two-stream approximation, wavelength-dependent stellar irradiation and the fully compressible Navier–Stokes equations. Section 2 presents the details on our modelling methodology including both the dynamics and radiation in the simulations and calculating observable signatures utilizing the models results. Section 3 presents results illustrating both the temperature and dynamical structure across the planet and detailed comparisons to observations. Section 4 concludes with a discussion of our results.

2 MODELLING METHODOLOGY

In this paper we present results from a numerical model coupling the fully compressible 3D Navier–Stokes equations to a two-stream, frequency-dependent radiative transfer calculation. Though we solve the Navier–Stokes equations in a manner similar to previous papers (Dobbs-Dixon & Lin 2008; Dobbs-Dixon et al. 2010, 2012) the radiative transfer methodology is presented here for the first time. Observables (phase curves, spectra, etc.) are calculated utilizing the results of the coupled radiative-hydrodynamical models. We describe all three components (hydrodynamics, radiative transfer and the calculation of observables) in detail below.

2.1 Hydrodynamics

The hydrodynamical portion of the model solves the fully compressible 3D Navier–Stokes equations given by

\begin{eqnarray} {{\mathrm{\partial} {\boldsymbol u}}\over {\mathrm{\partial} t}} + \left({\boldsymbol u}\cdot \nabla \right) {\boldsymbol u} &=& - \frac{\nabla {P}}{\rho } + {\boldsymbol g} -2{\boldsymbol \Omega \times {\boldsymbol u}} \nonumber \\ &&-\, {\boldsymbol \Omega \times }\left({\boldsymbol \Omega \times {\boldsymbol r}}\right) + \nu \nabla ^2{\boldsymbol u} +\frac{\nu }{3}\nabla \left(\nabla \cdot {\boldsymbol u}\right) \end{eqnarray}

(1)

\begin{equation} {{\mathrm{\partial} \rho }\over {\mathrm{\partial} t}} + \nabla \cdot \left(\rho {\boldsymbol u}\right) = 0. \end{equation}

(2)

\begin{equation} \left[ {{\mathrm{\partial} \epsilon }\over {\mathrm{\partial} t}} + ({\boldsymbol u}\cdot \nabla ) \epsilon \right] = - P \, \nabla \cdot {\boldsymbol u} - \frac{1}{r^2}\frac{{\rm d}}{{\rm d}r}\left(r^2 F_R\right) + S_{\star } + D_\nu , \end{equation}

(3)

representing the momentum, continuity and thermal energy equations, respectively. The 3D velocity is given by |${\boldsymbol u}$|⁠, ρ and P are the gas density and pressure, and |${\boldsymbol \Omega}$| is the rotation frequency. The gravitational acceleration, |${\boldsymbol g}$|⁠, is taken to be constant and purely radial. ν = η/ρ is the constant kinematic viscosity. We have neglected the coefficient for the bulk viscosity. The thermal energy equation contains heating and cooling contributions from the radial gradient of the radiative flux |${\boldsymbol F_{R}}$|⁠, the incident stellar energy _S_⋆ and the viscous heating _D_ν. The formalism for |${\boldsymbol F_{R}}$| and _S_⋆ is detailed in Section 2.2 and that for _D_ν can be found in Kley & Hensler (1987).

We solve these equations in spherical coordinates with a resolution of (N r, N_ϕ, N_θ) = (100, 160, 64), where ϕ is the longitude and θ is the latitude. At the equator this corresponds to cells that are 2| .circ_{.}^{\circ}.circ|25 by 3°. Flow over the pole is accounted for via the method described in Dobbs-Dixon et al. (2012).

2.2 Radiative transfer

In highly irradiated planets, transfer of energy via radiation plays a very important role. For HD 189733b this radiation can be broken into two distinct contributions: the irradiation the planet receives directly its host star and the local radiation comprised of reradiated stellar energy and the flux from the planets interior. Though we include a full-wavelength-dependent treatment, the incoming stellar irradiation can be largely characterized as short wavelength (visible) and the reradiated characterized as long wavelength (infrared).

2.2.1 Two-stream approximation

To address radiative transfer of the reradiated energy we have developed a frequency-dependent, two-stream approximation (Mihalas 1978) for the radial radiative flux, separating it into multiple upward (_F_↑, b) and downward (_F_↓, b) propagating channels. The governing equations in each wavelength bin can be written as

\begin{equation} \frac{{\rm d}F_{\uparrow ,b}}{{\rm d}\tau _{b}} = F_{\uparrow ,b}- S_{b}, \end{equation}

(4)

\begin{equation} \frac{{\rm d}F_{\downarrow ,b}}{{\rm d}\tau _{b}} = -F_{\downarrow ,b}+ S_{b}. \end{equation}

(5)

We have subdivided the full spectrum into 30 wavelength bins (identical to those in Showman et al. 2009; Dobbs-Dixon et al. 2012). Assuming that the gas emits as a blackbody, the source function S b in equations (4) and (5) in a given bin can be written as

\begin{equation} S_b = \pi \int _{\nu _1\left[b\right]}^{\nu _2\left[b\right]} B_{\nu }\left(T,\nu \right)\text{d}\nu . \end{equation}

(6)

These equations require two boundary conditions. For these we choose _F_↓, b(r = R p) = 0 at the surface and |$F_{\uparrow ,b}\left(r=R_{\text{bot}}\right)= S_b\left(T_{\text{base}}\right)$| at the interior, where we fix the total upward flux

\begin{equation} \displaystyle \sum \limits _{b} S_b\left(T_{\text{bot}}-\frac{{\rm d}T}{{\rm d}r}\Delta r\right)\Delta r = F_{\text{int}} \end{equation}

(7)

to the value required to match the observed radius. Note that though we set the downward flux to be zero at r = R p in equation (5), stellar heating is accounted for directly in the thermal energy equation (equation 3). This term is discussed further in the next subsection. The optical depth in equations (4) and (5) are calculated utilizing the averaged opacity within that bin and a diffusivity factor,

\begin{equation} \tau _b\left(r\right) = -1.66\int ^{r}_{R_p}\rho \kappa _b \,\mathrm{d}l. \end{equation}

(8)

The diffusivity factor of 1.66 is an approximation that accounts for an exponential integral that arises when taking the first moment of the intensity to calculate the flux (Elsasser 1942). The opacity κ_b_ is discussed in more detail below.

We solve equations (4) and (5) in the standard way, utilizing an integrating factor |$\text{e}^{\tau }$|⁠. Together with our boundary conditions we find

\begin{equation} F_{\downarrow ,b}\left(\tau _b\right) = \int _0^{\tau _b}S_b\text{e}^{-\left(\tau _b-\tau _b^{\prime }\right)} \,\mathrm{d}\tau _b^{\prime }. \end{equation}

(9)

Note that if we wished to include the stellar heating directly into these equations there would be an additional exponentially attenuated term |$S_{b,\star }\text{e}^{-\tau _b}$| added to the above equation for _F_↓, b, with

\begin{equation} S_{b,\star } = \left(\frac{R_{\star }}{a}\right)^2\pi \int _{\nu _1\left[b\right]}^{\nu _2\left[b\right]} B_{\nu }\left(T_{\star },\nu \right)\,\mathrm{d}\nu . \end{equation}

(10)

The equation for the upward flux does include an initial source term due to the flux from the interior and can be written as

\begin{equation} F_{\uparrow ,b}\left(\tau _b\right) = S_{b}\left(T_{bot}\right)\text{e}^{-\left(\tau _{b,bot}-\tau _b\right)} + \int _{\tau _{b,bot}}^{\tau _b}S_b\text{e}^{-\left(\tau _b^{\prime }-\tau _{b}\right)} \,\mathrm{d}\tau _b^{\prime }. \end{equation}

(11)

Finally, the net radial flux (defined with positive indicating outward flux) can be computed as

\begin{equation} F_{r}\left(\tau _b\right) = \displaystyle \sum \limits _{b} \left(F_{\uparrow ,b}\left(\tau _b\right)-F_{\downarrow ,b}\left(\tau _b\right)\right). \end{equation}

(12)

However, it is important to note that equation (12) only works in the upper atmosphere. As demonstrated by Rauscher & Menou (2012), in regions below the photosphere the temperature–pressure profile becomes isothermal. As the density increases, the integrated optical depth across a grid cell becomes very large, and the profile becomes dependent on the numerical resolution. The resolution required to compute an accurate temperature profile at depth using equation (12) far exceeds computational capabilities. Results from other groups that utilize some form of two-stream approximation (Showman et al. 2009; Heng et al. 2011) seem to contain similar isothermal regions at depth, though they do not discuss this issue. This effect is purely numerical and separate from the physically correct isothermal plateau region (e.g. Hubeny, Burrows & Sudarsky 2003) associated with the lower pressure of the infrared (IR) photosphere relative to the optical photosphere (note this is an expression of Sandstrom's theorem, which requires the heating be negatively correlated with height). To surmount the numerical issue, we follow Rauscher & Menou (2012) and transition to a diffusive scheme at τ_b_ = 2.5. This is done in each wavelength bin independently given that the photosphere can be at significantly different pressures across the spectrum, ranging from 10−1 to 5 bars. Given the larger optical depths, a diffusion approximation is entirely adequate in the deeper atmosphere. Below the active weather layer, radial radiative diffusion remains the dominant form of energy transport all the way down to the radiative-convective boundary. Finally, as seen in Dobbs-Dixon et al. (2010), non-radial radiative transfer does play a role in influencing the energetics. Though inconsequential near the sub-stellar point, the horizontal flux of radiation can play an important role. We have neglected the non-radial contribution here.

2.2.2 Stellar heating

The irradiation of the planets atmosphere from the central star is directly accounted for in equation (3) for the thermal energy of the gas. We choose this procedure rather than inserting an additional boundary term into equation (9) for _F_↓, b as it more accurately captures the slant-optical depth associated with the path of stellar photons through the atmosphere. The formalism here is identical to the procedure used in Dobbs-Dixon et al. (2012). The stellar heating term in the thermal energy equation can be written as

\begin{equation} S_{\star }= \left(\frac{R_{\star }}{a}\right)^2 \displaystyle \sum \limits _{b=1}^{nb}\frac{{\rm d}\tau _{b,\star }}{{\rm d}r} \text{e}^{-\tau _{b,\star }/\mu _{\star }}\int _{\nu _{b,1}}^{\nu _{b,2}}\pi B_{\nu }\left(T_{\star },\nu \right)\,\mathrm{d}\nu , \end{equation}

(13)

where _B_ν(_T_⋆, ν) is the stellar blackbody for HD 189733A taken from the Kurucz models.1_R_⋆ and a represent the stellar radius and semimajor axis. The slant-optical depth is accounted for by including μ⋆, the cosine of the angle between the normal and the incident stellar photons. The optical depth to stellar photons is calculated as

\begin{equation} \tau _{b,\star } = \int \rho \kappa _{b,\star }\,\mathrm{d}l. \end{equation}

(14)

2.2.3 Opacities

The details of the atmospheric opacity are a crucial parameter in determining the energy deposition and reradiation, subsequently playing an important role in the overall atmospheric dynamics. Dobbs-Dixon & Lin (2008) explored the role of changing the opacity within the context of a grey FLD approximation. The results of this study indicated that as opacity of the atmosphere is decreased, stellar energy is deposited at larger pressures, where the dynamics is more effective at redistributing it to the nightside. In this paper we are utilizing a much more sophisticated radiative transfer routine, but the general principle still holds, increasing opacities causes the energy to be deposited much higher in the atmosphere.

Here we are utilizing the frequency-dependent opacities κν of Sharp & Burrows (2007) which assume solar composition gas. The dominant molecules relevant for HD 189733b in these calculations are water, methane and carbon monoxide. The frequency-dependent opacities are averaged within each wavelength bin utilizing a Planck mean,

\begin{equation} \kappa _b = \frac{\int _{\nu _1\left[b\right]}^{\nu _2\left[b\right]} \kappa _{\nu }\left(\rho ,T\right)B_{\nu }\left(T\right)\text{d}\nu }{\int _{\nu _1\left[b\right]}^{\nu _2\left[b\right]}B_{\nu }\left(T\right)\text{d}\nu }. \end{equation}

(15)

However, we find that the exact form of the averaging method does not have a significant effect on the radiative solution. For the reradiated component of radiation (equations 9 and 11) the local gas temperature is used in _B_ν(T). When computing the contribution from the incident stellar irradiation (equation 14) we replace this with the stellar blackbody, _B_ν(_T_⋆).

After significant experimentation, we found that radiative-hydrodynamical models utilizing simply the opacities directly from Sharp & Burrows (2007) did not agree well with the observations. Comparison to emission spectra (Deming et al. 2006; Knutson et al. 2007, 2009; Charbonneau et al. 2008; Knutson et al. 2012) indicated that these initial models were too cool in the upper atmosphere, and comparison to transit spectra (Pont et al. 2008, 2013; Sing et al. 2011) indicated that the existing Rayleigh contribution in the opacities was woefully inadequate to explain the strongly varying transit radius with decreasing wavelength. To address this issue we have supplemented the opacities with two addition components, a uniform grey component |$\kappa _{\text{grey}}$| at all wavelengths and an component that scales as λ−4, as would be expected from Rayleigh scattering. This extra opacity can be written as

\begin{equation} \kappa _{e}\left(\lambda \right) = \kappa _{\text{grey}} + \kappa _{RS}\left(\frac{\lambda }{0.9\,\mu {\rm m}}\right)^{-4}. \end{equation}

(16)

We have run a number of simulations varying the magnitude of these two components. Setting |$\kappa _{\text{grey}}=0.035\,\text{cm}^2/g$| and |$\kappa _{RS}=0.6\,\text{cm}^2/g$| provides the best fits to the observations. Note that the magnitude of the grey component is the same as that required in the 1D models in Knutson et al. (2012). This extra opacity is added self-consistently to all the opacities utilized in the code and in the post-processing described in the next sub-section. This is the first multi-dimensional radiative-hydrodynamical model to include such a strong Rayleigh scattering component. Note however that we are treating κ_e_ solely as an absorptive opacity, ignoring its scattering properties. This should be included in future work.

2.3 Calculating observables from 3D models

To determine the transmission spectrum from the hydrodynamical models, we calculate the wavelength-dependent absorption of stellar light traversing through the limb of the planet. This allows us to determine the effective radius of the planet and a fractional reduction of stellar flux _F_⋆. Neglecting limb-darkening of the star, this can be expressed as

\begin{equation} \left(\frac{F_{\text{in transit}}}{F_{\star }}\right)_{\lambda } = \frac{\int \left(1-\text{e}^{\tau \left(b,\alpha ,\lambda \right)}\right)b\, \mathrm{d}b\, \mathrm{d}\alpha }{\pi R_{\star }^{2}}, \end{equation}

(17)

where τ(b, α, λ) is the total optical depth along a given chord with impact parameter b and polar angle α, defined on the observed planetary disc during transit.

To calculate the emission spectrum and phase curves, we integrate inwards along rays parallel to the line of sight. The emerging blackbody flux from the planet at each point on the surface is given by

\begin{equation} B\left(x,y,\lambda \right) = \int _{0}^{\infty }\frac{2\pi hc^2/\lambda ^5}{\text{exp}\left(\frac{hc}{\lambda kT}\right)-1}\text{e}^{-\tau }\,\mathrm{d}\tau . \end{equation}

(18)

To calculate the apparent dayside luminosity, we integrate over the observed disc

\begin{equation} L_p\left(\lambda \right) = \int I_{p}\mathrm{d}A, \end{equation}

(19)

where d_A_ is taken over the observer's plane and the emerging intensity is given by

\begin{equation} I_{p}\left(x,y,\lambda \right) = \frac{\lambda ^2}{c} \int B_{\lambda }\left(x,y,\tau \right)\text{e}^{-\tau _{\lambda }}\,\mathrm{d}\tau _{\lambda }. \end{equation}

(20)

For both transmission and emission calculations, the density and temperature needed to calculate τ at each location are interpolated from the values in the 3D models.

3 RESULTS

In this section we present the results of our 3D radiative-hydrodynamical simulations of the irradiated planet HD 189733b. We first discuss the wind and temperature structures throughout the atmosphere and then compare synthetic observations calculated using these results to a wide variety of observations.

3.1 Temperature and wind structure

As discussed in Section 1 the synchronous rotation and intense irradiation of the dayside by the host star HD 189733A drives winds throughout the upper portions of HD 189733b. These winds advect energy throughout the atmosphere resulting in significantly different temperature distributions than expected from purely radiative calculations. As is to be expected, the details of the dynamics depend on the input parameters, the model assumptions and the initial conditions (Thrastarson & Cho 2010); Table 1 lists a number of these parameters. The value of viscosity (⁠|$\nu =10^7\,\text{cm}^2\,\text{s}^{-1}$|⁠) was chosen to be slightly larger than the numerical viscosity. Viscosities below ∼106 produced no changes in the behaviour of the simulations, suggesting numerical diffusivity begins to dominate at this level.

Table 1.

Physical parameters chosen for simulations of HD 189733b presented in this paper taken from Torres, Winn & Holman (2008).

Table 1.

Physical parameters chosen for simulations of HD 189733b presented in this paper taken from Torres, Winn & Holman (2008).

The radial distribution of the energy deposited in the atmosphere from its host star extends from the upper regions of the atmosphere (10−5 bars) down to ∼10 bars due to the complicated wavelength-dependent opacities. Coupled with radial mixing, this leads to varying temperature distributions at different pressure levels within the atmosphere. Fig. 1 illustrates the temperature and zonal (i.e. longitudinal) velocity at constant pressure surfaces ranging from 10−4 to 1 bars. There are a number of features of interest in these plots. The hottest point in the atmosphere is near the sub-stellar point just above ∼0.1 bars, corresponding to the location where the majority of the stellar energy is deposited. The coolest points are high in the atmosphere in the nightside gyres found at mid-latitudes. Zonal jets, super-rotating at the equator and counter-rotating at mid-latitudes are also quite evident. Jet velocities increase with decreasing atmospheric pressure, reaching 6 km s−1 at pressures of 10−4 bars. The advective signature of the jets above 10−1 bars is seen near the terminators (90° and 270°) as tongues of hotter regions stretching on to the nightside. Deeper in the atmosphere, where the radiative time-scale is longer, advection is more efficient, equalizing temperature and moving the hottest point significantly downwind, east of the sub-stellar point. Also evident in Fig. 1 at large pressures (∼1 bar) is a westward extending hotter region at the equator near the western terminator (270°).

Figure 1.

Temperature in kelvin (left) and zonal speed in km s−1 (right) at 10−4, 10−3, 10−2, 10−1 and 1 bars from top to bottom. The images are projected on to a cos( latitude) map to highlight the importance of the equatorial features in determining observable properties. The sub-stellar point is located at 0° latitude and 0° longitude. The equator runs horizontally through the centre of each plot.

Fig. 2 illustrates the behaviour of the temperature and the zonal velocity with pressure at both the equator and mid-latitudes. Several separate temperature peaks are seen in these plots. The temperature inversion seen just below ∼10−2 bars near the sub-solar point again corresponds with the maximum energy deposition associated with the stellar heating. The larger temperatures at longitudes between 90° and 100° and pressures between 1 and 10−1 bars correspond to the heating caused when the westward mid-latitude jet deviates towards the equator and interacts with the eastward equatorial jet. At mid-latitude the heating between 1 and 10−1 bars is generated by the compression and shock-heating associated with converging flows. Much of this behaviour can also be seen in the zonal velocity shown in Fig. 1. Deeper in the planet temperatures converge on to a single curve that will continue to steepen until becoming unstable at the radiative-convective boundary.

Figure 2.

Temperature versus pressure (left-hand panels) and zonal velocity versus pressure (right-hand panels) at the equator and mid-latitudes. The inset colour wheel indicates the longitude of a given profile. The star is located to the right of the inset colour wheel; thus, the black line denotes the sub-stellar longitude. The eastward (positive) direction is counterclockwise.

The zonal velocity is almost entirely eastward (super-rotating) at the equator. However, unlike results presented in other studies, the zonal jets shown in Fig. 1 are significantly zonally asymmetric. This is potentially due to the combination of increased radial velocities allowed by solving the full momentum equation (i.e. the Navier–Stokes equations) and the fact that our simulation is compressible, a feature not included in other models. The largest velocities are found high in the atmosphere on the nightside of the planet just past the eastern terminator. The highest velocities are associated with the very large azimuthal pressure gradient found near the terminators. These flows have very high Mach numbers, reaching up to 3.5 across the nightside where the flow has cooled radiatively but still retains high velocities. With such high velocities it is likely that shocks play an important role in limiting the velocity. Our code utilizes artificial viscosity (Vonneumann & Richtmyer 1950; Stone & Norman 1992) to address the numerical difficulties associated with shocks. This approach calls for adding a non-linear viscous pressure term to the momentum and energy equations to spread the effect of the shock over multiple grid cells. Though the actual physical extent of the shock is much smaller, this method still results in the correct difference in entropy across the shock region. Thus, though our numerical scheme can capture much of the large-scale behaviour of shocks, it is likely that there are additional processes occurring on size scales smaller than our grid that contribute to the overall drag and heating of the flow. These velocities likely represent upper limits to the actual flow speeds. At mid-latitude, flow is in both the eastward and westward directions at low pressure. As the pressure increases, circumplanetary westward jets develop.

The details of jet formation (i.e. where they form and their strength) warrant an entire study in its own right that we feel is beyond the scope of this paper. Showman & Polvani (2011) suggest the interaction between the mean flow and standing Rossby waves act to pump momentum from high latitudes towards equator. Vertical motions induce opposite accelerations by transferring material with eastward (westward) momentum downward into the quiescent layer at the equator (mid-latitudes). However, accelerations due to horizontal eddies appear to dominate over vertical eddies. The result of this process is counter-rotating (westward) jets at mid-latitudes and super-rotating equatorial (eastward) jets seen in the simulations. Our initial calculations of momentum flux via eddies seem to corroborate the general idea, but work remains to be done to identify the eddy excitation mechanism.

3.2 Comparison to observations

HD 189733b is currently the most observed extrasolar planet given its relative proximity and favourable planet–star ratio. Therefore, in this section we perform detailed comparisons between these observations and synthetic observations calculated utilizing our numerical model. The methods for calculating observable metrics from the simulations are presented in Section 2.3.

In Fig. 3, we present the calculated transit spectrum from the UV through the IR. The calculated transit spectrum is dominated by Rayleigh scattering below ∼ 1 μm with molecular features becoming important at longer wavelengths. However, the uniform component of κ_e_ that we have added to the opacity suppresses the strength of many of these features. The observational data shown in Fig. 3 are a reanalysis of data from a number of groups by Pont et al. (2013). Absorption depths reported by a number of groups exhibit significant scatter (Ehrenreich et al. 2007; Knutson et al. 2007, 2009; Beaulieu et al. 2008; Pont et al. 2008; Désert et al. 2009, 2011; Sing et al. 2009, 2011; Gibson et al. 2012). Désert et al. (2011) and Pont et al. (2013) suggest that this may be due to stellar variability. Pont et al. (2013) have recently reanalyzed much of these data utilizing a consistent treatment for stellar spot corrections, transit properties and stellar-limb darkening. Our model agrees fairly well with these re-analysed data. We do see deviations at very short wavelengths ( < 0.6 μm) in the Rayleigh scattering portion of the spectrum. This is likely due to either a slight temperature inversion in the uppermost portion of the atmosphere that our model is not capturing or a change in grain size with height, physics not included in our current model. The very uppermost radius of our model corresponds to an absorption depth of ∼2.45 per cent and pressures of several micro-bars. We attempt to account for any additional absorption by adding an isothermal region above this pressure with the temperature within each radial column taken from the simulated region. The limits on the radial extent are a numerical constraint associated with very low density in this region and it may play a role in the deviations we see at wavelengths < 0.6 μm.

Figure 3.

The transit spectrum for the simulated planet HD 198733b. Overplotted are data from Pont et al. (2013) as solid points, and the band-averaged mode values as red triangles. Bands are taken from table 6 of Pont et al. (2013). Both the data and model exhibit strong Rayleigh scattering at short wavelengths and muted molecular features at longer wavelengths.

The next diagnostic we examine is the emerging spectrum on the dayside of the planet. Fig. 4 illustrates the calculated emission spectrum during the secondary eclipse. Again, we find general agreement throughout the IR, slightly overpredicting the newest measurements at 3.6 and 4.5 μm from Knutson et al. (2012), but doing very well at 5.8 and 8 μm. We underpredict the flux from Knutson et al. (2009) MIPS measurement at 24 μm. We do not reproduce the detailed shape seen by Swain et al. (2009) at ∼ 2 μm. As our overall goal is to utilize observations to validate our radiative-hydrodynamical model we plot the actual emerging intensity as λ_L_λ in Fig. 5. From this figure it is clear that the vast majority of the radiative energetics is occurring in the region from ≈ 1 to 10 μm. Outside this window the emerging intensity is orders of magnitude lower. The dynamics (jet formation, advection, meridional circulation, etc.) will overwhelmingly be driven by radiative energy transfer in this window. Therefore, the agreement between our models and observations in the Infrared Array Camera (IRAC) bandpasses gives us confidence that we are capturing a majority of the radiative energetics, despite some disagreement at longer wavelengths. To obtain a better sense of the planetary emission in the IRAC bandpasses, we plot the spatially resolved emission from the planet in Fig. 6 during the centre of the secondary eclipse. These figures illustrate a number of features including the offset of the hotspot, the advection of energy by the super-rotating equatorial jet and the varying level of asymmetry in different bandpasses, indicative of the varying pressure levels they probe. The IRAC 4 image compares favourably to the 2D secondary eclipse map of Majeau, Agol & Cowan (2012) and de Wit et al. (2012) utilizing 8 μm data.

Figure 4.

The dayside emission spectrum for HD 198733b shown as the black curve. Overplotted are data from Deming et al. (2006), Swain et al. (2009), Knutson et al. (2007), Charbonneau et al. (2008), Knutson et al. (2009) and Knutson et al. (2012).

Figure 5.

Apparent dayside luminosity λ_L_λ during secondary transit clearly illustrating the wavelengths where a vast majority of the radiative energy that ultimately drives the dynamics is found.

Figure 6.

The resolved emerging intensity given by equation (20) with units of erg cm−1 str−1 cm−2 s−1, during secondary eclipse as would it would be seen in the four IRAC bands across the entire face of the planet. The sub-stellar point lies in the centre of each image. The importance of the circumplanetary jet and the varying degrees of asymmetry are evident.

One of the most important diagnostics of atmospheric dynamics on extrasolar planets are multi-wavelength phase curves of the planet throughout its entire orbit. As varying phases of the planet come into view, we are able to directly compare the observed flux to that predicted by models. Deviations from radiative equilibrium indicate advective contributions to energy redistribution via the atmospheric dynamics. The presence of a strong circumplanetary jet was first robustly demonstrated for the planet HD 189733b (Knutson et al. 2007). Therefore, we examine the multi-wavelength phase curves. The amplitude, shape and the phases of the minimum and maximum of the flux are all important diagnostics. In Fig. 7, we show the simulated phase curves in the IRAC and MIPS bandpasses compared to observations from Knutson et al. (2009), Agol et al. (2009) and Knutson et al. (2012). In general we find remarkable agreement between our models and the data. We do very well in matching the observations of Agol et al. (2009) at 8 μm, including its distinct, non-symmetric shape. At 3.6 μm we slightly overpredict the flux at secondary eclipse (orbital phase of 0.5), and very slightly underpredict the flux from the nightside. At 4.5 μm we also overpredict the flux at secondary eclipse, but match the nightside flux fairly well. Finally, we significantly underpredict the flux at 24 μm. We saw this behaviour in Fig. 4 as well, but as discussed above the results at 24 μm likely have little influence on the overall atmospheric dynamics. However, in comparing our model to the phase curve observations it is important to note that Pont et al. (2013) claim that Knutson et al. (2012) overstate the ability of Gaussian Processing to correct for variability on the 1–2 d time-scale. Thus, the errors in their phasevcurve amplitude are likely ∼5 times too small. This suggests that our amplitudes may be consistent with the data given these larger uncertainties. Fig. 7 also shows the distinct offset of the phase of maximum and minimum fluxes from 0.5 and 0, respectively. Our models agree very well with the shape of the curve in the IRAC wavelengths of 3.6, 4.5 and 8 μm. The observations at 24 μm are too sparsely sampled to definitively compare the phase offsets. Table 2 lists the observed and simulated phase difference from transit and secondary eclipse for the minimum and maximum fluxes, respectively. Comparison between the observations and the models shows that we systematically underpredict the magnitude of the offset at all wavelengths, save an overprediction of the minimum at 4.5 μm. Taken at face-value this suggests that the jet strength may be somewhat underestimated in our models.

Figure 7.

The simulated IRAC phase curves for HD 189733b. Overplotted are observed phase curves from Knutson et al. (2009) at 24 μm, Agol et al. (2009) at 8 μm, and Knutson et al. (2012) at 4.5 and 3.6 μm.

Table 2.

Comparison between the calculated and observed offsets of the minimum and maximum fluxes (in hours) from the inferior and superior conjunction, respectively, in the IRAC and MIPS bandpasses. Observed values are taken from table 2 of Knutson et al. (2012). Note that the duration and accuracy of the 8 and 24 μm data was not sufficient define the minimum of the phase curve in these wavebands.

Table 2.

Comparison between the calculated and observed offsets of the minimum and maximum fluxes (in hours) from the inferior and superior conjunction, respectively, in the IRAC and MIPS bandpasses. Observed values are taken from table 2 of Knutson et al. (2012). Note that the duration and accuracy of the 8 and 24 μm data was not sufficient define the minimum of the phase curve in these wavebands.

Though phase curves provide an excellent diagnostic of the jet dynamics, Dobbs-Dixon et al. (2012) suggested a new observational diagnostic to determine the strength of the jet. In the presence of a strong jet, the temperature at the eastern terminator (in the direction of the jet) will be somewhat larger than that at the western terminator. The differential scale-heights between the two terminators will alter the timing of ingress and egress at a measurable level. This effect will change strength at different wavelengths as the day/night temperature differential depends on pressure and may be detectable with simultaneous observations at multiple wavelengths. In Fig. 8, we present the apparent offset in the time of central transit when the light curves computed from our model are fitted with a spherical planetary transit model (from Mandel & Agol 2002). We predict that the magnitude of this effect will be most visible in observations comparing the 3.6 and 4.5 μm regions of the spectra (with timing differences reaching up to 3 s), and also possibly at short wavelengths.

Figure 8.

The effective time-offset as a function of wavelength. Timing variations are primarily due to temperature differences between the eastern and western hemispheres.

The final diagnostic we explore is the stability of the atmospheric dynamics and thus the temperature distribution throughout the atmosphere. Agol et al. (2010) monitored secondary eclipse depths of seven events spread over 270 d. They find very little change in the observed depths, putting an upper limit on the variability of the dayside flux of 2.7 per cent. To test this, we have calculated the secondary eclipse depth variation over 30 orbits. Note that we do not actually have to simulate 30 orbits to derive this quantity. Because we do not have observational constraints for our models, we have simply continuously monitored the expected depth over this period. We find that the dayside emission from our models is extremely stable, with rms deviations of 0.1, 0.07, 0.07, 0.07 and 0.06 per cent at 3.6 4.5, 5.8, 8 and 24 μm, respectively. For comparison, as an unresolved source, Jupiter exhibits IR variability of the order of 20 per cent due to patchy clouds (Gelino & Marley 2000). However, a direct comparison between these planets may not be relevant, as Jupiter is largely convective, while HD 189733b has a radiative zone extending down to several thousand bars.

In summary, our model shows fairly good agreement with all the observations. It is possible that the jet is somewhat strongerthan we predict and thus more efficient at cooling the dayside and heating the nightside, as indicated primarily by the 3.6 μm phase curve. In addition, some fraction of the energy we see emitted at shorter wavelengths is actually escaping at longer wavelengths. Given the results presented in Fig. 5 only a very small fraction of the energy at short wavelengths must be transferred to longer wavelengths to explain the observations. There are a number of possible reasons that we are finding a discrepancy between the models and the observations. One likely culprit is the composition and thus the opacity of the atmospheric gas. As mentioned we are utilizing solar composition, molecular opacities from Sharp & Burrows (2007) augmented with an extra wavelength-dependent opacity of |$\kappa _{e}=0.035 + 0.6( \frac{\lambda }{0.9\,\mu \mathrm{m}})^{-4}$|⁠. In Fig. 9, we illustrate the effect of varying the strength of a simple grey component (κ_e_) on the phase curves. Indeed, it is clear that increasing κ_e_ causes the overall emission to shift from the 3.6 μm band to longer wavelengths. A small decrease at 3.6 μm is more than sufficient to cause large shifts at other wavelengths.

Figure 9.

The simulated phase curves for the four IRAC bands and the MIPS band for simulations with a varying strength of the extra grey absorber κ_e_.

4 DISCUSSION

We have presented a new set of radiative-hydrodynamical simulations for the irradiated gas giant planet HD 189733b. Simulations utilize a frequency-dependent two-stream approximation to calculate radiative transfer coupled to a 3D, fully compressible solution of the Navier–Stokes equations. Opacities, crucial for understanding the deposition of stellar energy and the subsequent reradiation and cooling of the atmosphere, are calculated assuming contributions from solar-composition molecules, a wavelength-independent grey component, and strong Rayleigh scattering. The latter two are assumed to come from the presence of grains in the atmosphere.

The atmospheric dynamics is characterized by supersonic winds that fairly efficiently advect energy from the dayside to the nightside. Super-rotating equatorial jets form for a wide range of pressures from 10−5 bars down to ∼10 bars. Counter-rotating jets form at higher latitudes flanking the equator setting up a pattern of three jets from north to south pole across the nightside. The westward, counter-rotating jets become slower at depth, but are also able to circumnavigate the planet below a pressure of ∼10−2 bars. Not only is the temperature structure modified by advection of energy in these jets, but the interaction between the jets also leads to compressional and shock heating, further modifying the temperature distribution across the planet.

We performed detailed comparisons between our model and a number of observations including transit spectra, emission spectra, eclipse maps, phase curves and variability. Model transit spectra agree fairly well with the data from the IR to the UV, though they underpredict the observations at wavelengths less than ∼ 0.6 μm. This may be due to an increase in temperature or a change in grain size, neither of which is captured by our models. Our predicted emission spectrum agrees remarkably well at 5.8 and 8 μm, but slightly overpredicts the emission at 3.6 and 4.5 μm when compared to the latest analysis by Knutson et al. (2012). Our measurement of the phase curve at 8 μm agrees very well, reproducing both the peak amplitude and the phase of maximum flux. As with the emission spectrum, we overpredict the peak amplitude at 3.6 and 4.5 μm, though agree fairly well with the phase offsets of minimum and maximum fluxes.

Among the other groups studying multi-dimensional radiation hydrodynamics of irradiated hot Jupiter atmospheres, Showman et al. (2009), Fortney et al. (2010) and Knutson et al. (2012) have also performed detailed comparisons to HD 189733b data utilizing the models of Showman et al. (2009). As both their model and opacities differ from ours it is not unexpected that we differ in some of the details. However, despite these differences we do appear to present a broadly consistent picture of the dynamics. Perhaps the most important difference between our results, with respect to the dynamics, is the location of the hotspot, indicative of the strength of the winds and the efficiency of advection. Showman et al. (2009) significantly overpredict this efficiency, finding the location of the hottest point somewhat further downwind than observed. As a result of increased advection, the nightside temperatures they predict are larger than observations. To address the discrepancy between the observed phase of the hotspot and that of the model, Showman et al. (2009) ran additional models with metallicity enhanced five times relative to solar (at odds with the measured sub-stellar metallicity; Torres et al. 2008) or a rotation period twice that of the orbital period (at odds with tidal synchronization theory; e.g. Rasio et al. 1996). Our fits do not rely on these mechanisms, but rather an additional opacity source most likely associated with condensates.

In an attempt to understand the discrepancy between their model predictions and the data at 4.5 μm, Knutson et al. (2012) suggest that advection of |$\text{CO}$| from the dayside to the nightside (where the carbon should convert to |$\text{CH}_4$| in chemical equilibrium) is important. Simulations of Cooper & Showman (2006) and Agundez et al. (2012), studying the carbon distribution, support this claim. Given the agreement between our simulations and the observations at 4.5 μm on the nightside, it appears that we do not require increased CO on the nightside. Moreover, the decreased advective efficiencies we find in our simulations would self-consistently predict less |$\text{CO}$| on the nightside as well. At the important shorter wavelengths of 3.6 and 4.5 μm, we overpredict dayside fluxes in comparison to observations of Knutson et al. (2012). If the amount of |$\text{CH}_4$| on the dayside was increased above chemical equilibrium values, as predicted by photochemical models (Visscher & Moses 2011), the flux at 3.6 μm would be suppressed bringing our models into better agreement. As noted in Knutson et al. (2012) this would also decrease the flux at 8 and 24 μm. However, it is important to remember that the addition of condensates may mask many details associated with non-equilibrium molecular distributions.

From our analysis of HD 189733b, it is clear that the hazes and/or grains play a significant, perhaps dominant, role in determining the overall opacity. Pont et al. (2013) suggest that molecular opacities may not be relevant for HD 189733b, and almost all the features of their data can be explained simply by a combination of a cloud deck and Rayleigh scattering dust grains. Thus, it is prudent to consider if the atmosphere can form grains and hazes. Ideally one would self-consistently calculate the formation, destruction, advection and settling of grains and hazes. This is not currently included in our model, but is an important direction for future research. The vertical velocities in the upper radiative zones are relatively modest compared to the zonal flows, but when coupled to the large scale-height of the atmosphere could lead to significant vertical mixing required to keep grains at low pressures.

We thank David Catling and Ty Robinson for discussions regarding the diffusivity factor, Heather Knutson for providing observational phase curves, and Adam Burrows for providing opacities. We would also like to thank the anonymous referee for comments that helped improve the paper. ID-D was partially supported by the Carl Sagan Postdoctoral program and the NASA Astrobiology Institute's Virtual Planetary Laboratory Lead Team, supported by NASA under solicitation NNH05ZDA001C. Addition support for this work was provided by NASA through an award issued by JPL/Caltech. We acknowledge support from NSF CAREER Grant AST-0645416. Finally, we would also like to acknowledge the use of NASA's High End Computing Program computer systems.

REFERENCES

,

Proc. IAU Symp. 253, Transiting Planets

,

2009

Dordrecht

Kluwer

pg.

209

,

ApJ

,

2010

, vol.

721

pg.

1861

,

A&A

,

2012

, vol.

548

pg.

A73

,

ApJ

,

2008

, vol.

677

pg.

1343

,

ApJ

,

2005

, vol.

618

pg.

512

,

ApJ

,

2008

, vol.

686

pg.

1341

,

ApJ

,

2003

, vol.

587

pg.

L117

,

ApJ

,

2008

, vol.

675

pg.

817

,

ApJ

,

2005

, vol.

629

pg.

L45

,

ApJ

,

2006

, vol.

649

pg.

1048

,

A&A

,

2012

, vol.

548

pg.

A128

,

ApJ

,

2006

, vol.

644

pg.

560

,

ApJ

,

2009

, vol.

699

pg.

478

et al. ,

A&A

,

2011

, vol.

526

pg.

A12

,

ApJ

,

2012

, vol.

751

pg.

87

,

ApJ

,

2010

, vol.

710

pg.

1395

,

ApJ

,

2008

, vol.

673

pg.

513

,

ApJ

,

2007

, vol.

668

pg.

L179

,

Heat Transfer by Infrared Radiation in the Atmosphere

,

1942

Milton, MA

Harvard University, Blue Hill Meteorological Observatory

,

ApJ

,

2010

, vol.

709

pg.

1396

,

ASP Conf. Ser. Vol. 212, From Giant Planets to Cool Stars

,

2000

San Francisco

Astron. Soc. Pac.

pg.

322

et al. ,

MNRAS

,

2012

, vol.

422

pg.

753

,

MNRAS

,

2011

, vol.

418

pg.

2669

,

ApJ

,

2003

, vol.

594

pg.

1011

,

A&A

,

1987

, vol.

172

pg.

124

et al. ,

Nat

,

2007

, vol.

447

pg.

183

et al. ,

ApJ

,

2009

, vol.

690

pg.

822

et al. ,

ApJ

,

2012

, vol.

754

pg.

22

,

ApJ

,

2007

, vol.

657

pg.

L113

,

ApJ

,

2008

, vol.

674

pg.

1106

,

ApJ

,

2012

, vol.

747

pg.

L20

,

ApJ

,

2002

, vol.

580

pg.

L171

,

ApJ

,

2009

, vol.

700

pg.

887

,

A Series of Books in Astronomy and Astrophysics: Stellar Atmospheres

,

1978

San Francisco

Freeman & Co.

,

ApJ

,

2010

, vol.

724

pg.

313

,

MNRAS

,

2012

, vol.

424

pg.

1307

,

MNRAS

,

2008

, vol.

385

pg.

109

,

MNRAS

,

2013

, vol.

432

pg.

2917

,

ApJ

,

1996

, vol.

470

pg.

1187

,

ApJ

,

2010

, vol.

714

pg.

1334

,

ApJ

,

2012

, vol.

750

pg.

96

,

ApJ

,

2013

, vol.

764

pg.

103

,

ApJ

,

2008

, vol.

681

pg.

1646

,

ApJS

,

2007

, vol.

168

pg.

140

,

ApJ

,

2008

, vol.

682

pg.

559

,

ApJ

,

2009

, vol.

699

pg.

564

,

A&A

,

2002

, vol.

385

pg.

166

,

ApJ

,

2011

, vol.

738

pg.

71

,

A&A

,

2009

, vol.

505

pg.

891

et al. ,

MNRAS

,

2011

, vol.

416

pg.

1443

,

ApJS

,

1992

, vol.

80

pg.

753

,

ApJ

,

2009

, vol.

690

pg.

L114

,

ApJ

,

2010

, vol.

716

pg.

144

,

ApJ

,

2011

, vol.

729

pg.

117

,

ApJ

,

2008

, vol.

677

pg.

1324

,

ApJ

,

2011

, vol.

738

pg.

72

,

J. Appl. Phys.

,

1950

, vol.

21

pg.

232

© 2013 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society

Citations

Views

Altmetric

Metrics

Total Views 813

485 Pageviews

328 PDF Downloads

Since 1/1/2017

Citations

99 Web of Science

×

Email alerts

Citing articles via

• Latest

• Most Read

• Most Cited

More from Oxford Academic