Confirming Variability in the Secondary Eclipse Depth of the Super-Earth 55 Cancri e
DDraft version April 16, 2018
Preprint typeset using L A TEX style emulateapj v. 01/23/15
CONFIRMING VARIABILITY IN THE SECONDARYECLIPSE DEPTH OF THE SUPER-EARTH 55 CANCRI E
P. Tamburo , A. Mandell , D. Deming , and E. Garhart Boston University, 1 Silber Way, Boston, MA 02215 NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, MD 20771 University of Maryland College Park, College Park, MD 20742 NASA Astrobiology Institute’s Virtual Planetary Laboratory Arizona State University, Tempe, AZ 85281
Draft version April 16, 2018
ABSTRACTWe present a reanalysis of five transit and eight eclipse observations of the ultra-short period super-Earth 55 Cancri e observed using the Spitzer Space Telescope during 2011-2013. We use pixel-leveldecorrelation to derive accurate transit and eclipse depths from the
Spitzer data, and we performan extensive error analysis. We focus on determining possible variability in the eclipse data, as wasreported in Demory et al. (2016a). From the transit data, we determine updated orbital parameters,yielding T0 = 2455733.0037 ± ± ± R p = 1.89 ± R ⊕ . Our transit results are consistent with a constant depth, and we conclude that theyare not variable. We find a significant amount of variability between the eight eclipse observations,and confirm agreement with Demory et al. (2016a) through a correlation analysis. We convert theeclipse measurements to brightness temperatures, and generate and discuss several heuristic modelsthat explain the evolution of the planet’s eclipse depth versus time. The eclipses are best modeled by ayear-to-year variability model, but variability on shorter timescales cannot be ruled out. The derivedrange of brightness temperatures can be achieved by a dark planet with inefficient heat redistributionintermittently covered over a large fraction of the sub-stellar hemisphere by reflective grains, possiblyindicating volcanic activity or cloud variability. This time-variable system should be observable withfuture space missions, both planned ( JWST ) and proposed (i.e.
ARIEL ). Keywords: planets and satellites: general — planets and satellites: individual: 55 Cancri e — tech-niques : photometric — occultations INTRODUCTION
Planets with extremely short orbital periods (P < e , a transiting super-Earth orbiting aSun-like star about 13 parsecs from Earth. Detected byMcArthur et al. (2004) and confirmed by Fischer et al.(2008), the planet was originally thought to have a 2.8day period and a minimum mass of 14 M ⊕ . Dawson &Fabrycky (2010), however, proposed that the detectedperiod and mass were due to aliasing in the radial veloc-ity data, and that the planet’s true period and mass were0.74 days and 8 M ⊕ . Winn et al. (2011) and Demory etal. (2011) made photometric observations that confirmedthe proposed values.A rapid orbit around a bright star that permits rel- [email protected] atively high signal-to-noise observations (von Braun etal. 2011) makes 55 Cnc e a favorable target for perform-ing frequent transit and eclipse observations that allowfor further characterization of the planet. Our desire tounderstand the nature of super-Earths has resulted in aplethora of data for the system (e.g. Dragomir et al. 2014;Demory et al. 2016a), and various studies have attemptedto characterize the planet’s structure and composition.Demory et al. (2016b), for example, created a day-nighttemperature map of the planet using phase-curve data,while Lopez (2017) found that the planet must have awater envelope of 8 ±
3% of the total planet mass to beconsistent with the mass and radius.One of the most notable results is that of Demory et al.(2016a) (hereafter D16), who found that the secondaryeclipse depth of the planet increased significantly overthe course of eight observations in the 4.5 µ m channelof the Spitzer Space Telescope (Werner et al. 2004) be-tween 2012 and 2013.
Spitzer has yielded high preci-sion photometry that has allowed for the characterizationof many diverse exoplanetary systems (e.g. Gillon 2017;Knutson et al. 2012). The increase in eclipse depth is in-dicative of an increase of the planet’s brightness temper-ature ( T B ) at 4.5 µ m from 2012 to 2013, and D16 foundthat T B increased from about 1600K to 2600K over thattime span. D16 posited that the variation could possiblybe accounted for by volcanism on the planet. Volcanicplumes could in principle raise the photosphere at 4.5 µ mhigher in the atmosphere where the temperature is lower, a r X i v : . [ a s t r o - ph . E P ] A p r Tamburo et al. leading to lower thermal emission when the plumes areactive (Demory et al. 2016a).Genuine detection of eclipse variability would thereforebe an extraordinary result, as it could indicate a volcanicprocess occurring on a body outside of our solar system.Because extraordinary results require extraordinary ev-idence, we present here a reassessment of the possiblevariability in eclipse and transit depths for this system.We not only use a new and powerful method to correctfor the instrumental signatures in the data, but we alsodevelop a new method for error analysis that providesan independent check on our quoted errors. Moreover,we evaluate a wider range of heuristic models for eclipsevariability than did D16. In Section 2, we describe theobservations of 55 Cnc e that were analyzed, as well asour data analysis procedures. In Section 3, we discussthe results for both the transit and eclipse observations.Section 4 discusses the implications of our analysis, andhow the implications of our results relate to the thermalvariability of 55 Cnc e . DATA ANALYSIS
Observations
We examined eight eclipses and six transits of55 Cnc e taken in the 4.5 µ m channel of the Spitzer Space Telescope’s
IRAC instrument (Fazioet al. 2004) using subarray mode. These data areidentical to those analyzed by D16, and are pub-licly available through the Spitzer Heritage Archive(http://sha.ipac.caltech.edu/applications/Spitzer/SHA/).We list the date, type, program ID, and AOR of eachobservation in Table 1.All observations were taken with 0.02-s exposures us-ing subarray mode. Observations taken after the 2012January 18 eclipse were all taken with PCRS peak-up,which is designed to limit stellar motion on the detec-tor and hence reduce the intra-pixel effect (Ingalls et al.2012). Table 1 gives the date, observation type, andAOR number for all the data used.Intra-pixel sensitivity fluctuations in the 4.5 µ m chan-nel can cause flux variations of about 1%, depending onhow the stellar PRF moves over the course of observa-tions (e.g. Charbonneau et al. 2005; Mighell et al. 2008).This effect is two orders of magnitude larger than theexpected eclipse depths, and has to be accurately re-moved in order to draw scientific conclusions from theobservations. Our process for removing the intrapixelfluctuations is detailed in Section 2.3. Photometry
The data comprise cubes of 64 frames, each having32x32 pixels. The star is very bright and our exposuresvery short, so background is negligible. Nevertheless, wedetermine the average background value by first fittinga Gaussian to a histogram of pixel intensities over theentire 32x32 pixel image to determine the centroid value.This value was then subtracted from our images.Because knowing the position of the star in the im-ages is necessary in order to place a numerical apertureand perform accurate photometry, we found the locationof the star in the images using two different methods.First, we determined the centroid of the stellar image us-ing a 2-D Gaussian fitting procedure in IDL. Second, we performed a center-of-light (COL) calculation, that sep-arately defined the centroid of the stellar image in eachcoordinate. We collapsed the image in the Y dimension,then calculated the centroid in X by intensity weightingthe X-coordinates, and vice-versa to calculate the imagecentroid in Y. Both of these fits were performed using thefull 32x32 pixel images. Having determined the center ofthe star, we performed aperture photometry for both the2-D Gaussian and COL methods.For both centering methods, we summed the stellarflux in 11 numerical circular apertures of various radii,ranging from 1.6 to 3.5 pixels. We used both constantand time-variable radii, the latter being determined usingthe “noise pixel” implementation of Lewis et al. (2013),to which we added constant increments in radius rangingfrom 0.1 to 2 pixels. The two approaches for identifyingthe location of the star along with the option for fixed andtime-variable apertures gave us four different photometryfiles that can be used in our fitting code.In practice, we found that it is advantageous to limitourselves to using only the photometry files that useGaussian centroiding and fixed apertures. In Figure 1,we show the performance of various fitting statistics forthe four photometry methods versus the amount that thedata were binned for one of the eclipses (data were binnedover a range of timescales in our fitting procedure to bestremove the intrapixel effect, as Deming et al. (2015) dis-cuss). The standard deviation of normalized residuals, orSDNR, is a measurement of the scatter in the unbinneddata after the fit has been calculated and applied, andthe Gaussian-fixed method achieved a smaller value overall bin sizes compared to the other three methods. Vari-able aperture photometry is advantageous when there issignificant variation in the shape of the apparent stellarimage, such as can be produced by pointing jitter duringan exposure. However, when the exposure time is veryshort (as it is here), pointing jitter within an exposureis minimal. In that case, variable aperture photometrycan produce inferior results because the determinationof the aperture radius is always affected by the photonnoise in the image, resulting in spurious and unnecessaryvariations in the aperture radius.All four photometry methods performed comparablywell with the Allan variance slope, which is shown inthe middle panel of Figure 1 and covered in more de-tail in Section 2.3.3. Lastly, the reduced- χ of fits (bot-tom panel of Figure 1) made using the Gaussian-fixedphotometry was better at all bin sizes than the otherapproaches. Similar results for other observations moti-vated us to restrict ourselves to only examining the pho-tometry made using fixed apertures and Gaussian cen-troiding.To clean up the photometry, we discarded pixels thatwere discrepant by more than 6 σ of the standard devi-ation of the intensity of 9 neighboring pixels. If a givenpixel was outside of this range, we replaced it with themedian value of the 9-pixel neighborhood. We investi-gated the effect of our treatment of discrepant pixels inmore depth. Strictly speaking, discrepant pixels shouldbe zero-weighted in the photometry, not replaced by amedian value (Horne 1986). However, zero-weightingwould be difficult to implement, given the structure ofour photometry code, so we investigated the magnitudeof possible bias introduced by median-filtering. We be- ariability in 55 Cnc e ’s Eclipse Depth Table 1
Date, type, program ID, and AOR number of observations analyzed in this work. We also list the percentage of frames discarded (
F r disc )and pixels corrected (
P ix corr ) for each observation, per our discussion in Section 2.2. σ phot are the photon-limited errors of the unbinneddata, while SDNRσ phot is the ratio of the standard deviation of unbinned residuals (photometry-fit) to the photon-limited errors. On average,the scatter of our residuals comes within 10% of the photon noise limit.Date Type Program ID AOR
F r disc (%)
P ix corr (%) σ phot (ppm) SDNRσ phot
Figure 1.
Standard deviation of normalized residuals, Allan vari-ance slope, and reduced χ versus bin size for our four differentphotometry methods for the 2012 January 31 eclipse. Similar re-sults for other observations motivated us to limit ourselves to us-ing the photometry with Gaussian centroiding and fixed aperturesexclusively. The different strands of similarly-colored lines corre-spond to different aperture sizes used to perform the photometry. gan with synthetic pixel data whose average value andscatter are the same as the real data. The frequency atwhich discrepant pixels occur is tracked by our photom-etry code, so we know on average how many pixels arereplaced in the real data (see Table 1). We chose thatnumber of pixels at random, and calculate the effect ofreplacing them with a median filter in the synthetic data.Doing that 10,000 times, we calculated the average effectover the duration of the eclipse. We found that the effectis less than 0.5 ppm, and therefore of negligible signifi-cance to our analysis. We conclude that median-filteringof discrepant pixels is a valid method to clean the dataprior to deriving eclipse depths.We also discarded images that exhibited a significantamount of jitter in image position using a two-pass tech-nique. The first pass made an absolute cut of images thatexhibited image motion greater than 0.3 pixels, whicheliminated large excursions. The second was a “soft” cutthat filtered images that showed jitter in disagreementwith the standard deviation of the residuals of image po-sition minus the median image position by more than 6 σ .The percentage of discarded frames and corrected pixelsare shown in Table 1.Because the current version of PLD is not designed tohandle large amounts of image motion, we performed fitsover a limited range of orbital phase around the eclipse.The size of this range is defined such that the amount ofout-of-eclipse data is approximately equal to the amountof in-eclipse data. The first-order implementation of PLDused in this paper is valid for less than ∼ Fitting for Eclipse Depths
Pixel-Level Decorrelation
To remove systematic intra-pixel sensitivity fluctua-tions, we utilized the PLD framework presented by Dem-ing et al. (2015). Other methods currently used to cor-rect for the intra-pixel effect depend on relating intensityfluctuations to the position of the stellar image on the de-tector (e.g. BLISS mapping [Stevenson et al. 2012], usedby D16). As Deming et al. (2015) discuss, the position of
Tamburo et al. the stellar image is a secondary data product, ultimatelyderived from the intensities of individual pixels. PLDutilizes the primary data, the intensities of individualpixels, to correct intra-pixel sensitivity variations.PLD is capable of removing red noise that can frus-trate other methods. Ingalls et al. (2016) comparedseven different Spitzer decorrelation methods, and foundthat PLD (as well as BLISS mapping used by D16) wasamong the three methods that came closest to finding thetrue eclipse depth when applied to synthetic Spitzer data(however, it should be noted that six of the seven meth-ods compared in Ingalls et al. (2016) produced resultsthat were indistinguishable within the expected measure-ment error when applied to real and synthetic data). ThePLD technique has been used successfullly in a numberof recent
Spitzer analyses (Kammer et al. 2015; Buhler etal. 2016; Fischer et al. 2016; Wong et al. 2016; Dittmannet al. 2017; Kilpatrick et al. 2017). PLD is also the con-ceptual basis for the highly successful EVEREST code(Luger et al. 2016) that achieves high precision in decor-relating data from the K2 mission.Other known systematic effects that PLD must cor-rect for include a temporal “ramp” effect, in which thereported intensities of detector pixels change rapidly un-til they reach a threshold value, and temporal baselineeffects, in which the intensity gradually increases or de-creases over the course of observations. To account forthe transient temporal ramp, we simply avoided datanear the start of observations (see Section 2.2). To re-move the second effect, we fit visit-long functions to thephotometry (linear, quadratic, or exponential) and sub-tracted them from the photometry. Through testing withthe Bayesian Information Criterion (covered in more de-tail in Section 4.2.2), we found that the best fits areachieved through use of visit-long linear baselines. Forthese observations, we found that not fitting a baseline atall typically results in failure to find a regression solution.As a result, we are confident that linear baselines are thebest option for removing temporal baseline effects.Previous applications of the code used 9 (Deming et al.2015) or 12 (Garhart et al. 2018) pixels, in accordancewith the number of pixels that captured a significant por-tion of the stellar flux in those works. Here, we haveupdated PLD to work with anywhere from 1-25 pixelsin a 5x5 grid surrounding the star. While in principle,the code can be can be used with any set of pixels inthis grid, we expect the starlight to be distributed sym-metrically about the central pixel upon which the star isplaced in the image. This expectation, coupled with thefact that edge pixels in the 5x5 box captured 15% of thetotal flux on average for these observations (see Figure2), motivated us to use a 5x5 pixel grid to decorrelatethe photometry.Figure 3 shows a typical application of the regressionfit on a 2011 transit, for the best binning and aperturesize (see Sec. 2.3.3). In the figure, we show the processof going from raw photometry, to binned fit, to transitmodel with systematics removed.
Determining Central Phase
In the original version of PLD, the central phase of theeclipse was determined through fitting an eclipse modelto unbinned data over a range of trial values, selecting theone that minimized χ . We found that this method was Figure 2.
The fraction of normalized flux captured by each pixelin the 5x5 grid averaged over all observations. While the central3x3 grid captures about the majority of the flux, we elect to usethe entire grid to decorrelate the photometry because edge pixelscapture about 15% of the flux in the images, on average. unreliable for these observations of 55 Cnc e , as the codecould have difficulty in finding the small transit/eclipseevents in the unbinned data. As a workaround, we setan initial guess for the central phase, based on litera-ture values for the ephemeris (Endl et al. 2012), and ranan instance of the MCMC on binned data to refine thecentral phase. This value was then used in subsequentinstances of the code.For several of the eclipses, we found that even thisworkaround could not find a definite central phase, dueto the shallow depths of the eclipses. We therefore electedto generate an ephemeris using the transit events, whichare several times deeper than the expected eclipse depths.We used this ephemeris to calculate the expected eclipsetimes/central phases, assuming a zero-eccentricity or-bit. Previous works (Nelson et al. 2014; Demory et al.2016a) found eccentricities consistent with zero, despitethe presence of four other planets in the system. We alsocalculated the tidal circularization timescale for 55 Cnc e , using τ e = 463 Q ( a GM ) . ( mM )( aR p ) (1)(Goldreich & Soter 1966), where Q is the planet’s spe-cific dissipation function. Assuming a Q = 100 (Glad-man et al. 1996), we find a tidal cirularization time ofabout 14,000 years, which justifies the assumption of azero-eccentricity orbit. Selecting Aperture and Bin Sizes
With the central phase provisionally determined, wevaried the photometric aperture size, and binned thephotometry over a range of time scales. Central to thePLD method is that it usually fits to data that are binnedin time, as discussed in Sec. 3.1 of Deming et al. (2015).Kammer et al. (2015) also discuss the rationale for bin- ariability in 55 Cnc e ’s Eclipse Depth Figure 3.
Application of PLD fit to photometry for the 2011June 20 transit.
Top:
Raw photometry, showing 239,443 unbinnedpoints.
Middle:
Binned photometry (grey and black circles) andoverlaid best fit from PLD regression (red). The binning of the greyand black points are about 1 and 10.5 minutes, respectively, forclarity. This fit is reapplied to the unbinned photometry to examinethe standard deviation of residuals on all time scales representedby the event.
Bottom:
Binned photometry with intrapixel effectremoved, with the best fit transit curve overlaid. Both grey andblack points have the same binning as the middle panel. ning in time. Briefly, if the scatter in the data is dom-inated by stochastic noise, it will be impossible to de-tect correlated noise and match it to our pixel basisvectors. By binning the data, we can diminish short-cadence white noise and expose long-cadence correlatednoise. Fitting to binned data involves a loss of informa-tion (e.g., due to high frequency image jitter). However,we find that fitting to binned data produces a bettermatch between the time scale of the fit and the dominanttime scale of the intra-pixel fluctuations. Other authors(e.g. Deming et al. 2015; Buhler et al. 2016) have alsoused 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 theintra-pixel effect and the eclipse of the planet. We used3 constant aperture radii (see Section 2.3.6), and 112 binsizes ranging from 2 to 2048 points-per-bin. For a giventrial of bin size and aperture radius, we found the bestfit by linear regression (solving Eq. 4 from Deming et al.2015), holding the central phase of the eclipse fixed atthe provisional value. We then applied this trial fit to the unbinned data, and we investigated the noise propertiesof those residuals as a function of time scale. We denote the standard deviation of the unbinned residuals as σ (1).We binned the residuals over 2, 4, 8, 16, etc. points, untilthe effect of binning would reduce the number of residualsto ≤
16 points. We calculated the standard deviation ofthe 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) shouldvary inversely with the square root of the bin size (a slopeof -0.5 in log space). We list the SDNR values of ourfits compared to the photon noise in Table 1. The rela-tion between standard deviation and bin size is knownas an Allan variance relation (Allan 1966). The origi-nal PLD method Deming et al. (2015) adopted the trialfit (i.e., bin size and aperture radius) that produced theminimum chi-squared in the Allan variance relation. Wehave found it desirable to slightly modify that originalcriterion, as we now explain.In choosing the best fit from different aperture radiiand (especially) bin sizes, PLD effectively makes a com-promise between mitigating the short-term scatter, rep-resented by the standard deviations of unbinned normal-ized residuals (SDNR), and the longer-term red noise,represented by the standard deviations of the residualsbinned to longer time scales. The standard deviation ofthe residuals itself has an uncertainty (i.e. the error onthe error). That uncertainty increases with bin size be-cause larger bins sample the error distribution with fewerpoints. Hence the chi-squared of the Allan variance rela-tion is dominated by the smallest bins sizes. For that rea-son, we found that the original criterion used by Deminget al. (2015) (minimum chi-squared in the Allan variancerelation) is virtually equivalent to choosing the solutionwith minimum unbinned scatter (i.e, minimizing SDNR).We therefore found it advantageous to modify the bestfit criterion to choose the solution that has the minimumraw scatter in the Allan variance relation ( σ AVR), afterforcing the relation to pass through σ (1). That is similarto the original criterion from Deming et al. (2015), butwith greater weighting for the larger bin sizes. Note thatthis fitting criterion is used (successfully) in Section 2.3.6wherein we derived null-test eclipse depths to check ourerrors.In order to quantify the improvement of using our newversion of the PLD code, we tested our new version versusthe original code described by Deming et al. (2015). Weused three representative eclipses (2012 Jan. 18, 2013Jun. 15, and 2013 Jun. 29) and derived eclipse depthsusing the original code and our new code applied to thesame photometry files, using the same degree of databinning. We also constrained the central phase of theeclipse to be the same in each code. We found that thenew code reduces the single-point photometric scatter inthe residuals (data minus fit) from 1.27 times the photonnoise to 1.11 times the photon noise. Given the largenumber of data points (239,000), that degree of reductionin the photometric scatter easily justifies (i.e, via a BICanalysis) the inclusion of the additional basis pixels in thePLD solutions. The average error on the eclipse depth,from the MCMC posterior distributions, is reduced from50 ppm produced by the original code to 39 ppm fromour new code. The average eclipse depth from the twocodes was the same to within 1.2 times the error of themean. Tamburo et al.
Stability of Eclipse Depth Fits
In testing these observations, we found that many com-binations of bin and aperture size produce comparablygood σ AVR scores for a given eclipse. The eclipse modelsassociated with those combinations, however, could haveeclipse depths that differed by 10s of parts-per-millionbetween each other (consistent with their errors), differ-ences that could be significant when searching for po-tential variability between such small events. The dis-tribution in eclipse depth for similar σ AVR scores is aresult of the fact that σ AVR is not a rigorous statisti-cal criterion. Instead, it’s an attempt to get a “broad-bandwidth” solution by trying to minimize both red andwhite noise through binning. Because of this, we selectedour best-fit eclipse depth for each observation by exam-ining a distribution of solutions, chosen by their σ AVRscores.After running the regression for all bin and aperturesizes, we fit a Gaussian to the distribution of σ AVR val-ues, which we refit after removing 3 σ outliers. We se-lected the center of this Gaussian as a cutoff value, andwe then focus on solutions with σ AVR scores better thanthe cutoff. There are typically ∼
100 solutions left af-ter removing the bad solutions, and the eclipse depths of these best solutions were approximately Gaussian-distributed for each eclipse. To determine the depth thatwe quote as our result, we fit a Gaussian to this distribu-tion, selecting the combination of aperture and bin sizethat produced the eclipse depth that is closest to thecentral value. This approach allowed us to find a solu-tion that is the most stable in the sense of representingthe best fits found by the code, and prevented us fromreporting an outlier as the best eclipse model.
MCMC
After finding the best combination of bin and aperturesize, we held the binning and aperture constant, and re-fined the fit using a Markov Chain Monte Carlo proce-dure (Ford 2005). While we report the depth associatedwith the regression solution, the MCMC allowed us toexplore a high-dimensional parameter space (our modelsuse ∼
30 parameters when including the pixel coefficients)which accounted for possible correlations between fittingparameters and defines the error on the eclipse depth viaa Gaussian fit to the posterior distribution.The MCMC utilized the Metropolis-Hastings algo-rithm with Gibbs sampling. At each step, we adjustedthe eclipse depth, select orbital parameters (for transitfits), temporal baseline and pixel coefficients. Each chainconsists of 500,000 steps, and we ran multiple chains so asto verify convergence through the Gelman-Rubin statis-tic, R (Gelman & Rubin 1992).Comparing our error estimates from the MCMC withthose found by the best regression solution chosen in Sec-tion 2.3.3, we found that the MCMC errors are on aver-age 13% higher than the regression errors for the eclipses.We also tested initializing the MCMC with differentsolutions than the one chosen by our method of selectinga representative fit from the distribution of best σ AVRscores. These fits are typically within 2 ppm of the “best”solution, but use different aperture and bin sizes. We findthat the different initializing solutions did not producesignificantly different errors from the solution chosen by our method. We therefore conclude that the errors arestable in the sense that they are not sensitive to the de-tails of the fit.
Independent Check on Errors
As stated above, our errors on the eclipse depths aredetermined using the MCMC simulations described inthe previous section. However, deriving realistic errors iscritical to the inference that the eclipse is variable, anda key feature of our analysis is an independent check onour derived errors.To make this check, we identified eight regions in the2013 eclipse data (two for each of the four eclipses) thatwere independent of each other and do not overlap withthe actual eclipse. In these regions, we expect to deriveeclipse depths that are consistent with zero within ourerror bars, because we know that they do not containthe actual eclipse. A critical test of the errors is that thederived “eclipse” depths for these regions should scatteraround zero with a dispersion that is consistent with theirerrors.We “fit” these non-eclipse regions following themethodology of Sections 2.3.1-2.3.5, locking the centralphase of the model in each case. From this test, wefound that if we restrict ourselves to using large photo-metric apertures (with radii of 3.1, 3.3, and 3.5 pixels, thesame used in our eclipse analysis), we indeed find eclipsedepths that vary around zero to a degree consistent withtheir derived errors (calculated through the standard de-viation of in-eclipse and out-of-eclipse residuals, added inquadrature). The restriction to the larger aperture sizesmakes sense given our choice of using all basis pixels in a5x5 pixel grid surrounding the star: these apertures arewell matched to the area of this basis pixel grid. Conse-quently, we limit ourselves to using these three apertureswhen fitting the actual transits and eclipses.The results of the fits to the non-event regions areshown in Figure 4, which plots the distribution of eclipsedepth divided by its derived error. If our errors are real-istic, that distribution should be consistent with a Gaus-sian having standard deviation of unity. The distributiononly contains eight samples, because we are limited bythe amount of independent data regions that are avail-able. Nevertheless, the distribution is approximately cen-tered on zero with a dispersion consistent with unity. Our1 σ error bars, which average to 34 ppm for these data,encompass all of the non-event results except two. Thefirst of these is within 1.3 σ of zero, which we expect witha frequency of once every five measurements. The secondis within 1.7 σ of zero, which is expected once every 11measurements.Moreover, performing an Anderson-Darling test fornormality (Stephens 1974), we find a p-value of 0.43,implying that our results for the non-event regions areconsistent with a Gaussian having a standard deviationof unity. We conclude that our fitting procedure pro-duces realistic errors, suitable for drawing inferences asto the variability of the eclipse. RESULTS
Transit Analysis
Early in our analysis, we found that the eclipse obser-vations poorly constrain the planet’s orbital parameters ariability in 55 Cnc e ’s Eclipse Depth Figure 4.
Results of testing the realism of our eclipse depth er-rors. The abscissa plots the distribution of each “eclipse” depthdivided by its error, with the histogram on the ordinate. If theerrors are realistic, then the histogram should be consistent with aGaussian distribution having a standard deviation of unity, over-plotted on the histogram. In spite of our limited number of samples(8), this test shows that our errors are realistically determined. which can directly affect the retrieved eclipse model. Asa result, we hold these parameters constant in our eclipsefits. However, we first needed the best possible estimatesof these values, and because using different sets of pub-lished orbital parameters (e.g. Demory et al. 2016a,b;Baluev 2015; Nelson et al. 2014) gave slightly differentvalues for inclination and central eclipse times, we electedto determine our own. The expected transit depths aresignificantly larger than the eclipse depths (by a factorof ∼ R p /R ∗ , or-bital inclination i , quadratic limb darkening coefficients u and u , the multiplicative coefficient of our linear Figure 5.
Raw photometry for the 2013-07-08 and 2013-07-11transits. The sudden shifts in the photometry values are due toshifts in the stellar image on the detector, and limit PLD’s abilityto fit the data accurately.
Table 2
Transit depths found in this work and D16. Our results share acorrelation coefficient of 0.67 with D16’s values.Date This Work (ppm) Demory et al. 2016 (ppm)2011 Jan 6 380 ±
36 484 ± ±
35 287 ± ±
38 325 ± ±
36 365 ± ±
77 406 ± ±
18 360 ± ramp, and mid-transit time. We adopted the publishedvalue of i from Demory et al. (2016b) as our startingvalue for inclination, and placed a Gaussian prior onthis value based their computed error. We also prioredthe limb darkening coefficients u and u about their ex-pected values from the tables of Claret & Bloemen (2011)for the published stellar parameters (von Braun et al.2011). The convergence of these chains was determinedby checking that the Gelman-Rubin statistic R (Gelman& Rubin 1992) for the various parameters was less than1.1.To determine the uncertainty on our regression tran-sit depths, we performed a Gaussian fit to the combinedMCMC posterior distribution of R p /R ∗ for each tran-sit. The resulting depths/errors can be found in Table 2,along with D16’s transit results. Models of each transitcan be found in Figure 3.1. We calculated the correla-tion coefficient between our transit depths and D16’s andfound a value of ∼ Tamburo et al.
Figure 6.
Best fit MCMC transit models. 2011 observations shown on top, 2013 observations shown on bottom.
Table 3
Self-consistent parameter values as determined by performingMCMC procedures on transit observations. R p is calculated usingthe weighted average transit depth (see Table 2), and assumes thepublished stellar radius from von Braun et al. (2011).Parameter ValueP (days) 0.7365454 ± T - 2,400,000 (HJD) 55733.0037 ± R p ( R ⊕ ) 1.89 ± erage of our transit depths, and found that it agrees wellwith D16’s weighted average.We tested the sensitivity of the transit depth errorsmeasured from the MCMC to the priors by doubling theprior widths on i , u , and u , and re-running three chainsfor each transit. We found that the errors measured fromthese runs differ by at most 2 ppm from those inferredfrom the original runs. We therefore conclude that theerrors we report on transit depth are insensitive to theimposed priors.We adopted a similar procedure for determining theplanet’s orbital parameters. We fit a Gaussian to thecombined posterior distributions for central phase foreach transit and use the value of each to perform a linearleast squares regression to determine P and T0. For in-clination, we combined the posterior distributions of allthe transit observations and performed a Gaussian fit.The values we determined for these parameters can befound in Table 3. We note that our values for P and i arewithin the 1 σ errors of the values published in Demoryet al. (2016b).Our value for T0, however, is a little over 1 σ awayfrom the value reported in Demory et al. (2016b). Thisis because our code finds that the first two transit eventsoccur 16 and 12 minutes before the times expected fromthat paper’s ephemeris, respectively. We tested initial-izing the MCMC with the value expected from theirephemeris, but the code still converged to the earlierevent times. We also checked our agreement for thesetwo observations with the ephemeris from Winn et al.(2011), made using observations with MOST. Those val- ues predict transits within 2.2 ± ± ∼
11 minutes and ∼ . ± . R ⊕ , a value in agreementwith Demory et al. (2016a,b). Eclipse Analysis
For our eclipse analysis, we followed the same proce-dure as the transits to determine our regression solutions.In the MCMC, we allowed the pixel coefficients, lin-ear slope coefficient, vertical offset, and eclipse depth tovary. Central phase, period, and inclination were lockedto their expected values from our transit analysis (3.1),because these elements were not strongly constrained bythe weak eclipses. a/R ∗ was locked to the value reportedby Demory et al. (2011), and we assumed an eccentricityof zero for the planet’s orbit. In Figure 7, we show cor-ner plots of the non-pixel coefficients for one run of theMCMC for all eight eclipses. The lightcurve parametersshow no correlation with individual pixel values, con-firming that no degeneracy exists between the lightcurvesolutions and the position of any specific pixel.The eclipse depths that we found are given in Table4, along with the reduced χ of the regression fits thatinitialized the MCMC. Lightcurve models for these fitsare shown in Figure 8. DISCUSSION
Transit Variability
D16 found a 25% variation in their transit depths incomparison to the mean value between 2011 and 2013 ariability in 55 Cnc e ’s Eclipse Depth Figure 7.
MCMC correlation plots for one 2012 and one 2013 eclipse. There is no correlation between any of the lightcurve parametersand any individual pixel position values.
Figure 8.
Eclipse light curve models. 2012 eclipses shown on top, 2013 eclipses on the bottom. Tamburo et al.
Table 4
Eclipse depths found and reduced χ of regression fits foundthrough PLD analysis. The depth given is our best regressionresult, while the error is derived from Gaussian fits to theposterior distribution of eclipse depths from the MCMC.Date Depth (ppm) χ ν ±
44 1.122012 Jan 21 -23 ±
37 1.282012 Jan 23 27 ±
35 1.242012 Jan 31 53 ±
45 1.262013 Jun 15 175 ±
38 1.132013 Jun 18 171 ±
35 1.192013 Jun 29 84 ±
39 1.202013 Jul 15 113 ±
38 1.25
Figure 9.
Our transit depth results versus time, with 2011 shownon the left and 2013 on the right. The 2 σ range of a flat line fit tothe depths is shown in red. at the 1 σ level, which they considered to be a marginaldetection of variability. With PLD’s improvement overBLISS mapping’s uncertainties for the transit observa-tions (a 37% reduction, on average), we examined ourown results to see if we find stronger support for transitdepth variability.We began by calculating the correlation coefficient ofD16’s transit depths and our own, finding a value of 0.67,indicative of agreement between the two sets of results.We then calculated a weighted average of our transitdepths, giving 336 ±
18 ppm. This calculation representsthe null hypothesis, that the transits can be representedby a single constant depth. We show our transit depthsversus time with the 2 σ range of this flat line fit over-plotted in Figure 9. The non-variable model gives a χ of6 . n = 4 degrees of freedom,gives a p-value of about 0 .
2. Our results for the first fivetransits, therefore, are consistent with a constant depth,as we would expect a higher χ value around 20% ofthe time due to random chance alone. On this basis,we conclude that our results do not provide evidence fordifferent transit depths for 55 Cnc e in 2011 comparedwith 2013. Therefore, both the results of D16 and thiswork find insufficient evidence to support a conclusion ofyear-to-year variability in 55 Cnc e ’s transit depth.However, the separate linear trends in the transit depthwithin each year found by D16 are still present in our re- sults, though the variation falls within our 3 σ uncertaintyenvelope. Due to the short baseline of the observationsfor each year, we cannot rule out intra-annual variability,and future higher-cadence and/or longer-baseline transitmeasurements over > Eclipse Variability
One way to determine the reality of statistical fluctua-tions in 55 Cnc e ’s eclipse depth in the presence of pos-sible analysis errors is to determine whether both D16’sapplication of BLISS mapping and our own use of PLDproduce correlated results when comparing eclipse-to-eclipse. We calculated the correlation coefficient betweenthe sets of eclipse depths found through these two meth-ods, and find a value of 0.90, which tells us that the twosets of results are significantly correlated. The high cor-relation between the results from PLD and BLISS map-ping proves that the variation in eclipse depth is present in the data , not in the analysis approach, and it reflectsthe agreement found between the two methods in Ingallset al. (2016).Moreover, fitting our eclipse depths with a flat linegave a χ p-value of 8 . e −
11, indicating that a con-stant eclipse depth is an extremely poor fit to the results.We therefore conclude that the variability of 55 Cnc e ’seclipse depth as reported in D16 is genuine.D16 reported an increase from an average 47 ±
21 ppmdepth in 2012 to a 176 ±
28 ppm depth in 2013 at the4 σ level. Having agreed that the eclipse depth is in-deed variable, we turned our attention to modeling thedepth versus time, to see if our results support the sameaverage-increase model as proposed by D16. Heuristic Models
We generated five different models to explain thepattern 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.The first model is a flat eclipse depth, i.e. one that isnot varying in time. This is the model that we wouldanticipate for a “normal” planet, as we would expect tosee a mostly constant thermal flux from the planet versustime. There is only one free parameter for the flat model,that being the weighted average eclipse depth. We finda depth of 84 ± ppm best fits all the eclipses in theflat model.The second model is similar to the one proposed byD16, with year-to-year variability between the 2012 and2013 eclipses. We fit constants to our depths from 2012and 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 oftransition. For 2012, we find a depth of 24 ± ppm , anda depth of 138 ± ppm in 2013.The third model, which we refer to as the “spike”model, uses a flat depth everywhere except for the firsttwo observations in 2013, which we found were the twodeepest eclipses. The four free parameters of the spikemodel are the baseline eclipse depth, the depth of the“spike”, and two transition times. For this model, wefind a baseline depth of 58 ± ppm , and a depth of173 ± ppm for the first two observations in 2013. ariability in 55 Cnc e ’s Eclipse Depth Figure 10.
Eclipse depths versus time with our five different models overplotted. By calculating the Bayesian Information Criterion foreach model (shown in the left panel of each model), we are able to conclude that a constant eclipse depth model fits the results the leastwell out of all five. The number of free parameters, N p , is shown in the right panel of each model. The fourth model, which we call the “sloped” model, issimilar to the third, but with a linearly decreasing slopein 2013. It also has four free parameters: the average2012 eclipse depth, a time of transition, and a peak andslope value for the 2013 data. It uses the same 2012value as the year-to-year model, and perform a linearfit to the 2013 data of the form y = mx + b . We find m = − . ± . ppmday , and b = 164 . ± . ppm .The final model is a simple sine wave, given by Asin (2 πf t + δ ) + c , where A is the amplitude of thewave, f is the frequency, δ is a phase shift and c is aconstant offset from 0 (i.e. four free parameters). Todetermine this model, we fit a sine function to our best-fit eclipse depth results using Levenberg-Marquardt fit-ting, stepping over periods that ranged from half a dayto about three years. At each frequency step, we fit forsine wave amplitude, phase offset, and constant offset,and recorded the χ of each fit. We plot the results inFigure 11, with frequencies converted to periods in days.This model represents a periodic cycling of the planet’seclipse depth over time. Physically, a sine wave could bea useful approximation of any process (e.g. volcanism)that has a characteristic timescale, i.e. it represents abandwidth-filtered stochastic process. Figure 11 shows that a number of different sine waveperiods fit the data similarly well. This is because thedata are heavily aliased, as our sampling rate does notpermit us to distinguish between sine waves of many dif-ferent frequencies. In the bottom panel of the figure, weshow the sine model generated using the period circledin red in the top panel. Fluctuations of 300 ppm over a2 day period seem unphysical, yet fitting a sine wave toour aliased data finds that this is the “best” model.While we again acknowledge that the data are heav-ily aliased, we continue forward with a sine wave modelbecause future data could better constrain the periodic-ity (if periodic fluctuation is actually the physical realityof the system). We pursue a sine model with a less-extreme period, electing to use the best fit with period-icity greater than 20 days. This fit is circled in blue in 11.We found the model’s parameters to be A = 76 . ± . ppm , f = 0 . days − , δ = 1 . ± .
30, and c = 103 . ± . ppm . Model BICs
As in other works (e.g. Gibson 2010; Haynes et al.2015), we calculated the Bayesian Information Criterion(Schwarz 1978) to assess the relative quality of fit be-tween these models. The BIC is calculated through the2
Tamburo et al.
Figure 11.
Top : χ versus period for our sine wave model fitsto the eclipse depths. A number of different periods give modelswith similar χ , a result of the data being heavily aliased. Bottom :Visualization of the “best” sine wave fit (circled in red in the toppanel), which clearly represents an unphysical model and heavilyoverfits the data. For the purposes of examining a sine modelfurther, we use the best model with a period greater than 20 days.This fit is circled in blue in the top panel, and shown in the fifthpanel of Figure 10. formula
BIC = χ + k · ln ( n ), where k is the number ofmodel parameters (1 for the flat model, 4 for the spikemodel, 3 for the year-to-year model, 4 for the slopedmodel, and 4 for the sine model) and n is the numberof data points (8 in our case). A lower BIC is preferred,with an improvement of 6-10 points representing strongevidence for one model over another, and more than10 points representing very strong evidence (see Kass &Raftery 1995). We list the BICs for each model in Figure10.Based on our analysis, we found that the flat modelis clearly ruled out by its BIC score, despite using thefewest number of parameters. The remaining models allhave BICs that are within 6 of each other, so we are un-able to confidently distinguish between them. Becausethe year-to-year model has the fewest free parametersout of the variable models, we conclude that it is thebest model we currently have for interpreting the planet’seclipse depth variability. While the sine model has the best BIC score, we again emphasize that the fit was per-formed to heavily aliased data, and more observationswould be required to determine if a sine model is realis-tic.
Brightness Temperature Calculation
Following the calculations of Crossfield (2012), we con-verted the maxima and minima of our five models intoplanetary brightness temperatures ( T B ) using an ob-served spectrum of 55 Cnc. That paper demonstratedthat the observed spectrum gives a significantly lowerbrightness temperature for 55 Cnc e than would be cal-culated when treating the star as a blackbody. Usingtheir measured flux density of 3 . ± .
006 Jy, we nu-
Model Minimum T B (K) Maximum T B (K)Flat 1820 +142 − +142 − Spike 1550 +175 − +240 − Year-to-Year 1115 +266 − +183 − Sloped 1115 +266 − +242 − Sine 1163 +330 − +249 − Table 5
Maximum and minimum brightness temperatures for each of fivemodels. Temperatures were calculated using an observed stellarspectrum of 55 Cnc from Crossfield (2012). merically evaluated the following integral for values of T B : (cid:90) B λ ( T B ) S x ( λ ) dλ = δ occ R ∗ R p (cid:90) F ∗ λ S x ( λ ) dλ (2)...where B λ ( T B ) is the wavelength version of Planck’slaw as a function of the planet’s brightness temperature, S x ( λ ) is the relative system response in Spitzer IRACChannel 2 , F ∗ λ is the stellar flux density (converted to W m − sr − ), and δ occ is the planet’s measured eclipsedepth. The maximum and minimum values of T B foreach model can be found in Table 5.To examine the plausibility of our different models,we are motivated to explore what physical mechanismscould account for the brightness temperature variationssuggested by the year-to-year and sine models. Variability in Stellar Flux
The most likely mechanism forvariations in stellar flux on day-to-week timescales is theinfluence of star spots. Analyzing 11 years of photomet-ric observations of 55 Cnc, Fischer et al. (2008) foundmeasurements covering two consecutive rotational peri-ods ( ∼
80 days) that showed evidence of low-amplitude,short term-variability. This variability was attributed tostar spots covering less than 1% of the star’s surface, re-sulting in a brightness amplitude variation of 0.006 mag.While this variation could produce a signal of ∼
200 ppmin transit/eclipse observations (Demory et al. 2016b), thefact that Fischer et al. (2008) found 55 Cnc to be a quietstar makes this an unlikely scenario. Furthermore, thestar would have to be variable on the time scale of theeclipse itself in order to affect our results. On this ba-sis, we conclude that the star is not responsible for theobserved eclipse depth variability.
Asynchronous Rotation with a Hot Spot
Because the mea-sured variability is not well accounted for by changes inthe stellar flux, the variation must be due to changes inthe planet’s flux at 4.5 µ m. If the variability is truly peri-odic, one of the easiest physical explanations to visualizeis that 55 Cnc e is rotating asynchronously while main-taining a hotspot (e.g., from volcanism) on one side of theplanet. In this scenario, we would observe the occulta-tion of the planet’s hot, bright side in some observations,and its cooler, darker side in others. This would lead tomore significant drops in relative flux during some obser-vations than in others in a periodic fashion that matchesthe planet’s spin. While we cannot decisively constrain http://irsa.ipac.caltech.edu/data/SPITZER/docs/irac/ cali-brationfiles/spectralresponse/ ariability in 55 Cnc e ’s Eclipse Depth e , whichwould cause us to see the same face of the planet be-fore/after eclipse in every observation. We can estimatethis timescale from Eq. 9 given in Gladman et al. (1996): τ lock = ωa IQ Gm s k R p (3)where ω is the spin rate, a is the radial distance fromplanet to star, I is the planet’s moment of inertia aboutthe spin axis, Q is the planet’s specific dissipation func-tion, G is the gravitational constant, m S is the mass ofthe star, k is the tidal Love number of the planet, and R p is the mean radius of the planet. We use the estimateof Q = 100 from Gladman et al. (1996) for satellites inour solar system, and calculate k by k = 3 /
21 + µ ρgR p (4)(Burns et al. 1977) where µ is estimated as 3x10 Nm for rocky bodies (Gladman et al. 1996), ρ is the averageplanetary density (using M p from Demory et al. (2016b)and R p from this work), g is the surface gravity of theplanet, and R p is the planetary radius. Due to the smallsemi-major axis and high bulk rigidity of 55 Cnc e , de-spin times are on the order of hundreds of years, whichmakes unsynchronized rotation an unlikely explanationfor variability.Alternatively, it is possible that the planet could be-come captured in a higher-order spin-orbit resonance asit spun down. Mercury, with its 3:2 resonance (Pettengill& Dyce 1965), is an example of such a configuration inour own solar system. Further observations of 55 Cnc e will be required to resolve this possibility. Equilibrating Refractory Particulates
If we do not expectsignificant flux changes from the host star and we assumethat the planet is tidally locked, then thermal variabilitymust be occurring on the face of 55 Cnc e that we ob-serve before/after eclipse. Due to its close-in orbit andproximity to its host star, the planet most likely has amolten surface and may experience intense tidal forces,both of which could drive volatile loss through either sur-face evaporation and/or volcanic activity (Schaefer & Fe-gley 2009). Volcanic eruptions or a volatile atmospherecould result in particulates either lofted or condensedhigh above the surface (Mahapatra et al. 2017). D16 ex-amined this as a source of thermal variability, with vol-canic plumes raising regions of the planet’s photosphereto higher, cooler regions in 2012, which lead to smallereclipse 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 whichthis doesn’t apply because the depth is constant), we as-sume that our measured value of T b,max represents thetemperature of the radiating planetary surface, and wecheck to see if the planet’s equilibrium temperature (dueto stellar insolation alone) is consistent with this mea-surement. If the planet has a molten surface or is vol-canically active, it may outgas material that cools at a height above the planet, blocking out the hot surface andleading to lower temperatures. For each model, we checkto see if T b,min can be achieved by this obscuring volcanicejecta.To start, we determined the equilibrium temperatureof 55 Cnc e due to stellar insolation alone over a rangeof planetary albedos. We did this both for an efficientcooling case, in which the planet is able to radiate heataway over its whole surface, and an inefficient case, inwhich it can only radiate over half the surface.Next, for each model, we calculated the equilibriumtemperature of refractory particulates at a height of 100km above the planet’s surface due to both solar andplanetary radiation over a range of grain albedos. Ineach case, we assumed that the planet is radiating atits T b,max temperature. Our assumption of the heightof the grain layer above the surface is about twice themaximum theoretical height of volcanic plumes on earth(Wilson et al. 1978), which we feel is not unrealistic for aplanet of ∼ ⊕ that may experience intense tidal heat-ing. In practice, however, we found that our results forthis calculation are not strongly dependent on the loca-tion of the layer of grains. The results of these calcu-lations are shown in Figure 12. The top panel showsthe equilibrium temperature for efficient and inefficientcooling models compared to T b,max of the model, whilethe bottom shows the temperature of grains comparedto T b,min .While results vary for the different models, they all re-quire a dark planetary surface (albedo < > σ ranges of T b,max and T b,min basedon planetary equilibrium temperature alone. They alsorequire the inefficient cooling case to agree with T b,max ,which may be indicative of a lack of an atmosphere forefficient heat redistribution.However, it is possible that the planet’s thermal ra-diation budget is not limited to incoming flux from itsstar. If the four other planets in the system (Fischer etal. 2008) help maintain even a slightly eccentric orbit,55 Cnc e may dissipate a significant amount of radiationin the form of tidal heating. Demory et al. (2016b) ex-amined tidal dissipation as a possible explanation for a3100 K region in their longitudinal temperature map of55 Cnc e , using the Mercury-T code of Bolmont et al.(2015). They found that tidal heat fluxes ranging from10 − to 10 W m − are achievable for the planet, cor-responding to temperature increases of a few Kelvin to ∼ T b,max tobe consistent with the efficient cooling model.Budaj et al. (2015) tabulated the opacities and equilib-rium temperatures of several species of dust grains thatmay be present in exoplanet atmospheres over a wave-length range from 0.2 to 500 µ m and grain radii from0.01 to 100 µ m. Extracting their opacity results at 4.5 µ m from their tables , we calculated the single scattering ∼ budaj/dust/ Tamburo et al.
Figure 12.
Refractory particulate temperature calculations for our four variability models. The top panel in each shows planetaryequilibrium temperatures for efficient (blue) and inefficient (red) cooling cases over a range of planetary albedos. The 2 σ range of T b,max for each model is shown in orange. The bottom panel shows the equilibrium temperature of refractory grains (green) at a height of 100km above the planet’s surface, due to radiation from the planet (at T b,max ) and the star. The 2 σ range of T b,min is shown in blue. For allfour models, we find that T b,max can be achieved by a dark planet with inefficient heat redistribution, and T b,min could be achieved withreflective grains obscuring the hot planetary surface. Values for T b,max and T b,min are given in Table 5. albedo (Eq. 23 from their paper) for all grain sizes fromall species they examined. We used their temperaturetables for a star of T eff = 5000 K, and solid angle of0.01 and 0.03 steradians (55 Cnc subtends a solid angleof ∼ e ’s sky).Cross-referencing these results, we identified the com-position and size of grains that have the requisite albedosand equilibrium temperatures to match our brightnesstemperatures. These include alumina, olivine, and py-roxene, with grain radii ranging from around 0.01 to 0.6 µ m. If the planet regularly condenses or outgasses thesematerials in an opaque cloud layer, its composition couldachieve the appropriate albedo and equilibrium temper-ature to match T b,min from our models.Mahapatra et al. (2017) applied a kinetic, non-equilibrium cloud formation model to the atmosphere of55 Cnc e using several different abundance cases. UsingD16’s derived T-P profile, they found that 55 Cnc e cansupport a highly opaque cloud layer consisting of mainlyMg-silicates and Si-oxides. The local gas temperature ofthis cloud layer (see panel 1 of Fig. 10 in their paper) ranges from about 850-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 por-tions of the planetary surface to be blocked by a layerof cooler refractory grains, depending on the local struc-ture of planetary surface temperature and the tempera-ture of occulting material. Determining whether or notthe maintenance of such a large layer of condensed orejected material is feasible is beyond the scope of thissimple calculation. Future Work
As we’ve discussed, more closely-spaced observationsof 55 Cnc e ’s eclipse will be needed to constrain thetimescale of thermal emission variability on the planet.The sine model created in Section 4.2.1 represents theshort end of this timescale, with a period of ∼
35 days.This model can be ruled in or out with Nyquist samplingof the planet’s eclipse depth at 4.5 µ m, by taking obser-vations spaced at most 17.5 days apart over multiple pe- ariability in 55 Cnc e ’s Eclipse Depth e . As Kaltenegger & Sasselov (2010) discuss, it may bepossible to detect features of geochemical cycles indica-tive of volcanic activity on super-Earths through spec-troscopy with the James Webb Space Telescope (Gardneret al. 2006).Unfortunately, 55 Cnc’s brightness (mag. 4.0 in K s band, Skrutskie (2006)) will saturate many observingmodes of JWST. Recent performance tests by Greeneet al. (2016), however, indicate that spectroscopic ob-servations of the system will be achievable using NIR-Cam’s F444W filter. Tests performed with PandExo(Batalha et al. 2017), a tool that simulates observationswith JWST , confirm that observations of 55 Cnc are pos-sible with NIRCam without saturating. The F444W fil-ter, which offers wavelength coverage from ∼ µ m,will not likely constrain the presence of silicate materials,which have the strongest spectral features between 8 and12 µ m (Hu et al. 2012). Nevertheless, observations withNIRCam will be useful for further probing the variabilityof 55 Cnc e and for searching for features indicative ofthe presence of an atmosphere.Additionally, a mode proposed by Schlawin (2017)using NIRCAM’s Dispersed Hartmann Sensor (DHS)may be able to take observations of comparably brightobjects, as it spreads light over 10 separate spectra.Whether this proposed mode will be approved, andwhether or not it will be able to image 55 Cnc withoutsaturating, remain to be seen.If selected to go forward, the 55 Cnc system may alsobe observable with the proposed ARIEL mission (seeTinetti et al. 2016), which will be designed to monitorthe atmospheres of exoplanets, some of which will be inthe super-Earth size range. Such a mission would alsoallow us to better constrain the presence or lack of anatmosphere on 55 Cnc e . CONCLUSIONS