PHOTOMETRIC AND SPECTRAL SIGNATURES OF THREE-DIMENSIONAL MODELS OF TRANSITING GIANT EXOPLANETS

Circa 2010 April, more than 70 exoplanets had been discovered transiting their primary stars.⁴ Transits are useful in the characterization of exoplanets because they break the mass / orbital-inclination degeneracy of radial-velocity measurements. Moreover, a planet's radius can be directly measured, yielding mass (M_p) and radius (R_p ) correlations with which theorists can extract useful information about bulk composition and structure and can attempt to fit radius evolution models (Guillot et al. 1996; Burrows et al. 2000; Brown et al. 2001; Hubbard et al. 2001; Baraffe et al. 2003; Chabrier et al. 2004; Charbonneau et al. 2007).

However, since atmospheric opacity is a function of wavelength, a planet's transit radius is also a function of wavelength (Seager & Sasselov 2000; Brown 2001; Spiegel et al. 2007; Fortney et al. 2003, 2010). The variation with wavelength in the measured radius (actually an impact parameter) provides an ersatz spectrum directly related to the planet's atmospheric composition near the terminator. This radius spectrum, unlike the planet's spectrum at secondary eclipse, is more dependent upon composition than upon the temperature profile, and so provides complementary information to a planet's own direct emissions. The latter can not only be measured during secondary eclipse, but also during the traverse by the planet of its orbit as it traces out its phase light curve.

The transit radius spectrum and the planet emission spectrum as a function of phase can, therefore, together help constrain planet properties (Fortney et al. 2006, 2010). However, the stellar irradiation of such an exoplanet severely breaks what would otherwise be quasi-spherical symmetry, producing large day–night contrasts in thermal structure, zonal flows and banding, and violent atmospheric dynamics (Showman & Guillot 2002; Guillot & Showman 2002; Cho et al. 2003, 2008; Burkert et al. 2005; Cooper & Showman 2005, 2006; Showman et al. 2008, 2009, 2010; Dobbs-Dixon & Lin 2008; Dobbs-Dixon et al. 2010; Rauscher et al. 2007, 2008; Menou & Rauscher 2009; Rauscher & Menou 2010). Importantly, such behavior may have photometric and spectral signatures. It is the exploration of such signatures that motivates this paper. The state-of-the-art in the three-dimensional (3D) modeling of exoplanet atmospheres is still evolving and has not yet reached a level of maturity where detailed predictions for each known exoplanet are robust. Hence, in this paper we focus on a few generic ideas and conjectures that emerge from our modeling efforts. The procedure we have pursued is the following: first, using a 3D general circulation model (GCM) we derive dynamic model atmospheres at various epochs after they have achieved a steady state. Second, we post-process these 3D model atmospheres with a spectral atmosphere code to obtain transit radius spectra at various epochs during ingress, egress, and total transit. We assume equilibrium molecular compositions and solar metallicity. Third, using the same 3D model atmosphere, we integrate over the "visible" disk and calculate phase light curves as a function of wavelength and for various standard photometric bands. The results are generic predictions which, though not expected to be quantitatively precise and constraining on any particular giant exoplanet, nevertheless contain qualitatively interesting features that should inform future measurements. For specificity, we focus on two models of HD 209458b (Charbonneau et al. 2000; Henry et al. 2000), one with and one without an "extra absorber" and a thermal inversion (Hubeny et al. 2003; Burrows et al. 2007; Fortney et al. 2008; Spiegel et al. 2009; Knutson et al. 2008), but suggest our qualitative results are generic beyond this planet and this modeling paradigm (see also Fortney et al. 2006, 2010).

In Section 2, we summarize our methodology and techniques. Then, in Section 3, we describe the 3D models and various of their salient characteristics. We go on in Section 4 to present our results for the wavelength-dependent transit radius and ingress–egress asymmetries, and in Section 5 we turn to a full discussion of the derived light curves. This section contains not only wavelength-dependent planet–star flux ratios as a function of phase, but the phase variation of several photometric planet band fluxes. In Section 6, we summarize our general conclusions.

2.1. 3D GCM Models

We use the same modified version of the University of Reading's Intermediate GCM (Hoskins & Simmons 1975) as presented in Menou & Rauscher (2009) and Rauscher & Menou (2010), where a detailed description of the numerical implementation can be found. We use a horizontal spectral resolution of T31, roughly equivalent to a 4° resolution on the sphere, and 45 vertical levels which are logarithmically spaced in pressure from 100 bar to 10 microbar. We employ a simplified Newtonian relaxation scheme for the radiative forcing, so that everywhere the atmosphere heats or cools toward a prescribed 3D equilibrium temperature profile, on some representative radiative timescale (Iro et al. 2005). The equilibrium temperature profiles are chosen to have a horizontal dependence such that on the dayside the fourth power of the equilibrium temperature goes linearly with the cosine of the angle from the substellar point and on the nightside it is constant. The amplitude of the temperature difference between the substellar point and the nightside is taken from one-dimensional (1D) radiative transfer models (Burrows et al. 2008; see Figure 1), as are the dependences with pressure and depth of the equilibrium temperatures needed in the Newtonian scheme.

Zoom In Zoom Out Reset image size

Our choice for the angular dependence of the dayside equilibrium temperatures is loosely motivated by the fact that the equilibrium black body temperature (T_eq) of an irradiated surface at a slant angle is determined by setting the incident flux on the angled surface equal to the emitted flux and by the σT⁴_eq dependence of the latter. We scale the entire vertical profile of internal equilibrium temperatures to which the Newtonian cooling scheme is dynamically driving the actual temperature by this same angular factor. In this way, the expectation that the annuli on the planet away from the substellar point will be heated less by the star than the gas at the substellar point is realized. For the nightside, which by definition is not irradiated by the star, we employ the simplifying assumption that in equilibrium without nightside dynamics the planet would resemble more an isolated brown dwarf with uniform emission and profiles. Note that these assumptions about the solid-angle distribution of the equilibrium temperature profiles still allow the 3D GCM dynamics to redistribute heat in all regions of the planet's atmosphere and that the actual temperatures perforce deviate from these equilibrium temperatures due to such advection, sometimes to a significant degree. However, clearly our ansatz could be improved by doing the 3D radiative transfer that the problem will eventually demand. For now, we believe our approach allows us to capture the essence, if not the detail, of the thermal character of the dynamical atmospheres we are studying.

Thus, the only difference between the inputs for our two 3D GCM models of HD 209458b is the different 1D profiles of the substellar equilibrium temperatures, where one was calculated including an extra absorber high in the atmosphere and the other was calculated without. The radiative timescales are constant on a given pressure level and they increase with pressure at depth. They are taken from the profile in Iro et al. (2005), except for the deepest model levels. Figure 1 provides a detailed profile of the radiative times adopted. All other planetary and numerical parameters match those described in Rauscher & Menou (2010). Each model was run for 500 planetary days (≡ orbits), after which a statistically steady state has been reached for all levels at pressures less than 1 bar. Levels deeper than this are not directly probed by the types of observations considered here. To summarize, the key differences between the present models and those in Rauscher & Menou (2010) are (1) the different relaxation temperature profiles used and (2) the greater vertical extent of pressures modeled, up to 10 microbar in the present study.

2.2. Method for Deriving Transit Curves and Their Wavelength Dependence

2.3. Spectral Post-Processing to Derive Wavelength-Dependent Light Curves

Although the relaxation temperature profiles are very different for our two atmospheric models with and without a high-altitude absorber (Figure 1), we nevertheless find that these models have rather similar flow patterns. This is perhaps even more striking since these circulation patterns are also qualitatively similar to that described in Rauscher & Menou (2010), which has a rather different temperature relaxation profile. This suggests that the circulation pattern studied here is more strongly governed by the choice of the T⁴-behaves-linearly-with-cos θ horizontal dependence than by the amplitude or the detailed vertical structure of the relaxation temperature profile adopted. This supports the notion that the main qualitative results shown in the present study may be extended to other close-in giant planets.

While details differ, key flow features are common to both model atmospheres: transonic wind speeds, with peak wind speeds of 10 and 20 km s⁻¹ high up in the models with and without an extra absorber, respectively; a super-rotating (eastward) equatorial jet that extends across many pressure scale heights (see Figure 2); a pronounced chevron-shaped shock-like feature at a longitude of ∼135° at intermediate pressure levels; and large stationary vortices in each hemisphere that extend almost from the equator to the poles and are centered at a longitude of ∼−120° (Rauscher & Menou 2010). High in the atmosphere, where heat advection is not dominant, the day–night temperature difference is greater in the model run with an extra absorber because of the stronger forcing profile adopted. In both models, the horizontal temperature field is close to equilibrium, with a hot dayside and a cold nightside at pressures less than 1 mbar. The hottest atmospheric region is advected eastward by ∼45° at the 100 mbar level and by ∼90° at the 1 bar level, while the temperature field is fairly well longitudinally homogenized at deeper pressure levels.

Zoom In Zoom Out Reset image size

The general trend in these models is to have the hottest atmospheric regions advected farther away from the substellar point at deeper pressure levels, with potentially observable consequences as discussed in Section 5.2. For the model with the extra absorber, however, there is an additional complication. Although both models show localized heating where a shock-like feature develops in the flow, this heating is significantly more pronounced in the model with absorber. We attribute this difference to faster winds driven by a stronger day–night forcing in the model with absorber. This, rather than direct heat advection from the dayside, determines the location of the hottest atmospheric region at pressure levels 15–45 mbar in the model with absorber. As a result, the peak temperature region at these levels does not simply transition away from the substellar point gradually with depth, but it is significantly affected by a systematic shift to ∼135°, where the shock-like feature is located.⁵

Figure 3 depicts temperature maps with and without dynamical redistribution and with and without an extra optical absorber. To emphasize the contrasts, the atmospheric scale heights are exaggerated by a factor of 10. The longitudinal flows distribute heat to the nightside, and thereby partially even out the scale height bulge differential between the dayside and the nightside. Additionally, Figure 4 depicts the distribution of methane (CH₄) at a pressure level of 5.7 mbar, if in chemical equilibrium, over the surface of our 3D GCM models of HD 209458b, with and without an extra high-altitude optical absorber. If not in methane, carbon would reside in carbon monoxide in these models, which predominates on the daysides (particularly, in this model set, for the model with the extra optical absorber). In the model without the extra optical absorber, we note the presence of the tongue of methane-depleted (CO-rich) material dragged from the dayside to the nightside in the equatorial belt (see also Cooper & Showman 2006). The degree of such advection is a function of vertical height/pressure level. In principle, this feature could be detected in a precision comparison of the ingress and egress transit spectra (Fortney et al. 2010). However, we find, as did Fortney et al. (2010), that this differential effect can be quite small, less than one percent of the already small variation with wavelength expected for the transit radius. Be that as it may, these two figures partially represent the model context in which we have performed our post-processing exercise.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The top panel of Figure 5 shows the normalized fractional area of the atmosphere occulting the star as a function of wavelength as the planet enters and leaves transit for the model without an extra optical absorber at altitude. The fractions are all normalized to the instantaneous average value in the optical wavelength range from 0.5 to 0.7 μm. Hence, the average value of this fraction in that wavelength range (optical) is zero. For this model, the base position in the atmosphere where the chord optical depth is of order unity for optical photons (which defines the transit radius) is deep inside the atmosphere (near 0.1 bar). The bottom panel of Figure 5 is the same as the top panel, but for the model with the extra absorber at altitude. The action of the extra absorber is clearly seen in the optical region, and the normalization is as performed for the top panel. However, for this model, the chord optical depth is of order unity for a base pressure near a millibar, much higher in the atmosphere than for the model without the extra absorber. This difference results in perceptible differences between the two classes of models.

Zoom In Zoom Out Reset image size

Water is ubiquitous in our planet model atmospheres and its features dominate these spectra. For our models, the egress and ingress are rather symmetrical, differing by no more than ∼0.8% at any wavelength, due predominantly to slight differences in the scale heights near the respective terminators. These small ingress/egress differences are most prominent in the water bands, but there are slight differential signatures in the methane and carbon monoxide bands at ∼1.6 μm and ∼3.3 μm and at ∼2.3 μm and ∼4.6 μm, respectively, as might be expected from the advection of gas across the terminator depicted in Figure 4.

Figure 6 compares the transit radii during full transit as a function of wavelength for both models depicted in Figure 5. The normalization for both models is to the measured radius in the optical (∼1.32 R_J; Knutson et al. 2007). We believe this to be the most physically and observationally sensible normalization. Note that when so normalized the two models, which differ only in whether there is an extra optical absorber at altitude, predict very different radii, though the relative differences when comparing radii at different wavelengths outside the optical are similar. Therefore, comparing the radius as inferred from optical observations to that inferred from infrared observations is diagnostic of the presence of an extra absorber. Unfortunately, the infrared transit radius data are not consistent with either the extra-absorber model or the no-extra-absorber model. A comparison with the Spitzer/IRAC data from Beaulieu et al. (2010) and the 24 μm MIPS data from Richardson et al. (2006) suggests that neither model fits both sets of data simultaneously, but a major discrepancy is also the wide spread between IRAC 3/4 and IRAC 1/2. The same discrepancy between model and data was noted by Fortney et al. (2010) for their models and is as yet unexplained.

Zoom In Zoom Out Reset image size

Our implementation of the extra absorber puts it throughout the upper atmosphere, at the lowest pressures. Given this, and as Figure 6 demonstrates, the corresponding model does not evince much contrast in and out of the Na-D doublet. This is in contradiction with the observations of Charbonneau et al. (2002) and may hint at an upper boundary to the extent of the optical absorber (if it exists) if the Na-D data are to be accommodated.

The planet fluxes themselves, and their phase variation, should speak volumes about the atmospheres of hot Jupiters, in particular their thermal and compositional profiles and the longitudinal dependence of their circulation regimes (Cowan & Agol 2008). Using the methodologies described in Section 2.3, we now present results for the planet–star flux ratios and photometric variations for our two models of HD 209458b as it traverses its orbit.

5.1. Planet–Star Flux Ratio as a Function of Wavelength and Phase

Figure 7 portrays the planet–star flux ratios versus wavelength, for the models with and without the extra optical absorber and as a function of phase in the planet's (HD 209458b's) orbit. Here, zero phase is at secondary eclipse. Superposed are various data sets and the corresponding "1D" spectral models using the method of Burrows et al. (2007, 2008). The slight differences between the zero-phase 3D calculations and the "1D" predictions are due to the different methods of determining the average effect of emission from the dayside hemisphere and to the dynamics inherent in 3D models. We note many things about these curves. First, the data points at secondary eclipse are much more consistent with the zero-phase prediction(s) for the model with an extra absorber. This was pointed out in Burrows et al. (2007) and Knutson et al. (2008). Second, the variation with phase from the dayside to the nightside is much larger for the model with an extra absorber and inversion. This is a generic feature of predictions with strong optical absorption at altitude on the dayside. Third, there are interesting differences in the phase dependence at different wavelengths. The variations longward of ∼15 μm are generally largest. The variation at ∼10 μm is muted. The water feature near ∼6.2 μm becomes relatively more prominent on the nightside, while modest on the dayside. The phase variation near ∼4 μm is slight. The phase variation from ∼8 to ∼10 μm is much larger for the hot dayside model with the extra absorber. In sum, there are important variations in infrared colors with orbital phase that may be diagnostic of the atmosphere models.

Zoom In Zoom Out Reset image size

The planet–star flux ratios depicted in Figure 7 are integrals over the surface of the planet of brightness maps such as are shown in Figure 8. The associated panels in Figure 8 are representative planet brightness maps in the IRAC 3 band, the I band, and the J band, the latter for the models with and without the extra absorber. All the maps are at full phase (secondary eclipse). In principle, these are what the planet HD 209458b would look like in "glasses" for the various bands, given the various relative color maps. The color mappings depict the magnitude difference between the brightest point on the planet and are different for each given rendering and band. Note that the hot spot in each band is shifted by a different angle relative to the substellar point, reflecting the different approximate "photospheric" pressures of the various bands and the variation in the degree and strength of zonal heat advection with depth (Showman et al. 2009). In particular, the IRAC band fluxes, formed as they are at altitude, are shifted the least, by no more than ∼20° but mostly close to ∼ 0°. These differential phase shifts are useful diagnostics of the zonal flow regimes and models (Cowan & Agol 2008).

Zoom In Zoom Out Reset image size

5.2. Phase Dependence of the Photometric Band Fluxes

This is particularly clear from Figure 9, which depicts for the two different atmosphere models and for our fiducial HD 209458b assumptions integrated relative light curves versus orbital phase in the R, I, J, H, K, IRAC 1, IRAC 2, IRAC 3, and IRAC 4 bands. These data are derived by integrating maps such as are shown in Figure 8. The different phase shifts of the hot spots suggested in Figure 8 are manifest in these plots. While the shifts in the IRAC bands are smallest (generally being formed at altitude at lower pressures), the shift in IRAC 1 being the largest among these (∼20° for the no-absorber model and ∼10° for the absorber model), for the model with no extra absorber the shifts in the I and J bands can be ∼40° and ∼45°, respectively. For the absorber model, the shift in the J band is comparable, while the shift in the I band is only ∼10°. The fact that many of the curves are not flat on the far right of these plots near 180° is another indication of the skews in the phase light curves and of the (differential) dynamical advection of heat by planetary zonal winds.

Zoom In Zoom Out Reset image size

Ignoring the IRAC bands, the contrasts are generally largest for bands at the shortest wavelengths. In Figure 9, the differences between the models with and without the high-altitude absorber are clear, with the brightness variations with phase being larger for the absorber model in the IRAC, R, and I bands. The two models are comparable in the near-infrared bands. The specific numbers derived are not as important as the general trends revealed by this representative modeling exercise. In principle, by measuring such phase curves in different bands, one could back out information on the thermal and wind profiles as a function of depth and longitude. This is in part because the different bands are formed and have approximate photospheric positions at different pressures in the atmosphere. Such remote sensing of a hot Jupiter atmosphere would enable a new level of scrutiny for exoplanets and may be possible using James Webb Space Telescope (JWST).

Using a 3D GCM with Newtonian cooling and day- and nightside equilibrium temperature profiles based on two 1D spectral models (with and without a thermal inversion and a hot upper atmosphere), we created dynamical models of the transiting giant exoplanet HD 209458b. We post-processed these 3D model atmospheres with a detailed opacity code to obtain transit radius spectra at various epochs during the primary transit. Then, using a spectral atmosphere code, we integrated over the face of the planet seen by an observer at various orbital phases and calculated light curves as a function of wavelength and for nine different photometric bands from R (∼0.6 μm) through IRAC 4 (∼8 μm). The products of this study are generic predictions for the character of the expected phase variations both of a giant planet's transit spectrum and of its light curves. Since our GCM employed Newtonian cooling (and not radiative transfer using opacities that corresponded with those used in the post-processing), the calculations are slightly inconsistent. Nevertheless, the results have interesting features that should inform future measurements. We found that the temporal variations in the derived integral quantities for this model suite are small (less than 1%) and that the ingress/egress contrasts due to zonal flows, while also small (<0.5–1.0%), are most manifest due to the scale height asymmetries near the terminators and in the water bands. While models we have generated with and without the methane and carbon monoxide bands indicate that ingress/egress differences in the CO and CH₄ bands might be diagnostic of both carbon chemistry and abundance asymmetries due to zonal flows near the terminators, the differential effects due to different carbon-species abundances at the different terminators are small (less than a few tenths of a percent). A similar conclusion was reached by Fortney et al. (2010), who nevertheless found a slightly larger quantitative effect.

Since the transit radius of HD 209458b has been measured in the optical, it should not be allowed to differ when comparing models with and without an upper-atmosphere absorber and the predicted radii at other wavelengths must be determined relative to it. When this is done, the predicted radii in the near- and mid-infrared are very different (by as much as 0.05 R_J), though our current models do a poor job of fitting all the transit radius data for HD 209458b simultaneously (e.g., Knutson et al. 2007; Beaulieu et al. 2010; Richardson et al. 2006). The reason for this discrepancy is currently unknown.

We determined that the angular phase shifts of the brightness peaks relative to the orbital ephemeris are functions of photometric band. The J, H, and K bands are shifted most, by as much as ∼45°, while the IRAC bands are shifted least. This is because different bands are spectrally formed at different pressure depths and the zonal winds that advect heat are depth dependent. The IRAC band photospheres are generally at altitude, while those in the near-infrared are at deeper pressures near and above ∼0.1 bar. Therefore, the question of the magnitude of the downwind shift in the planetary "hot spot" due to equatorial winds is nuanced, depending upon the waveband in which measurements are made. This wavelength dependence of the phase shift in the brightness light curve, and the associated day/night contrasts, can in principle be used to constrain the circulation regimes of irradiated giant planets and to probe different pressure levels of a hot Jupiter atmosphere. Though our calculations focused on models of HD 209458b, similar calculations for other transiting hot Jupiters in low-eccentricity orbits should yield transit spectra and light curves of a similar character. As the subject of comparative exoplanetology matures and JWST comes online, such observational manifestations of global circulation may well become possible for a wide range of irradiated planets. Our models have been constructed to help prepare the way, however imperfectly, for that era.

We acknowledge useful conversations with Ian Dobbs-Dixon, Brad Hansen, and Ivan Hubeny for his general support of the COOLTLUSTY code. A.B. was supported by NASA grant NNX07AG80G and under JPL/Spitzer Agreements 1328092, 1348668, and 1312647. E.R. was supported by a NASA Graduate Student Research Program Fellowship, contract NNX08AT35H; and K.M. was supported by the Spitzer Space Telescope Program under contract JPLCIT 1366188. The authors are pleased to note that part of this work was performed while in residence at the Kavli Institute for Theoretical Physics, funded by the NSF through grant No. PHY05-51164.