A Case for an Atmosphere on Super-Earth 55 Cancri e

For our analysis, we use the same photometric data as Demory et al. (2016a), taken by the Spitzer Space Telescope Infrared Camera (IRAC 2) at 4.5 μm. Observations during the primary transit and occultation of 55 Cancri e were taken between 2013 June 15 and July 15. A total of 4,981,760 frames were obtained with an integration time of 0.02 s. The frames were taken during eight observing sessions of 9 hr each (half of the planet's orbital period), yielding a total observation time of 75 hr.

Once the observations were obtained from Spitzer, they needed to be converted into a photometric time series from which we can perform phase-curve analysis and modeling. The raw photometric time series was computed from the individual frames by Demory et al. (2016a) using methods outlined in Demory et al. (2011). These raw data were subject to noise from the IRAC detector itself as well as correlated noise from Spitzer as it moves during observation. This noise was removed from our data set by Demory et al. (2016a) to produce a final time series of 30 s bins with an average error that we calculate to be ∼360 parts per million.

After the raw frames were converted into photometric time-series data, we then binned the data to  200 bins for improved visual inspection. The flux value F_p/F_⋆ relative to the in-occultation stellar flux for each bin was computed by averaging the 30 s flux values for each sample. Uncertainties in each bin were computed by dividing the standard deviation of the sample by the square root of the total sample size per bin, yielding an average bin error of ∼56 parts per million. We then phase-folded our time-series data and removed data taken during primary transits that are not accounted for in the phase-curve model we implemented. While not fitting the primary transits, we use the constraints from the primary transits on planetary radius, semimajor axis, and impact parameter in the subsequent phase-curve modeling (see Section 2.3). Our photometric time series is shown in Figure 1.

Zoom In Zoom Out Reset image size

Several effects beyond those of planetary atmospheric processes can determine the properties of an occultation phase curve. For example, stellar variability brought about by magnetic field interactions with the planet can cause star spots with a period similar to the planet's orbital period (Shkolnik et al. 2008). Additionally, gravitational interactions between a star and its planetary companion can produce effects such as relativistic beaming and tidal ellipsoidal distortion modulations that can alter the phase curve (e.g., Shporer & Hu 2015; Shporer 2017). These processes do not produce significant effects on the infrared light curve for 55 Cancri e (e.g., beaming modulation ∼1 ppm, ellipsoidal modulation ∼0.6 ppm; Demory et al. 2016a). We therefore make no corrections to our sample to account for these effects. Furthermore, the orbit of 55 Cancri e can be approximated as circular (e ≈ 0.040 ± 0.027, see Baluev 2015), allowing us to use a model fitting method in line with that of Hu et al. (2015).

One must be cautious when binning time-series photometric data for phase-curve analysis to prevent loss or distortion of information (e.g., Kipping 2010). We have additionally computed model fits for the same data with a much smaller bin size (the native 30 s bins), and with a larger bin size similar to the phase curve in Demory et al. (2016a; a total of ≈70 bins of in- and out-of-transit data). The constraints on the jump parameters are identical to the 200-bin data set to within uncertainty limits, confirming that our model fit does not depend on the binning of the data. We therefore proceed with our analysis using only results from the ≈200-bin data set described above.

To perform our model fit to the photometric data, we use a semianalytical method outlined in Hu et al. (2015). The general model performs a fit to the occultation light curve of a planet and also computes three independent light curves representing flux contributions from symmetric reflection, asymmetric reflection (e.g., due to patchy clouds), and thermal emission. It also uses an existing model of secondary transit (Mandel & Agol 2002) from which it can derive the phase-curve occultation depth, phase amplitude, and phase offset.

As shown in Hu et al. (2015), the scaled total reflective component (i.e., symmetric + asymmetric) of the transiting planet $\tfrac{{F}_{R}}{{F}_{\star }}$ is computed as follows

$$\begin{eqnarray}&&\displaystyle \frac{{F}_{R}}{{F}_{\star }}={A}_{G}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where A_G is the uniform geometric albedo of the planet, R_p is the planet radius, and a is the semimajor axis of the planet's orbit. This means that the contribution of the reflective components to the planet's phase curve is largely governed by the magnitude of ${\left(\tfrac{{R}_{p}}{a}\right)}^{2}$. For 55 Cancri e, using values for R_p and a found in Demory et al. (2016a), we find that this value is 2.77 × 10⁻⁵, or ≈28 parts per million (p.p.m.). This value, even with a geometric albedo of 1, would thus contribute no more than 28 p.p.m. to our final model fit, which takes on much larger flux fluctuations ∼200 p.p.m. The data, with a dispersion of 56 ppm, are too noisy to detect any reflected light component. It is also not possible to distinguish between uniform and patchy reflective clouds or surfaces, because the asymmetric reflection component would make a contribution much smaller than 28 ppm to the light curve. We therefore assume for the purposes of this paper a purely symmetric reflection component and take the asymmetric component of the reflection contribution to be zero at all phase values.

The general model fit computes the posterior distributions of three to five parameters from the phase curve. The first is the Bond albedo (A_B), which is a measurement of the fraction of incident light on the planet that is scattered back into space. The second, referred to as the heat redistribution efficiency (), is defined as

$$\begin{eqnarray}&&\epsilon =\displaystyle \frac{{\tau }_{\mathrm{rad}}}{{\tau }_{\mathrm{adv}}},\end{eqnarray} \tag{ 2 }$$

where τ_rad is the radiative timescale and τ_adv is the advective timescale (Cowan & Agol 2011; Hu et al. 2015). This parameter describes how well heat is transported from the planet's dayside to nightside. When $| \epsilon | \gg 1$ heat transport is more efficient than radiative cooling and the longitudinal variation of temperature will be small, corresponding to giant planets in the solar system; when $| \epsilon | \ll 1$ transport is minimal and the planet will be in local thermal equilibrium, leading to large day–night contrast (Cowan & Agol 2011). The sign of is the same as the sign of τ_adv, which indicates the direction of transport and the thermal phase shift: when  > 0, the transport is eastward and the peak of the thermal emission appears prior to the occultation; when  < 0, the transport is westward and the peak of the thermal emission appears after the occultation. The magnitude of the radiative timescale depends on the heat capacity and the temperature, and the magnitude of the advective timescale is determined by the speed of fluid movement (atmosphere or molten lava), which is discussed further in Section 4.

The third parameter is the greenhouse factor (f), a measurement of the extent to which infrared radiation is absorbed by greenhouse gases in the atmosphere. It is defined as the ratio between the brightness temperature of the photosphere and the equilibrium temperature. The last two parameters are the cloud condensation temperature (T_c), under which the atmosphere precipitates and forms patchy clouds, and the reflectivity boosting factor (κ), which describes the proportional increase in the reflectivity of the clouds in the atmosphere. The first three parameters (A_B, , f) are computed for each model, while the last two (T_c and κ) are only relevant in the scenario in which the planet has a patchy cloud covering (see Hu et al. 2015 for more details). Because we are assuming the surface of the planet to be symmetrically reflective (i.e., without a patchy cloud or lava lake), T_c and κ are not used by the phase-curve model. We are thus performing a model fit that depends on three free parameters: Bond albedo, heat redistribution efficiency, and greenhouse factor.

The methods developed and outlined in Hu et al. (2015) use a Markov Chain Monte Carlo (MCMC) to compute both the fitted parameters and their posterior distributions. We first input the allowed ranges within which the MCMC can sample for each parameter. The Bond albedo is defined as a fraction that naturally ranges in [0, 1]. The heat redistribution is allowed to range between −100 and 100, where the sign is determined by the direction of the photosphere rotation with respect to the planet's orbit. For the atmospheric greenhouse factor f, we allowed a range of [1, 2] with 1 corresponding to no greenhouse effect and 2 corresponding to a thermal photosphere temperature two times greater than the planet's equilibrium temperature. Allowing f > 2 is not necessary because the temperature of the planet does not have to double even for A_B approaching 1 (see Section 3.2). These input ranges are wide enough and do not cause distortion in the parameters' posterior distributions; we confirm this by inspecting the output of the Markov chains, which shows that the posterior distributions are narrower than the prior ranges (Figures 3 and 4).

In addition to the allowed ranges for the fitted parameters, we input values for system parameters that the algorithm uses to compute the model phase curves. We use an empirical and absolutely calibrated infrared spectrum of 55 Cancri (Crossfield 2012) to compute the stellar flux in the IRAC 4.5 μm bandpass, and also use the scaled semimajor axis (a/R_⋆), the ratio between the planet's radius and the star's radius (R_p/R_⋆), and the impact parameter (b) constrained by primary transit observations (Demory et al. 2011, 2016a; von Braun et al. 2011). We use the value of R_p/R_⋆ measured from the same observation as the phase curve (Demory et al. 2016b), which is 365 ± 25 ppm, corresponding to R_p/R_⋆ = 0.0191 ± 0.0007, in our analysis. The stellar luminosities and orbital parameters are known to high precisions, and they make negligible contributions to the uncertainties of our fitted parameters. The planetary radius is more uncertain and we thus propagate its uncertainty to the fitted parameters in Section 3.2.

Next we computed two separate Markov chains, each with 100,000 steps representing different combinations of the fitted parameters. The first half of each chain, which we refer to as the "burn-in" period, was discarded from our final results. We then confirmed the robustness of our final results by verifying their convergence (i.e., assuring R < 1.01 for all parameters, as described in Gelman & Rubin 1992).

The fitted and derived system parameters for 55 Cancri e are listed in Table 1. The values for these parameters were computed by taking the median of the posterior distribution generated by the MCMC analysis (see Section 2.3). The final model fit and residuals computed from these parameters are plotted against our occultation data in Figure 2. From the model fit, we can see that the phase curve of 55 Cancri e is dominated by thermal emission with no significant contribution from the reflected light. We can also see that the fit to the light curve is not symmetric about the secondary transit. Instead, there is a phase curve maximum that occurs pre-occultation at an offset of ≈−34°.

Zoom In Zoom Out Reset image size

Table 1.  Fitted and Derived Parameters for 55 Cancri e^a

Parameter	Value
Derived Parameters
Phase amplitude (ppm)	${127}_{-14}^{+12}$
Eclipse depth (ppm)	${181}_{-19}^{+15}$
Phase offset (degree)b	$-{34}_{-5}^{+4}$
Fitted Parameters
Bond albedo AB	${0.43}_{-0.29}^{+0.32}$
Heat redistribution c	${1.47}_{-0.25}^{+0.30}$
Greenhouse factor f	${1.31}_{-0.13}^{+0.28}$
Calculated Surface Temperatures
Maximum hemisphere-averaged temperature (K)	${2709}_{-159}^{+129}$
Minimum hemisphere-averaged temperature (K)	${1613}_{-131}^{+118}$
Average dayside temperature (K)	${2573}_{-153}^{+120}$

Notes.

^aThe best-fit model has χ²/dof of 1.14. ^bDefined to be positive for post-occultation maximum. ^cDefined to be positive for eastward-traveling winds in a synchronously rotating reference frame.

Download table as:  ASCIITypeset image

The phase-curve model from Hu et al. (2015) generates the amplitude, occultation depth, and offset based on the fit to the observed data (see Table 1). Errors for all values in Table 1 were computed using the values 1 standard deviation above and below the mean, corresponding to z scores of ±1 for a standard normal distribution. Here we measure our phase-curve amplitude to be ${127}_{-14}^{+12}$ p.p.m. This value is ≈2σ smaller but better constrained than the amplitude found in Demory et al. (2016a). For our eclipse depth, we derive a value of ${181}_{-19}^{+15}$ p.p.m., which is 1σ greater than the value derived in Demory et al. (2016a). Upon inspection of the model phase curves in Demory et al. (2016a), the best-fit models look similar at the peak ($\tfrac{{F}_{p}}{{F}_{\star }}$ around ∼200 p.p.m.), but our best-fit model indicates a minimum out-of-transit value of ∼70 p.p.m., as opposed to ∼50 p.p.m. of Demory et al. (2016a). The difference is mostly due to the different models used in the phase-curve analysis. Demory et al. (2016a) used a single-longitude-band model and a three-longitude-band model to fit the phase curve, and their models have sharp temperature discontinuities between the bands. Our model does not allow for these discontinuities and instead calculates a smooth longitudinal distribution of temperature by solving an energy transport equation (see Hu et al. 2015). Our physically motivated model provides a good fit to the phase curve, and indicates a smaller phase amplitude and a higher nightside temperature than Demory et al. (2016a). This analysis thus highlights the need to use a physically motivated model in the phase-curve analysis compared to the longitudinal band models.

Our third derived parameter is the phase-curve offset, defined to be positive when the phase-curve maximum occurs post-occultation. We find a pre-occultation phase curve offset of $-{34}_{-5}^{+4}$ degrees, which agrees with the offset from Demory et al. (2016a) to within error estimates. Because the phase curve is dominated by contributions from thermal emission (Figure 2), the offset is likely due to a hotspot located east of the substellar point (see Section 4 for more details).

The posterior probability distributions for our fitted model parameters are shown in the top three panels of Figure 3. As we can see, the three derived parameters discussed above are tightly constrained in the scenario to which we fit our model curve. In other words, the phase-curve amplitude, offset, and occultation depth are well constrained to the values listed in Table 1.

Zoom In Zoom Out Reset image size

We also used our derived phase-curve parameters to compute a series of surface temperatures for the planet (Table 1). We used the eclipse depth to compute an average dayside temperature of ${2573}_{-153}^{+120}$ K. From the phase-curve amplitude we calculate a minimum and maximum hemisphere-averaged temperature of ${2709}_{-159}^{+129}$ K and ${1613}_{-131}^{+118}$ K respectively. All of these values agree with those found in Demory et al. (2016a) to within 1σ, but we compute a higher nightside temperature and a smaller day–night temperature contrast of ∼950 K that is more in line with the presence of a convective envelope. The fact that our maximum hemisphere-average temperature is greater than the average dayside temperature is consistent with our phase-curve offset of $-{34}_{-5}^{+4}$ and is likely due to an eastward-shifted hotspot. These implications are discussed further in Section 4.

In addition to the derived parameters discussed in Section 3.1, the model phase curve contains enough information to compute a set of fitted planetary parameters (see Section 2.3). In the case of a homogeneous reflective layer, the model outputs estimates for the planet's Bond albedo A_B, heat redistribution efficiency , and greenhouse factor f. The fitted values for these parameters can be found in Table 1.

The posterior probability distributions of the fitted planetary parameters for 55 Cancri e can be found in the three bottom panels of Figure 3. As we can see, the value for the heat redistribution efficiency (middle panel) is well constrained between 1 and 3 at an estimated value of  = ${1.47}_{-0.25}^{+0.30}$. A positive value for at almost 5σ above zero indicates that the advective frequency of the planet's envelope is nonzero and consequently that material at photospheric pressures above the planet's surface travel eastward in a synchronously rotating frame. The magnitude of is inconsistent with a lava ocean and thus suggestive of a thick atmosphere. This scenario is discussed in more detail in Section 4.

As we can tell by inspection of Figure 3, the Bond albedo and greenhouse factor do not have the same narrow distribution as , and they are not constrained independently. The planet's Bond albedo A_B can take any value smaller than 0.9. Similarly, the atmospheric greenhouse factor f is typically closer to 1.2 but also takes essentially all of the allowed values between 1 and 2. The two parameters are correlated, as shown in Figure 4. This is due to the proportionality from Hu et al. (2015), viz.

$$\begin{eqnarray}&&T\propto f{(1-{A}_{B})}^{1/4},\end{eqnarray} \tag{ 3 }$$

which our model uses to compute the temperature distribution on the planet's surface. As we can see from Equation (2), raising f requires an additional increase in A_B for a given temperature and vice versa. This is in line with what we see in Figure 4, where a lower A_B ≈ 0.1 corresponds to lower f ≤ 1.3 and f slowly leveling off at 2 as A_B approaches 1. Additionally, the fact that f > 2 corresponds to a forbidden Bond albedo greater than 1 further motivates our allowed ranges for f in our model fit of [1, 2].

Zoom In Zoom Out Reset image size

The uncertainty in the planetary radius (8% in terms of (R_p/R_⋆)²) has an impact on the derived temperatures, and thus the greenhouse factor f. We propagate this uncertainty by calculating the brightness temperatures corresponding to the 1σ upper and the lower bounds of the planetary radius, and find that this would increase the uncertainty in the dayside and nightside temperatures by 50 K and 20 K respectively. This in turn causes the uncertainty in the greenhouse factor to increase by 0.04, in addition to 0.28 as tabulated in Table 1. The uncertainty in the planetary radius thus has no impact on the general results.