Confirming Variability in the Secondary Eclipse Depth of the Super-Earth 55 Cancri e

To remove systematic intrapixel sensitivity fluctuations, we utilized the PLD framework presented by Deming et al. (2015). Other methods currently used to correct for the intrapixel effect depend on relating intensity fluctuations to the position of the stellar image on the detector (e.g., BLISS mapping Stevenson et al. 2012, used by D16). As Deming et al. (2015) discuss, the position of the stellar image is a secondary data product, ultimately derived from the intensities of individual pixels. PLD utilizes the primary data, the intensities of individual pixels to correct intrapixel sensitivity variations.

PLD is capable of removing red noise that can frustrate other methods. Ingalls et al. (2016) compared seven different Spitzer decorrelation methods and found that PLD (as well as the BLISS mapping used by D16) was among the three methods that came closest to finding the true eclipse depth when applied to synthetic Spitzer data (however, it should be noted that six of the seven methods compared in Ingalls et al. 2016 produced results that were indistinguishable within the expected measurement error when applied to real and synthetic data). The PLD technique has been used successfully in a number of recent Spitzer analyses (Kammer et al. 2015; Buhler et al. 2016; Fischer et al. 2016; Wong et al. 2016; Dittmann et al. 2017; Kilpatrick et al. 2017). PLD is also the conceptual basis for the highly successful EVEREST code (Luger et al. 2016) that achieves high precision in decorrelating data from the K2 mission.

Other known systematic effects that PLD must correct for include a temporal "ramp" effect, in which the reported intensities of detector pixels change rapidly until they reach a threshold value, and temporal baseline effects, in which the intensity gradually increases or decreases over the course of observations. To account for the transient temporal ramp, we simply avoided data near the start of observations (see Section 2.2). To remove the second effect, we fit visit-long functions to the photometry (linear, quadratic, or exponential) and subtracted them from the photometry. Through testing with the Bayesian Information Criterion (BIC; covered in more detail in Section 4.2.2), we found that the best fits are achieved through the use of visit-long linear baselines. For these observations, we found that not fitting a baseline at all typically results in failure to find a regression solution. As a result, we are confident that linear baselines are the best option for removing temporal baseline effects.

Previous applications of the code used either 9 (Deming et al. 2015) or 12 (Garhart et al. 2018) pixels, in accordance with the number of pixels that captured a significant portion of the stellar flux in those works. Here, we have updated PLD to work with anywhere from 1 to 25 pixels in a 5 × 5 grid surrounding the star. While in principle the code can be used with any set of pixels in this grid, we expect the starlight to be distributed symmetrically about the central pixel upon which the star is placed in the image. This expectation, coupled with the fact that edge pixels in the 5 × 5 box captured  15% of the total flux on average for these observations (see Figure 2), motivated us to use a 5 × 5 pixel grid to decorrelate the photometry.

Zoom In Zoom Out Reset image size

Figure 3 shows a typical application of the regression fit on a 2011 transit for the best binning and aperture size (see Section 2.3.3). In the figure, we show the process of going from raw photometry, to binned fit, to transit model with systematics removed.

Zoom In Zoom Out Reset image size

In the original version of PLD, the central phase of the eclipse was determined through fitting an eclipse model to unbinned data over a range of trial values and selecting the one that minimized χ². We found that this method was unreliable for these observations of 55 Cnc e as the code could have difficulty in finding the small transit/eclipse events in the unbinned data. As a workaround, we set an initial guess for the central phase, based on literature values for the ephemeris (Endl et al. 2012), and ran an instance of the Markov Chain Monte Carlo procedure (MCMC) on binned data to refine the central phase. This value was then used in subsequent instances of the code.

For several of the eclipses, we found that even this workaround could not find a definite central phase, due to the shallow depths of the eclipses. We therefore elected to generate an ephemeris using the transit events, which are several times deeper than the expected eclipse depths. We used this ephemeris to calculate the expected eclipse times/central phases, assuming a zero-eccentricity orbit. Previous works (Nelson et al. 2014; Demory et al. 2016b) found eccentricities consistent with zero, despite the presence of four other planets in the system. We also calculated the tidal circularization timescale for 55 Cnc e using

$$\begin{eqnarray}&&{\tau }_{e}=\displaystyle \frac{4}{63}Q{\left(\displaystyle \frac{{a}^{3}}{{GM}}\right)}^{0.5}\left(\displaystyle \frac{m}{M}\right){\left(\displaystyle \frac{a}{{R}_{p}}\right)}^{5}\end{eqnarray} \tag{ 1 }$$

(Goldreich & Soter 1966), where Q is the planet's specific dissipation function. Assuming Q = 100 (Gladman et al. 1996), we find a tidal circularization time of about 14,000 years, which justifies the assumption of a zero-eccentricity orbit.

With the central phase provisionally determined, we varied the photometric aperture size and binned the photometry over a range of timescales. Central to the PLD method is that it usually fits to data that are binned in time, as discussed in Section 3.1 of Deming et al. (2015). Kammer et al. (2015) also discuss the rationale for binning in time. Briefly, if the scatter in the data is dominated by stochastic noise, it will be impossible to detect correlated noise and match it to our pixel basis vectors. By binning the data, we can diminish short-cadence white noise and expose long-cadence correlated noise. Fitting to binned data involves a loss of information (e.g., due to high-frequency image jitter). However, we find that fitting to binned data produces a better match between the timescale of the fit and the dominant timescale of the intrapixel fluctuations. Other authors (e.g., Deming et al. 2015; Buhler et al. 2016) have also used binned PLD fits very successfully.

As a function of both aperture radius and bin size, we searched for an optimum solution to the sum of the intrapixel effect and the eclipse of the planet. We used three constant aperture radii (see Section 2.3.6) and 112 bin sizes ranging from 2 to 2048 points per bin. For a given trial of bin size and aperture radius, we found the best fit by linear regression (solving Equation (4) from Deming et al. 2015), holding the central phase of the eclipse fixed at the provisional value. We then applied this trial fit to the unbinned data, and we investigated the noise properties of those residuals as a function of timescale. We denote the standard deviation of the unbinned residuals as σ(1). We binned the residuals over 2, 4, 8, 16, etc. points, until the effect of binning would reduce the number of residuals to ≤16 points. We calculated the standard deviation of the residuals for each of these bin sizes, denoted as σ(N).

For a perfect fit that removes all but the photon noise, the standard deviation of the residuals (SDNR) should vary inversely with the square root of the bin size (a slope of −0.5 in log₁₀ space). We list the SDNR values of our fits compared to the photon noise in Table 1. The relation between standard deviation and bin size is known as the Allan variance relation (Allan 1966). In the original PLD method, Deming et al. (2015) adopted the trial fit (i.e., bin size and aperture radius) that produced the minimum chi-squared in the Allan variance relation. We have found it desirable to slightly modify that original criterion, as we now explain.

In choosing the best fit from different aperture radii and (especially) bin sizes, PLD effectively makes a compromise between mitigating the short-term scatter, represented by the standard deviations of unbinned normalized residuals (SDNR), and the longer-term red noise, represented by the standard deviations of the residuals binned to longer timescales. The standard deviation of the residuals itself has an uncertainty (i.e., the error on the error). That uncertainty increases with bin size because larger bins sample the error distribution with fewer points. Hence, the chi-squared of the Allan variance relation is dominated by the smallest bin sizes. For that reason, we found that the original criterion used by Deming et al. (2015; minimum chi-squared in the Allan variance relation) is virtually equivalent to choosing the solution with minimum unbinned scatter (i.e., minimizing SDNR). We therefore found it advantageous to modify the best-fit criterion to choose the solution that has the minimum raw scatter in the Allan variance relation (σAVR), after forcing the relation to pass through σ(1). That is similar to the original criterion from Deming et al. (2015), but with greater weighting for the larger bin sizes. Note that this fitting criterion is used (successfully) in Section 2.3.6 wherein we derived null-test eclipse depths to check our errors.

In order to quantify the improvement of using our new version of the PLD code, we tested our new version versus the original code described by Deming et al. (2015). We used three representative eclipses (2012 January 18 and 2013 June 15 and 29) and derived eclipse depths using the original code and our new code applied to the same photometry files, using the same degree of data binning. We also constrained the central phase of the eclipse to be the same in each code. We found that the new code reduces the single-point photometric scatter in the residuals (data minus fit) from 1.27 times the photon noise to 1.11 times the photon noise. Given the large number of data points (239,000), that degree of reduction in the photometric scatter easily justifies (i.e., via a BIC analysis) the inclusion of the additional basis pixels in the PLD solutions. The average error on the eclipse depth, from the MCMC posterior distributions, is reduced from the 50 ppm produced by the original code to the 39 ppm from our new code. The average eclipse depth from the two codes was the same to within 1.2 times the error of the mean.

In testing these observations, we found that many combinations of bin and aperture size produce comparably good σAVR scores for a given eclipse. The eclipse models associated with those combinations, however, could have eclipse depths that differed by tens of parts per million between each other (consistent with their errors), differences that could be significant when searching for potential variability between such small events. The distribution in eclipse depth for similar σAVR scores is a result of the fact that σAVR is not a rigorous statistical criterion. Instead, it is an attempt to get a "broad-bandwidth" solution by trying to minimize both red and white noise through binning. Because of this, we selected our best-fit eclipse depth for each observation by examining a distribution of solutions, chosen by their σAVR scores.

After running the regression for all bin and aperture sizes, we fit a Gaussian to the distribution of σAVR values, which we refit after removing 3σ outliers. We selected the center of this Gaussian as a cutoff value, and we then focus on solutions with σAVR scores better than the cutoff. There are typically ∼100 solutions left after removing the bad solutions, and the eclipse depths of these best solutions were approximately Gaussian-distributed for each eclipse. To determine the depth that we quote as our result, we fit a Gaussian to this distribution, selecting the combination of aperture and bin size that produced the eclipse depth that is closest to the central value. This approach allowed us to find a solution that is the most stable in the sense of representing the best fits found by the code, and prevented us from reporting an outlier as the best eclipse model.

After finding the best combination of bin and aperture size, we held the binning and aperture constant, and refined the fit using an MCMC procedure (Ford 2005). While we report the depth associated with the regression solution, the MCMC allowed us to explore a high-dimensional parameter space (our models use ∼30 parameters when including the pixel coefficients) which accounted for possible correlations between fitting parameters and defines the error on the eclipse depth via a Gaussian fit to the posterior distribution.

The MCMC utilized the Metropolis–Hastings algorithm with Gibbs sampling. At each step, we adjusted the eclipse depth and selected orbital parameters (for transit fits), temporal baseline, and pixel coefficients. Each chain consists of 500,000 steps, and we ran multiple chains so as to verify convergence through the Gelman–Rubin statistic, R (Gelman & Rubin 1992).

Comparing our error estimates from the MCMC with those found by the best regression solution chosen in Section 2.3.3, we found that the MCMC errors are on average 13% higher than the regression errors for the eclipses.

We also tested initializing the MCMC with different solutions than the one chosen by our method of selecting a representative fit from the distribution of best σAVR scores. These fits are typically within 2 ppm of the "best" solution, but use different aperture and bin sizes. We find that the different initializing solutions did not produce significantly different errors from the solution chosen by our method. We therefore conclude that the errors are stable in the sense that they are not sensitive to the details of the fit.

As stated above, our errors on the eclipse depths are determined using the MCMC simulations described in the previous section. However, deriving realistic errors is critical to the inference that the eclipse is variable, and a key feature of our analysis is an independent check of our derived errors.

To make this check, we identified eight regions in the 2013 eclipse data (two for each of the four eclipses) that were independent of each other and do not overlap with the actual eclipse. In these regions, we expect to derive eclipse depths that are consistent with zero within our error bars, because we know that they do not contain the actual eclipse. A critical test of the errors is that the derived "eclipse" depths for these regions should scatter around zero with a dispersion that is consistent with their errors.

We "fit" these non-eclipse regions following the methodology in Sections 2.3.1–2.3.5, locking the central phase of the model in each case. From this test, we found that if we restrict ourselves to using large photometric apertures (with radii of 3.1, 3.3, and 3.5 pixels, the same used in our eclipse analysis), we indeed find eclipse depths that vary around zero to a degree consistent with their derived errors (calculated through the standard deviation of in-eclipse and out-of-eclipse residuals, added in quadrature). The restriction to the larger aperture sizes makes sense given our choice of using all basis pixels in a 5 × 5 pixel grid surrounding the star: these apertures are well matched to the area of this basis pixel grid. Consequently, we limit ourselves to using these three apertures when fitting the actual transits and eclipses.

The results of the fits to the non-event regions are shown in Figure 4, which plots the distribution of eclipse depth divided by its derived error. If our errors are realistic, that distribution should be consistent with a Gaussian having a standard deviation of unity. The distribution only contains eight samples because we are limited by the amount of independent data regions that are available. Nevertheless, the distribution is approximately centered on zero with a dispersion consistent with unity. Our 1σ error bars, which average to 34 ppm for these data, encompass all of the non-event results except for two. The first of these is within 1.3σ of zero, which we expect with a frequency of once every five measurements. The second is within 1.7σ of zero, which is expected once every 11 measurements.

Zoom In Zoom Out Reset image size

Moreover, performing an Anderson–Darling test for normality (Stephens 1974), we find a p-value of 0.43, implying that our results for the non-event regions are consistent with a Gaussian having a standard deviation of unity. We conclude that our fitting procedure produces realistic errors, suitable for drawing inferences as to the variability of the eclipse.

We generated five different models to explain the pattern of eclipse depths versus time using Levenberg–Marquardt fitting in IDL, which we show in Figure 10. We now describe each of these models in detail.

Zoom In Zoom Out Reset image size

The first model is a flat eclipse depth, i.e., one that is not varying in time. This is the model that we would anticipate for a "normal" planet, as we would expect to see a mostly constant thermal flux from the planet versus time. There is only one free parameter for the flat model, that being the weighted average eclipse depth. We find that a depth of 84 ± 14 ppm best fits all of the eclipses in the flat model.

The second model is similar to the one proposed by D16, with year-to-year variability between the 2012 and 2013 eclipses. We fit constants to our depths from 2012 and 2013 separately to determine this model. The year-to-year model is characterized by three free parameters: an eclipse depth for both 2012 and 2013, and a time of transition. For 2012, we find a depth of 24 ± 20 ppm, and a depth of 138 ± 18 ppm in 2013.

The third model, which we refer to as the "spike" model, uses a flat depth everywhere except for the first two observations in 2013, which we found were the two deepest eclipses. The four free parameters of the spike model are the baseline eclipse depth, the depth of the "spike", and two transition times. For this model, we find a baseline depth of 58 ± 16 ppm, and a depth of 173 ± 26 ppm for the first two observations in 2013.

The fourth model, which we call the "sloped" model, is similar to the third, but with a linearly decreasing slope in 2013. It also has four free parameters: the average 2012 eclipse depth, a time of transition, and a peak and slope value for the 2013 data. It uses the same 2012 value as the year-to-year model, and we perform a linear fit to the 2013 data of the form y = mx + b. We find m = −2.36 ± 1.59 ppm day^–1 and b = 164.7 ± 26.0 ppm.

The final model is a simple sine wave, given by A sin(2πft+δ) + c, where A is the amplitude of the wave, f is the frequency, δ is a phase shift, and c is a constant offset from 0 (i.e., four free parameters). To determine this model, we fit a sine function to our best-fit eclipse depth results using Levenberg–Marquardt fitting, stepping over periods that ranged from half a day to about three years. At each frequency step, we fit for sine wave amplitude, phase offset, and constant offset, and recorded the χ² of each fit. We plot the results in Figure 11, with frequencies converted to periods in days. This model represents a periodic cycling of the planet's eclipse depth over time. Physically, a sine wave could be a useful approximation of any process (e.g., volcanism) that has a characteristic timescale, i.e., it represents a bandwidth-filtered stochastic process.

Zoom In Zoom Out Reset image size

Figure 11 shows that a number of different sine wave periods fit the data similarly well. This is because the data are heavily aliased, as our sampling rate does not permit us to distinguish between sine waves of many different frequencies. In the bottom panel of the figure, we show the sine model generated using the period circled in red in the top panel. Fluctuations of 300 ppm over a 2 day period seem unphysical, yet fitting a sine wave to our aliased data finds that this is the "best" model.

While we again acknowledge that the data are heavily aliased, we continue forward with a sine wave model because future data could better constrain the periodicity (if periodic fluctuation is actually the physical reality of the system). We pursue a sine model with a less extreme period, electing to use the best fit with periodicity greater than 20 days. This fit is circled in blue in Figure 11. We found the model's parameters to be A = 76.46 ± 18.99 ppm, f = 0.179057 days⁻¹, δ = 1.28 ± 0.30, and c = 103.54 ± 14.16 ppm.

As in other works (e.g., Gibson et al. 2010; Haynes et al. 2015), we calculated the BIC (Schwarz 1978) to assess the relative quality of fit between these models. The BIC is calculated through the formula BIC = χ² + k · ln(n), where k is the number of model parameters (one for the flat model, four for the spike model, three for the year-to-year model, four for the sloped model, and four for the sine model) and n is the number of data points (eight in our case). A lower BIC is preferred, with an improvement of 6–10 points representing strong evidence for one model over another, and more than 10 points representing very strong evidence (see Kass & Raftery 1995). We list the BICs for each model in Figure 10.

Based on our analysis, we found that the flat model is clearly ruled out by its BIC score, despite using the fewest number of parameters. The remaining models all have BICs that are within 6 of each other, so we are unable to confidently distinguish between them. Because the year-to-year model has the fewest free parameters out of the variable models, we conclude that it is the best model we currently have for interpreting the planet's eclipse depth variability. While the sine model has the best BIC score, we again emphasize that the fit was performed to heavily aliased data, and more observations would be required to determine if a sine model is realistic.

Following the calculations of Crossfield (2012), we converted the maxima and minima of our five models into planetary brightness temperatures (T_B) using an observed spectrum of 55 Cnc. That paper demonstrated that the observed spectrum gives a significantly lower brightness temperature for 55 Cnc e than would be calculated when treating the star as a blackbody. Using their measured flux density of 3.444 ± 0.006 Jy, we numerically evaluated the following integral for values of T_B:

$$\begin{eqnarray}&&\int {B}_{\lambda }({T}_{B}){S}_{x}(\lambda )d\lambda ={\delta }_{\mathrm{occ}}\displaystyle \frac{{R}_{* }}{{R}_{p}}\int {F}_{\lambda }^{* }{S}_{x}(\lambda )d\lambda ,\end{eqnarray} \tag{ 2 }$$

where B_λ(T_B) is the wavelength version of Planck's law as a function of the planet's brightness temperature, S_x(λ) is the relative system response in Spitzer IRAC Channel 2,⁶ ${F}_{\lambda }^{* }$ is the stellar flux density (converted to W m⁻³ sr⁻¹), and δ_occ is the planet's measured eclipse depth. The maximum and minimum values of T_B for each model can be found in Table 5.

Table 5.  Maximum and Minimum Brightness Temperatures for Each of the Five Models

Model	Minimum TB (K)	Maximum TB (K)
Flat	${1820}_{-147}^{+142}$	${1820}_{-147}^{+142}$
Spike	${1550}_{-188}^{+175}$	${2642}_{-242}^{+240}$
Year-to-Year	${1115}_{-415}^{+266}$	${2331}_{-185}^{+183}$
Sloped	${1115}_{-415}^{+266}$	${2571}_{-245}^{+242}$
Sine	${1163}_{-628}^{+330}$	${2705}_{-252}^{+249}$

Note. Temperatures were calculated using an observed stellar spectrum of 55 Cnc from Crossfield (2012).

Download table as:  ASCIITypeset image

To examine the plausibility of our different models, we are motivated to explore what physical mechanisms could account for the brightness temperature variations suggested by the year-to-year and sine models.

Variability in Stellar Flux. The most likely mechanism for variations in stellar flux on day-to-week timescales is the influence of star spots. Analyzing 11 years of photometric observations of 55 Cnc, Fischer et al. (2008) found measurements covering two consecutive rotational periods (∼80 days) that showed evidence of low-amplitude, short-term variability. This variability was attributed to star spots covering less than 1% of the star's surface, resulting in a brightness amplitude variation of 0.006 mag. While this variation could produce a signal of ∼200 ppm in transit/eclipse observations (Demory et al. 2016a), the fact that Fischer et al. (2008) found 55 Cnc to be a quiet star makes this an unlikely scenario. Furthermore, the star would have to be variable on the timescale of the eclipse itself in order to affect our results. On this basis, we conclude that the star is not responsible for the observed eclipse depth variability.

Asynchronous Rotation with a Hotspot. Because the measured variability is not well accounted for by changes in the stellar flux, the variation must be due to changes in the planet's flux at 4.5 μm. If the variability is truly periodic, one of the easiest physical explanations to visualize is that 55 Cnc e is rotating asynchronously while maintaining a hotspot (e.g., from volcanism) on one side of the planet. In this scenario, we would observe the occultation of the planet's hot, bright side in some observations and its cooler, darker side in others. This would lead to more significant drops in relative flux during some observations than in others in a periodic fashion that matches the planet's spin. While we cannot decisively constrain the planet's rotation rate from our results, Figure 11 suggests periods between 15 and 80 days.

A valid objection to this scenario would be the expected timescale for tidal locking for 55 Cnc e, which would cause us to see the same face of the planet before/after eclipse in every observation. We can estimate this timescale from Equation (9) given in Gladman et al. (1996):

$$\begin{eqnarray}&&{\tau }_{\mathrm{lock}}=\displaystyle \frac{\omega {a}^{6}{IQ}}{3{{Gm}}_{s}^{2}{k}_{2}{R}_{p}^{5}},\end{eqnarray} \tag{ 3 }$$

where ω is the spin rate, a is the radial distance from planet to star, I is the planet's moment of inertia about the spin axis, Q is the planet's specific dissipation function, G is the gravitational constant, m_S is the mass of the star, k₂ is the tidal Love number of the planet, and R_p is the mean radius of the planet. We use the estimate of Q = 100 from Gladman et al. (1996) for satellites in our solar system, and calculate k₂ as

$$\begin{eqnarray}&&{k}_{2}=\displaystyle \frac{3/2}{1+\tfrac{19\mu }{2\rho {{gR}}_{p}}}\end{eqnarray} \tag{ 4 }$$

(Burns 1977), where μ is estimated as 3 × ${10}^{10}\tfrac{N}{m}$ for rocky bodies (Gladman et al. 1996), ρ is the average planetary density (using M_p from Demory et al. 2016a and R_p from this work), g is the surface gravity of the planet, and R_p is the planetary radius. Due to the small semimajor axis and high bulk rigidity of 55 Cnc e, de-spin times are on the order of hundreds of years, which makes unsynchronized rotation an unlikely explanation for variability.

Alternatively, it is possible that the planet could have become captured in a higher-order spin–orbit resonance as it spun down. Mercury, with its 3:2 resonance (Pettengill & Dyce 1965), is an example of such a configuration in our own solar system. Further observations of 55 Cnc e will be required to resolve this possibility.

Equilibrating Refractory Particulates. If we do not expect significant flux changes from the host star and we assume that the planet is tidally locked, then thermal variability that we observe before/after eclipse must be occurring on the face of 55 Cnc e. Due to its close-in orbit and proximity to its host star, the planet most likely has a molten surface and may experience intense tidal forces, both of which could drive volatile loss through either surface evaporation and/or volcanic activity (Schaefer & Fegley 2009). Volcanic eruptions or a volatile atmosphere could result in particulates either lofted or condensed high above the surface (Mahapatra et al. 2017). D16 examined this as a source of thermal variability, with volcanic plumes raising regions of the planet's photosphere to higher, cooler regions in 2012, which led to smaller eclipse depths (i.e., less thermal emission) that year.

Here, we explore a similar but more general concept. For each of our models (except the flat model, for which this does not apply because the depth is constant), we assume that our measured value of T_b,max represents the temperature of the radiating planetary surface, and we check to see if the planet's equilibrium temperature (due to stellar insolation alone) is consistent with this measurement. If the planet has a molten surface or is volcanically active, it may outgas material that cools at a height above the planet, blocking out the hot surface and leading to lower temperatures. For each model, we check to see if T_b,min can be achieved by this obscuring volcanic ejecta.

To start, we determined the equilibrium temperature of 55 Cnc e due to stellar insolation alone over a range of planetary albedos. We did this both for an efficient cooling case, in which the planet is able to radiate heat away over its whole surface, and an inefficient case, in which it can only radiate over half the surface.

Next, for each model, we calculated the equilibrium temperature of refractory particulates at a height of 100 km above the planet's surface due to both solar and planetary radiation over a range of grain albedos. In each case, we assumed that the planet is radiating at its T_b,max temperature. Our assumption of the height of the grain layer above the surface is about twice the maximum theoretical height of volcanic plumes on earth (Wilson et al. 1978), which we feel is not unrealistic for a planet of ∼2 R_⊕ that may experience intense tidal heating. In practice, however, we found that our results for this calculation are not strongly dependent on the location of the layer of grains. The results of these calculations are shown in Figure 12. The top panel shows the equilibrium temperature for efficient and inefficient cooling models compared to the T_b,max of the model, while the bottom shows the temperature of grains compared to T_b,min.

Zoom In Zoom Out Reset image size

While results vary for the different models, they all require a dark planetary surface (albedo <0.5) and grains that are reflective (albedo >0.4) at visible wavelengths where the star releases most of its energy in order to be consistent with the 2σ ranges of T_b,max and T_b,min based on planetary equilibrium temperature alone. They also require the inefficient cooling case to agree with T_b,max, which may be indicative of a lack of an atmosphere for efficient heat redistribution.

However, it is possible that the planet's thermal radiation budget is not limited to incoming flux from its star. If the four other planets in the system (Fischer et al. 2008) help maintain even a slightly eccentric orbit, 55 Cnc e may dissipate a significant amount of radiation in the form of tidal heating. Demory et al. (2016a) examined tidal dissipation as a possible explanation for a 3100 K region in their longitudinal temperature map of 55 Cnc e, using the Mercury-T code of Bolmont et al. (2015). They found that tidal heat fluxes ranging from 10⁻³ to 10⁶ W m⁻² are achievable for the planet, corresponding to temperature increases of a few Kelvin to ∼2000 K. While the amount of tidal heating is largely unconstrained by current observations, it could represent a significant fraction of the planet's thermal radiation budget. Coupled with a low planetary albedo, tidal heating could significantly increase the surface temperature of the planet, which in turn could allow T_b,max to be consistent with the efficient cooling model.

Budaj et al. (2015) tabulated the opacities and equilibrium temperatures of several species of dust grains that may be present in exoplanet atmospheres over a wavelength range from 0.2 to 500 μm and grain radii from 0.01 to 100 μm. Extracting their opacity results at 4.5 μm from their tables,⁷ we calculated the single scattering albedo (Equation (23) from their paper) for all grain sizes from all species they examined. We used their temperature tables for a star of T_eff = 5000 K and solid angle of 0.01 and 0.03 steradians (55 Cnc subtends a solid angle of ∼0.02 sr in 55 Cnc e's sky).

Cross-referencing these results, we identified the composition and size of grains that have the requisite albedos and equilibrium temperatures to match our brightness temperatures. These include alumina, olivine, and pyroxene, with grain radii ranging from around 0.01 to 0.6 μm. If the planet regularly condenses or outgasses these materials in an opaque cloud layer, its composition could achieve the appropriate albedo and equilibrium temperature to match T_b,min from our models.

Mahapatra et al. (2017) applied a kinetic, non-equilibrium cloud formation model to the atmosphere of 55 Cnc e using several different abundance cases. Using D16's derived T–P profile, they found that 55 Cnc e can support a highly opaque cloud layer consisting of mainly Mg silicates and Si oxides. The local gas temperature of this cloud layer (see panel 1 of Figure 10 in their paper) ranges from about 850 to 950 K as a function of pressure, which is within the 1σ errors of T_b,min for our year-to-year, slope, and sine models.

As a caveat, our model would require significant portions of the planetary surface to be blocked by a layer of cooler refractory grains, depending on the local structure of planetary surface temperature and the temperature of occulting material. Determining whether or not the maintenance of such a large layer of condensed or ejected material is feasible is beyond the scope of this simple calculation.