SPITZER SECONDARY ECLIPSE OBSERVATIONS OF FIVE COOL GAS GIANT PLANETS AND EMPIRICAL TRENDS IN COOL PLANET EMISSION SPECTRA

These observations were obtained in the 3.6 and 4.5 μm bandpasses using the Infra-Red Array Camera (IRAC) on the Spitzer Space Telescopes. HAT-P-19 was observed in full array mode with no peak-up for the 2011 4.5 μm eclipse, and in subarray mode with peak-up pointing (this pointing adjustment corrects the initial telescope pointing in order to place the star near the center of the pixel) for all subsequent visits. The 2012 3.6 μm eclipse was observed with a single peak-up at the start of the observation while the 2014 observations used the current standard approach, which includes both an initial peak-up of the position of the star on the array and a second peak-up 30 minutes into the observations in order to correct the initial pointing drift of the telescope. WASP-6 was observed in full array mode with no peak-up pointing, as this mode was not available when these data were taken. WASP-10 was observed in full array mode with no peak-up pointing in the initial 2010–2011 observations, while the 2014 observations were obtained in subarray mode with the standard dual peak-up pointing adjustments. WASP-39 was observed in subarray mode with no peak-up pointing, and WASP-67b was observed in subarray mode with the standard dual peak-up pointing method. Additional observation details can be found in Table 2.

For each data set, we calculate ${\mathrm{BJD}}_{\mathrm{UTC}}$ mid-exposure times using information in the BCD files provided by the Spitzer pipeline. We estimate and subtract the sky background, calculate the flux-weighted centroid position of the star on the array, and calculate the total flux in a circular aperture for each individual image as described in previous studies (Knutson et al. 2012; Lewis et al. 2013; Todorov et al. 2013; Kammer et al. 2014). We show the resulting light curves in Figures 1 and 2. We consider both fixed and time varying photometric aperture sizes in our fits. For the fixed apertures we consider radii between 2.0 and 5.0 pixels, where we step in increments of 0.1 pixels between 2.0 and 3.0 pixels and in 0.5 pixel increments for larger radii. The time varying apertures were calculated based on a scaling of the noise pixel parameter (Mighell 2005), which is proportional to the square of the full width half max of the stellar point spread function, and described by Equation (1) below:

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{(\sum _{n}{I}_{n})}^{2}}{\sum _{n}{I}_{n}^{2}}\end{eqnarray} \tag{ 1 }$$

where I_n is the measured intensity of the nth pixel. We iteratively re-scale the noise pixel aperture radii as $r=a\sqrt{\beta }+C$, where a is a scaling factor with a value of [0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 1.0, 1.05, 1.1, 1.15, 1.2], and C is a constant with a value between $-0.8$ and $+0.4$ stepping in 0.1 pixel increments.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We account for variations in intrapixel sensitivity using a pixel-level decorrelation (PLD) method, which Deming et al. (2015) found to produce results that are superior to those of a simple polynomial decorrelation or pixel mapping method. Details of this PLD technique are described in Deming et al. but the basic methodology is as follows. First, we obtain the raw flux values for a 3 × 3 grid of pixels centered on the position of the star, and then normalize these individual pixel values by dividing by the total flux in each 3 × 3 postage stamp. This normalization effectively removes any astrophysical signals, and we expect that the remaining variations in the fluxes of individual pixels as a function of time are primarily due to changes in the position of the star on the array. Equation (2) shows the resulting forward model, which we fit simultaneously with the eclipse model:

$$\begin{eqnarray}&&{F}_{\mathrm{model}}(t)=\displaystyle \frac{{\displaystyle \sum }_{i}{w}_{i}{F}_{i}(t)}{{\displaystyle \sum }_{i}{F}_{i}(t)}\end{eqnarray} \tag{ 2 }$$

where F_model is the predicted stellar flux in an individual image, F_i is the measured flux in the ith individual pixel, and w_i is the weight associated with that pixel. We leave these weights as free parameters in our fit, and solve for the values that best match our observed light curves.

We consider versions of the photometry with varying bin sizes, aperture sizes, and with up to two hours of data trimmed from the start of the time series in our fits. By fitting the binned light curves we are able to identify solutions with significantly less red noise on longer time scales and slightly higher scatter on the shortest time scales. The eclipses of the planets in this study typically have durations of several hours, and we are therefore most sensitive to noise on these time scales in our fits. For each fit we determine the best-fit model solution using the binned photometry, and then apply our best-fit model to the unbinned light curves in order to generate the plots shown here. We calculate our best-fit eclipse model and corresponding measurement uncertainties using the version of the photometry that minimizes the amount of time-correlated ("red") noise in the residuals after the best-fit model has been removed. As discussed in Deming et al. (2015), we evaluate the noise properties of a given set of residuals by calculating the root mean variance of the residuals as a function of bin size, where we step through bin sizes in intervals of powers of two. If the noise is Gaussian and there is no correlation between adjacent data points, we would expect the rms of the binned data to decrease as $1/\sqrt{n}$ where n is the number of points in each bin. We evaluate the relative amount of red noise in each version of the photometry by taking the difference between a model with the ideal $1/\sqrt{n}$ scaling and the observed rms at each bin size, squaring the difference, and summing over all bins. By minimizing this least-squares metric, we can select the version of the photometry which has the least amount of red noise across all time scales (see Figures 3 and 4).

Zoom In Zoom Out Reset image size

When considering the optimal choice of aperture and bin size, we first identify the version of the photometry with the lowest overall rms variance in the best-fit residuals and discard all of the alternative versions of the photometry with a rms variance that is a factor of 1.2 greater than this value. This ensures that we avoid solutions that minimize the red noise component by significantly increasing the amount of white noise (these are usually the largest apertures, where the sky background contributes additional noise). We then pick the bin size and photometric aperture that has the smallest amount of red noise as measured by our least squares metric.

We find that there is an exponential ramp visible at the start of some visits, which is a well-known feature of the IRAC 3.6 and 4.5 μm arrays and has a shape that varies depending on the illumination history of the array immediately prior to each observation as well as the brightness of the target star (e.g., Lewis et al. 2013; Zellem et al. 2014). We remove this ramp by trimming data from the start of each light curve. After identifying the optimal bin size and aperture choice, we determine the correct amount of data to trim by examining the best light curve after detector effects have been removed and trimming until no ramp is visible at the start of the observations. We then revisit our choice of bin size and aperture with the new trim duration in order to make sure that our previous choice is still the best one with the new trim interval. As shown in Figures 3 and 4, this results in light curves with no detectable red noise component on the time scales spanned by our observations. We also calculate a Lomb–Scargle periodogram for the best-fit residuals from our nominal light curve for each visit and find that for periods between five minutes and five hours there are no peaks in the power spectrum with false alarm probabilities lower than 10%. As a final step, we check to make sure that we obtain eclipse depths that are consistent at the 1σ level across different versions of the photometry that meet our 20% excess noise threshold and have minimal red noise components.

Zoom In Zoom Out Reset image size

We use the routines described in Mandel & Agol (2002) in order to generate our eclipse light curves, where we set the planet–star radius ratio, orbital inclination, and the ratio of the orbital semi-major axis to the stellar radius equal to their best-fit values from previous studies (see Table 1 for a complete list of references for each planet). We assume a circular orbit for the purpose of calculating the eclipse shape and duration, as there is no evidence for an orbital eccentricity large enough to appreciably alter the eclipse duration (e sin $\omega \gt 0.1$, where e is the orbital eccentricity and ω is the longitude of periastron) in any of these systems. We allow the individual eclipse depths and center of eclipse times to vary as free parameters in our fits. For the second 4.5 μm eclipse observations of HAT-P-19b and both 3.6 μm observations of WASP-67b the eclipse is detected with marginal significance and we therefore place a Gaussian prior on the phase of the secondary eclipse in order to reduce the corresponding uncertainties on the eclipse depth. We do not fix the center of eclipse time in these fits, but instead implement this prior as a penalty in ${\chi }^{2}$ proportional to the deviation from the error-weighted mean center-of-eclipse phase and corresponding uncertainty from the two 3.6 (HAT-P-19b) or 4.5 (WASP-67b) μm eclipses.

We find that the use of this prior changes the best-fit eclipse depth for the second 4.5 μm observation of HAT-P-19b by less than 0.1σ, while reducing the measurement error on this depth by 40% as compared to a fit with no constraint on the center of eclipse phase. For consistency, we apply the same prior to the first 4.5 μm eclipse observation of HAT-P-19b, which is detected with a significance of 3.2σ in fits with the standard uniform prior on the eclipse phase. We find that in this fit the use of the prior has no effect on the best-fit eclipse depth and reduces the measurement error in this parameter by 5%, as expected for a case where the data themselves provide good constraints on the eclipse time. For WASP-67b our 3.6 μm fits with a uniform prior prefer an eclipse solution at a much earlier phase with a best-fit eclipse depth that is consistent with zero at the 1σ level. We therefore conclude that this is likely a local minimum in ${\chi }^{2}$ caused by a trade-off with the instrumental noise model, and place a prior on the eclipse phase based on the measured 4.5 μm eclipse phase. This allows us to obtain new constraints on the eclipse depth at 3.6 μm for a range of phases consistent with the measured 4.5 μm eclipses.

We determine the best-fit model for each individual eclipse using a Levenberg–Marquardt minimization routine, where each fit has 13 free parameters. In addition to the nine pixel weight parameters and two eclipse parameters, we also fit for a linear function of time in order to account for long-term instrumental and stellar trends. Figures 1 and 2 show the raw photometry, while Figures 5 and 6 show the corresponding normalized photometry after dividing out the best-fit PLD model and linear function of time. Table 3 lists the best-fit eclipse depths and center of eclipse times for each observation. We estimate the corresponding uncertainties on our model parameters using a Markov chain Monte Carlo (MCMC) analysis with 10⁶ steps per chain, where we define the 1σ uncertainties as the symmetric interval around the median parameter value containing 68% of the total probability. We assume uniform priors on all parameters with the exception of the 3.6 μm eclipse phase of WASP-67b and the 4.5 μm eclipse phase for HAT-P-19b as noted in the previous paragraph.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We initialize our chains using the best-fit solution from the Levenberg–Marquardt minimization in order to minimize the burn-in time. We then check to see where the ${\chi }^{2}$ value of the chain first drops below the median value over the entire chain, and trim all points prior to this step in order to remove any residual burn-in phase. We also examine our resulting posterior probability distributions in order to make sure that our choice of prior bounds is broad enough that it does not affect our final results, and re-run our chain with a broader prior range if required. Although we find that individual pixel weights are often degenerate with each other, our probability distributions for the best-fit eclipse depths and center of eclipse times are Gaussian and uncorrelated with the other fit parameters, with the exception of the 3.6 μm eclipse of WASP-6b which we discuss in more detail below. The correlation between pixel weights could be mitigated by a transformation to an independent set of basis vectors (e.g., Morello et al. 2014, 2015), but as these correlations do not affect our final best-fit eclipse parameters and their corresponding uncertainties we do not pursue this option here.

We find that our measured center of eclipse time for WASP-6b's 3.6 μm eclipse has an uncertainty that is significantly larger (15 minutes) than is typical for our other observations with statistically significant eclipse detections. We examined the results from our MCMC analysis and identified a correlation between the slope of the linear trend across the 3.6 micron eclipse and the measured center of eclipse time in this visit (see Figure 7). As discussed earlier, we found it necessary to trim the first 1.5 hr of data in order to remove an unusually long-lived ramp at the start of these observations. This trimming means that we have a relatively short baseline prior to the start of the eclipse, and it is therefore not surprising that there might be a correlation between the eclipse time and the linear trend in our fits. Our best-fit eclipse depth for this visit does not appear to be correlated with either of these parameters, and we find no evidence for a similar correlation between eclipse time and linear trend slope in other visits with short pre-eclipse baselines.

Zoom In Zoom Out Reset image size

WASP-10b is the only planet in our study that was previously known to have an eccentric orbit, with $e=0.0473\pm 0.0032$ based on a total of sixteen radial velocity measurements from Keck HIRES and SOPHIE (Christian et al. 2009; Knutson et al. 2014b). Our error-weighted average secondary eclipse phase of 0.4649 ± 0.004 confirms this planet's non-zero orbital eccentricity, and we combine our measurement with published radial velocity data from Knutson et al. (2014b) in order to obtain improved constraints on the orbital eccentricity e and longitude of periastron ω (see Figure 8). We also add two new Keck HIRES measurements that extend the baseline of these data by approximately a year and a half beyond the last measurement reported in our previous paper. These new measurements clearly establish the presence of curvature in the observed radial velocity trend, and we have adapted our fits accordingly.

Zoom In Zoom Out Reset image size

We give the updated orbital parameters from a joint fit to the radial velocity measurements and using the measured transit ephemeris and secondary eclipse times in Table 4 (see Knutson et al. 2014b, for additional information on the fitting method). Our new orbital eccentricity of $e={0.0589}_{-0.0026}^{+0.0033}$ differs from our previous value by 4σ; this change can be attributed to the inclusion of our measured secondary eclipse times in the new fits as well as the switch to a curved radial velocity trend model. We note that a previous study reported a detection of the secondary eclipse for this planet in the Ks band at an orbital phase of 0.4972 ± 0.0005 (Cruz et al. 2015); this offset differs significantly from that of our Spitzer secondary eclipse observations and is also inconsistent with the current radial velocity measurements for this system. The authors of this study note that their reported error is likely an underestimate as it does not account for additional uncertainties contributed by the removal of systematic noise sources from their light curve. See Section 3.2.1 for additional discussion of this observation.

Published radial velocity observations for the remaining four planets in this study suggest that their orbits are likely circular, in good agreement with our measured secondary eclipse times. WASP-6b has the most extensive radial velocity data set, with a total of 79 measurements of which 38 were obtained in a single night in order to measure the Rossiter–McLaughlin effect as the planet passed in front of its host star (Gillon et al. 2009). Although the initial discovery paper reported a non-zero orbital eccentricity of 0.054 ± 0.017, this finding was subsequently disputed by Husnoo et al. (2012). These authors reanalyzed the same data while excluding the single night corresponding to the Rossiter measurement and found a best-fit orbital eccentricity of 0.04 ± 0.02. For HAT-P-19b the radial velocity data constrain the orbital eccentricity to a value of 0.07 ± 0.04, and for WASP-39b the orbital eccentricity was fixed to zero in the published fits (Faedi et al. 2011). Spitzer secondary eclipse observations provide a significantly tighter constraint on the value of $e\mathrm{cos}\omega$, where e is the orbital eccentricity and ω is the longitude of periastron. These eclipse times can be used to further refine the global orbital eccentricity when fitting radial velocity data (e.g., Knutson et al. 2014b). We find that the observed secondary eclipse times for these four planets are consistent with circular orbits, and list the corresponding $e\mathrm{cos}\omega$ values in Table 3.

Table 3.  Best Fit Eclipse Parameters

Target	Band (μm)	${F}_{p}/{F}_{*}$ (%)	${F}_{p}/{F}_{*,\mathrm{avg}}$ (%)a	T${}_{\mathrm{bright}}$ (K)a	Tsb	O–C (d)c	e cos ω
HAT-P-19b	3.6	0.071 ± 0.020	0.062% ± 0.014%	${1090}_{-69}^{+61}$	5999.5293 ± 0.0074	−0.0062 ± 0.0075	−0.0015 ± 0.0021
	3.6	0.053 ± 0.021			6933.5697 ± 0.0071	−0.0013 ± 0.0076
	4.5	0.066 ± 0.019	0.062% ± 0.016%	${914}_{-71}^{+62}$	5815.1239d
	4.5	0.053 ± 0.028			6929.5614d
WASP-6b	3.6	0.094 ± 0.019	⋯	${1235}_{-77}^{+70}$	5196.3606 ± 0.0106	−0.0110 ± 0.0106	−0.0001 ± 0.0020
	4.5	0.115 ± 0.022		${1118}_{-74}^{+68}$	5219.9007 ± 0.0044	−0.0021 ± 0.0044
WASP-10b	3.6	0.103 ± 0.015	0.100% ± 0.011%	${1152}_{-36}^{+34}$	5580.9226 ± 0.0025	−0.1095 ± 0.0025	−0.0552 ± 0.0007
	3.6	0.097 ± 0.017			6910.7917 ± 0.0027	−0.1141 ± 0.0027
	4.5	0.126 ± 0.025	$0.146\%\pm 0.016\%$	${1086}_{-39}^{+38}$	5457.2249 ± 0.0040	−0.0980 ± 0.0040
	4.5	0.158 ± 0.020			6923.1693 ± 0.0020	−0.1077 ± 0.0020
WASP-39b	3.6	0.088 ± 0.015	⋯	${1213}_{-62}^{+58}$	6022.2356 ± 0.0066	−0.0033 ± 0.0068	−0.0007 ± 0.0017
	4.5	0.096 ± 0.018		${1055}_{-65}^{+60}$	6026.2878 ± 0.0058	−0.0005 ± 0.0060
WASP-67b	3.6	0.014 ± 0.018	0.022% ± 0.013%	<1046e	6860.3124d
	3.6	0.032 ± 0.020			6837.2416d
	4.5	0.088 ± 0.024	0.080% ± 0.018%	${1015}_{-74}^{+67}$	6846.4689 ± 0.0035	−0.0026 ± 0.0041	−0.0007 ± 0.0012
	4.5	0.070 ± 0.026			6869.5390 ± 0.0064	−0.0007 ± 0.0068

Notes.

^aWe report the error-weighted mean eclipse depths at 3.6 and 4.5 μm for planets with multiple observations in each bandpass. We note that the averaged 3.6 μm eclipse depth for WASP-67b is formally a non-detection with a significance of 1.7σ. Brightness temperatures are calculated using a PHOENIX stellar model interpolated to match the published temperature and surface gravity for each star. ^bBJD${}_{\mathrm{UTC}}$—2,450,000. To convert from UTC to TDB time standards, simply add 66.2 s to eclipse times obtained between 2009–2012, and 67.2 s to eclipse times from 2014. ^cObserved minus calculated eclipse times, where we have accounted for the uncertainties in both the measured and predicted eclipse times as well as the added delay from the light travel time in the system. This delay is 61 s, 54 s, 49 s, 63 s, and 66 s for HAT-P-19b, WASP-6b, WASP-10b, WASP-39b, and WASP-67b, respectively. ^dWe allow the eclipse times in this bandpass to vary as free parameters in our fit, but we use the error-weighted mean orbital phase and corresponding uncertainty from the other bandpass to place a prior constraint on their values. This reduces the uncertainty in the measured eclipse depths for visits where the eclipse depth is not detected at a statistically significant level (>3σ; see Section 2.2 for more details). ^e2σ upper limit based on the error-weighted average of the two 3.6 μm eclipse measurements.

Download table as:  ASCIITypeset image

We next compare the measured Spitzer secondary eclipse depths for each planet to the predictions of standard one-dimensional atmosphere models. We consider two classes of models, based on Burrows et al. (2008) and Fortney et al. (2008). Both models assume that the chemistry of these atmospheres is in local thermodynamic equilibrium and parameterize the unknown recirculation of energy to the night side.

Table 4.  Updated Orbital Solution for WASP-10 from Radial Velocity Fit

Parameter	Value	Units
Pb	3.0927278 $\pm \ 2.8e-06$	days
${T}_{{\rm{conj}},b}$	2454664.03809 $\pm \ 4.7e-05$	${\mathrm{BJD}}_{\mathrm{TDB}}$
Kb	556.7 ${}_{-5.6}^{+5.4}$	m s−1
$e\mathrm{cos}\omega$	−0.054 ± 0.001	⋯
$e\mathrm{sin}\omega$	0.027 ${}_{-0.007}^{+0.0071}$	⋯
eb	0.0608 ${}_{-0.003}^{+0.0036}$	⋯
${\omega }_{b}$	153.6 ${}_{-5.7}^{+6.2}$	degrees
${\gamma }_{1}$	15.3 ± 2.7	m s−1
${\gamma }_{2}$	−64 ± 41	m s−1
$\dot{\gamma }$	−0.037 ± 0.013	m s−1 day−1
$\ddot{\gamma }$	−5.94e-05 ${}_{-9.9e-06}^{+1e-05}$	m s−1 day−2
jitter	6.3 ${}_{-1.3}^{+1.7}$	m s−1

Note. These fits used the transit ephemeris from Barros et al. (2013) and the secondary eclipse times reported in this paper as priors on the fit in order to obtain improved constraints on the orbital solution. Radial velocity data taken from Christian et al. (2009) and Knutson et al. (2014b), with the addition of two new unpublished RV measurements from Keck/HIRES. We allow the zero points of each radial velocity data set (${\gamma }_{1}$ and ${\gamma }_{2}$ to vary independently in our fits, and describe the radial velocity trend with a linear ($\dot{\gamma }$) and quadratic ($\ddot{\gamma }$) function of time. The reference epoch for γ, $\dot{\gamma }$, $\ddot{\gamma }$ is 2455615.0.

Download table as:  ASCIITypeset image

3.2.1. Constraints on Atmospheric Metallicity and Day–Night Circulation

We use the Fortney et al. (2008) models to explore a range of atmospheric metallicities from 1–1000× solar. For planets at these temperatures, an increase in the atmospheric metallicity results in a decrease in the relative amount of methane as compared to CO and CO₂. Because methane absorbs in the 3.6 μm Spitzer bandpass while CO and CO₂ absorb in the 4.5 μm bandpass, metal-rich models will have correspondingly strong emission at 3.6 μm and weak emission at 4.5 μm as compared to their solar metallicity counterparts. At 1× solar, the pressure–temperature profiles in the models presented here are in radiative equilibrium and are self-consistent between the predicted equilibrium chemistry mixing ratios, the opacity of these molecules, and the the absorption and emission of flux. At 1000× solar, the models use the radiative-equilibrium pressure–temperature profiles from a 50× solar model where we have recalculated the chemical equilibrium and spectra at 1000× solar metallicity. Although not fully self-consistent, these models are a reasonable approximation for the appearance of a very metal-enhanced atmosphere. We test the validity of these models in the high metallicity regime by comparing them to the measured secondary eclipse depths for GJ 436b (Lanotte et al. 2014). This planet appears to be best-matched by a model with 1000x solar metallicity in good agreement with previous studies (e.g., Moses et al. 2013), but we are not able to fully match the steep slope across the 3.6 and 4.5 μm bands. We expect that models including disequillibrum chemical process and/or non-solar carbon to oxygen ratios (e.g., Moses et al. 2011; Visscher & Moses 2011) could provide an improved match to these observations, and leave this as an exercise for future modeling studies.

Figures 9 and 10 compare our measured secondary eclipse depths for the five planets in this study to 1× and 1000× solar metallicity models with either efficient or inefficient redistribution of energy to the night side. The relative efficiency of recirculation is implemented by varying the amount of flux incident at the top of the one-dimensional atmospheric column, and is therefore equivalent to a change in the planet's dayside albedo. For the inefficient case, the incident flux is set equal to the dayside average value, while for the efficient case the incident flux is calculated as the average over the entire surface. For four of the five new planets, our observations are best-matched by models with efficient recirculation of energy to the night side but we are unable to distinguish between low and high metallicity atmosphere models with high statistical significance. WASP-10b appears to be the exception to this trend, and is best matched by models with inefficient recirculation and a relatively low atmospheric metallicity. It is also by far the most massive planet in our sample, a fact that we discuss in more detail below. Although we do not show it in these plots, there is a published Ks (2.1 μm) secondary eclipse measurement available for this planet with a value of 0.137% ± 0.016% (Cruz et al. 2015). The measured planet–star flux ratio in this bandpass is 3–10× larger than the predictions of our models, and the reported center of eclipse phase differs from our measurement by 2.4 hr (54σ), comparable to the eclipse duration. It is therefore unlikely that this ground-based measurement corresponds to the same eclipse signal detected in our Spitzer photometry, and in fact the authors note that their reported uncertainties are likely underestimated due to the unknown contribution of systematic noise sources to their error budget.

Zoom In Zoom Out Reset image size

Given that we obtain relatively tight constraints on the atmospheric properties of WASP-10b, it is interesting to ask whether or not it is representative of other Jovian planets in this temperature range. A search of the literature indicates that there are two other planets (HAT-P-20b and WASP-8b) with measured Spitzer secondary eclipse depths at 3.6 and 4.5 μm, masses greater than that of Jupiter, and predicted equilibrium temperatures less than 1200 K. We compare the observations for these two planets to the same set of models as before, and plot the results in Figures 9 and 10. For HAT-P-20b the data are best matched by the solar metallicity model with inefficient circulation of energy to the night side, consistent with the conclusions from Deming et al. (2015) and Triaud et al. (2015). For WASP-8b (Cubillos et al. 2013; Deming et al. 2015) we also prefer models with inefficient energy recirculation and solar metallicity, but this model still under-predicts the measured flux in the 3.6 μm band. The unusual nature of WASP-8b's dayside emission spectrum is discussed in detail in Cubillos et al. (2013), and may be related to its eccentric orbit (e = 0.30). The only other eccentric planets in this sample are GJ 436b and HAT-P-20b, which have orbital eccentricities of 0.15 and 0.016, respectively (Knutson et al. 2014b). Despite this discrepancy, the data for these three planets suggest a clear pattern: planets with masses greater than that of Jupiter are best described by models with solar metallicity and relatively inefficient day–night circulation.

Zoom In Zoom Out Reset image size

We next consider whether or not the smaller planets in our sample show a different behavior. As noted earlier, HAT-P-19b, WASP-6b, WASP-39b, and WASP-67b are all best-matched by models with relatively efficient day–night circulation and are consistent with a range of atmospheric metallicities. We identify two additional planets in the literature with measured Spitzer secondary eclipse depths at 3.6 and 4.5 μm, masses greater than that of Jupiter, and predicted equilibrium temperatures less than 1200 K. As before, we compare the observations for these planets to our standard set of models, and plot the results in Figures 9 and 10. For WASP-80b (Triaud et al. 2015), we find the data require a model with relatively efficient day–night circulation and an atmosphere with at most a moderately enhanced metallicity. For GJ 436b we obtain the best match using a model with efficient atmospheric circulation (this may be related to the planet's non-zero orbital eccentricity) and a high atmospheric metallicity as discussed above.

3.2.2. Temperature Inversions, Clouds, and Other Caveats

Although the models described above suggest the presence of some intriguing patterns, it is important to note that with measurements in just two bandpasses we cannot provide unique constraints on the atmospheric chemistries of these planets. In addition to the potential effects of disequilibrium chemistry and non-solar C/O ratios, the vertical thermal profiles of these planets can also alter the shapes of their dayside emission spectra. We explore the potential effects of temperature inversions using a suite of models from Burrows et al. (2008), which allow us to artificially create atmospheric temperature inversions with varying strengths. Temperature inversions are created by adding a generalized gray absorber in the stratosphere above 0.03 bars with an absorption coefficient κ of 0.1 cm² g⁻¹. A dimensionless parameter P_n describes the efficiency of energy redistribution, with P_n = 0.0 indicating recirculation on the dayside only, and P_n = 0.5 indicating energy uniformly distributed between the two hemispheres. In these models the redistribution is implemented via a heat sink located at pressures between 0.1 and 0.01 bars.

As shown in Figure 11, our measured eclipse depths for for four of the five planets presented in this paper require models with relatively efficient day–night circulation but we are unable to distinguish between models with and without a temperature inversion at a statistically significant level. As with the previous set of models, WASP-10b is best described by models with relatively inefficient day–night circulation. We note that the solar metallicity models with no inversion, which should be identical to the Fortney solar metallicity models, in fact display significantly stronger molecular absorption bands. This difference has been consistently observed in models for planets across a wide range of temperatures (e.g., Deming et al. 2011; Todorov et al. 2012, 2013; O'Rourke et al. 2014; Shporer et al. 2014), and may be due to either the use of different opacity tables, a change in their pressure–temperature profiles due to differences in their method for incorporating heat redistribution to the planet's night side, or some combination of the two.

Zoom In Zoom Out Reset image size

We also consider the possibility that the dayside emission spectra of these planets may be affected by the presence of high altitude clouds or hazes, which have been observed in the transmission spectra of several cool gas giant planets (Jórdan et al. 2013; Line et al. 2013; Knutson et al. 2014a; Nikolov et al. 2015). A range of condensate clouds are possible at these temperatures and pressures (for a recent review see Marley et al. 2013), some of which may occur high enough up in the atmosphere to have a detectable effect. If the atmospheres of these planets have significant methane, then it may also be possible to form a layer of hydrocarbon soot via photochemistry in the upper layers of the atmosphere (Miller-Ricci Kempton et al. 2012; Morley et al. 2013). Published optical transmission spectra for WASP-6b indicate that it does indeed possess a high altitude cloud layer (Jórdan et al. 2013; Nikolov et al. 2015), but there are no published transmission spectra available for HAT-P-19b, WASP-10b, WASP-39b, or WASP-67b. If these planets possess a cloud layer that is optically thick when viewed during secondary eclipse, we would expect such a cloud layer to suppress the strengths of the observed molecular absorption or emission features, leading to a more blackbody-like spectrum. HAT-P-19b has a dayside emission spectrum that is moderately inconsistent with that of a pure blackbody, which should produce a planet–star flux ratio that increases consistently toward longer wavelengths. WASP-6b, WASP-10b, WASP-39b, and WASP-67b are closer to a blackbody shape, and we therefore cannot exclude the possibility that their dayside emission spectra are affected by clouds. However, we note that the transit technique is sensitive to more tenuous clouds and hazes than the secondary eclipse technique, as the slant optical path length traveled by light transmitted through the planet's atmosphere is longer than the radial path traveled by light emitted by the planet (Fortney 2005). GJ 436b serves as a useful case in point, as it has an apparently featureless infrared transmission spectrum (e.g., Knutson et al. 2014a) but still exhibits strong spectral features in its dayside emission spectrum (e.g., Stevenson et al. 2010; Lanotte et al. 2014).

Although our constraints on the atmospheric properties of planets in individual cases are limited by the signal to noise of our measurements, it is possible that the aggregate sample may display trends that are not apparent on a case-by-case basis. We therefore examine the properties of the planets in our sample using a series of model-independent metrics that allow us to construct a view of the global properties of this population. We begin by calculating the ratio of the 4.5–3.6 μm brightness temperatures for each planet; this ratio is independent of both the planet–star radius ratio and the effective temperatures of the planet and star, and instead simply reflects the relative shape of each planet's emission spectrum as determined by the presence of any molecular absorption or emission features. If we assume that molecular features are observed in absorption (i.e., none of the planets have dayside temperature inversions), we would expect the ratio of the 4.5–3.6 μm brightness temperatures to vary primarily as a function of atmospheric metallicity. In this scenario metal-rich planets will have relatively small brightness temperature ratios due to strong absorption from CO at 4.5 μm and weak absorption from methane at 3.6 μm, while metal-poor planets will have correspondingly large ratios due to strong methane absorption and weak CO absorption.

We plot the measured brightness temperature ratios of the nine planets with published Spitzer secondary eclipse depths and effective temperatures cooler than 1200 K as a function of planet mass and equilibrium temperature in Figure 12. Interestingly, we see no evidence for a correlation with planet equilibrium temperature, but instead see tentative evidence for an increase in this ratio as planet mass increases from that of Neptune up to Jupiter. This is consistent with the hypothesis that lower mass planets are more likely to have metal-rich atmospheres, while planets with masses comparable to or larger than that of Jupiter should have envelopes that reflect the (near-solar) metallicity of their local region of the protoplanetary disk. We obtain a nearly identical plot if we simply take the ratio of the measured eclipse depths in these two bands, indicating that the varying effective temperatures of the planets and their host stars have a relatively small effect on the measured eclipse depths as compared to changes in the observed molecular absorption or emission features.

Zoom In Zoom Out Reset image size

We evaluate the statistical significance of this trend by taking the error-weighted average brightness ratio for planets with masses less than 0.3 ${M}_{\mathrm{Jup}}$ and comparing this ratio to the error-weighted average ratio for planets with masses greater than 0.5 ${M}_{\mathrm{Jup}}$. We exclude GJ 436b and WASP-67b from this calculation, as we only have lower and upper limits on their brightness ratios, respectively. We also exclude WASP-8 in light of its large orbital eccentricity, which may affect the observed properties of its dayside atmosphere. We find that planets with masses less than 0.3 ${M}_{\mathrm{Jup}}$ have an average brightness temperature ratio that is 2.1σ smaller than the average ratio for the more massive planets in our sample. We conclude more observations will be required in order to determine whether or not the observed trend is statistically significant.

We next evaluate the relative efficiency of circulation between the day and night sides of these planets by comparing their observed versus predicted dayside emission temperatures as described in Cowan & Agol (2011) and Schwartz & Cowan (2015). For this exercise we consider both the subset of planets with relatively low irradiation levels including the five planets from this study, and also the full set of fifty transiting gas giant planets with broadband infrared eclipse depths detected with greater than 3σ significance. We calculate the irradiation temperature, ${T}_{0}={T}_{*}\sqrt{{R}_{*}/a}$, for each planet based on published values for the stellar effective temperature and the planet's scaled semi-major axis. This temperature corresponds to the predicted temperature at the planet's substellar point assuming an albedo of zero and no redistribution of energy. The maximum hemisphere-averaged dayside effective temperature is then ${(2/3)}^{1/4}$ of this value, and the temperature corresponding to complete redistribution of energy across the entire planet surface would be ${(1/4)}^{1/4}$ of this value.

We calculate the effective dayside temperature measured for each planet using all available infrared ($\lambda \gt 1\;\mu {\rm{m}}$) eclipse depth measurements. We convert the set of individual dayside brightness temperatures across all bandpasses into a combined dayside effective temperature using two different schemes: error-weighted mean (wavelengths with small uncertainties contribute most to the estimate) and power-weighted mean (wavelengths that capture large fractions of the planetary emission contribute most to the estimate). The uncertainty in the averaged dayside effective temperature and in its ratio with the planet's irradiation temperature is then estimated using a 1000-step Monte Carlo method with equal numbers of trials drawn from each scheme (see Schwartz & Cowan 2015, for more details). We use the ratio of the measured vs predicted (irradiation) dayside temperatures to determine whether or not the correlation between planet mass and day–night circulation efficiency inferred in our cool planet sample is consistent with the larger sample of published secondary eclipse data for transiting gas giant planets with a range of incident flux levels. We plot this ratio as a function of planet mass in Figure 13, and find that our new observations clearly indicate the existence of a correlation between planet mass and day–night circulation efficiency for planets with lower irradiation levels. This correlation also appears to be consistent with the published measurements for more highly irradiated planets, although there are fewer low-mass planets in this sample. We conclude that lower-mass planets may be more reflective and/or better at transporting heat to their night sides. Although the current data do not allow us to distinguish between these two hypotheses, the combination of visible-light secondary eclipse observations and infrared phase curve measurements could provide a definitive answer to this question in future studies.

Zoom In Zoom Out Reset image size

Lastly, we create a color magnitude plot for the five planets in this study that allows us to compare them to ultracool brown dwarfs and directly imaged planets spanning the same range of effective temperatures using the same procedure as described in Triaud (2014) and Triaud et al. (2014). In order to maintain consistency with the other planets in this diagram, we obtain stellar parameters for the five planets presented in this paper from TEPCAT (Southworth 2011) and compute their absolute magnitudes using relations by Flower (1996) and Torres & Giménez (2010). Distance moduli were calculated by comparing them to visual magnitudes, which we take from the UCAC4/APASS catalog (Zacharias et al. 2013). Distance moduli, derived distance, visual and absolute magnitudes are provided in Table 5 and visually represented in Figure 14. We calculate these quantities using WISE 1 and WISE 2 apparent magnitudes (Cutri et al. 2013), which are a good match for the Spitzer 3.6 and 4.5 μm bands (Triaud et al. 2014), and list the resulting values in Table 5.

Zoom In Zoom Out Reset image size

Table 5.  Distance Modulus, Distance, Apparent and Absolute Magnitudes for HAT-P-19b, WASP-6b, WASP-10b, WASP-39b, and WASP-67b

	Distance Modulus	Distance	Apparent Magnitude		Absolute Magnitude
	(m − M)	(pc)	${m}_{[3.6]}$	${m}_{[4.5]}$	${M}_{[3.6]}$	${M}_{[4.5]}$
HAT-P-19b	6.73 ± 0.25	222 ± 28	18.51 ± 0.23	18.57 ± 0.31	11.79 ± 0.53	11.85 ± 0.55
WASP-6b	6.65 ± 0.12	214 ± 12	17.84 ± 0.22	17.69 ± 0.21	11.19 ± 0.38	11.04 ± 0.36
WASP-10b	5.58 ± 0.22	131 ± 14	17.43 ± 0.12	17.09 ± 0.12	11.85 ± 0.34	11.50 ± 0.33
WASP-39b	6.67 ± 0.12	216 ± 12	17.80 ± 0.19	17.73 ± 0.21	11.12 ± 0.33	11.06 ± 0.35
WASP-67b	6.64 ± 0.14	213 ± 14	>18.32	17.82 ± 0.25	>11.68	11.18 ± 0.37

Download table as:  ASCIITypeset image

As shown in Figure 14, all five planets have similar magnitudes in this diagram, albeit with varying colors. Along with GJ 436b (Lanotte et al. 2014) and WASP-80b (Triaud et al. 2015), these five planets have the smallest absolute magnitudes of any extrasolar planets observed to date and demonstrate the benefit of extending the current Spitzer observations to smaller and cooler planets. Three of the five planets appear to be bluer than a blackbody with the same absolute magnitude, placing them near the late L dwarfs in this diagram. WASP-6b and WASP-67b lie on or near the blackbody line, where the 2σ upper limit on the 3.6 μm eclipse depth for WASP-67b indicates that it is almost certainly redder than a blackbody. As discussed in Triaud et al. (2014), irradiated exoplanets have optical properties that appear to be more diverse than those of self-luminous brown dwarfs, as demonstrated by the varying colors of the planets in our sample. Although part of this variation likely reflects intrinsic variations in surface gravities and metallicities (e.g., Zahnle & Marley 2014), our analysis in Section 3.2 suggests that disequilibrium chemistry, cloud formation, and variations in atmospheric circulation patterns may also contribute to the observed diversity of dayside emission spectra for the transiting planet sample.