Variability in the Atmosphere of the Hot Jupiter Kepler-76b

TL;DR

This study analyzes over 1,000 days of Kepler data to reveal atmospheric variability in the hot Jupiter Kepler-76b, including changes in emission, reflectivity, and phase curve offset, linked to atmospheric dynamics.

Contribution

It provides the first detailed evidence of atmospheric variability in Kepler-76b using long-term Kepler observations, highlighting dynamic thermal and cloud processes.

Findings

Detected variability in phase curve amplitude between 35 and 70 ppm

Observed a consistent eastward phase offset of 9 degrees

Identified potential thermal and cloud structure changes over tens of days

Abstract

Phase curves and secondary eclipses of gaseous exoplanets are diagnostic of atmospheric composition and meteorology, and the long observational baseline and high photometric precision from the Kepler Mission make its dataset well-suited for exploring phase curve variability, which provides additional insights into atmospheric dynamics. Observations of the hot Jupiter Kepler-76b span more than 1,000 days, providing an ideal dataset to search for atmospheric variability. In this study, we find that Kepler-76b's secondary eclipse, with a depth of $87 \pm 6$ parts-per-million (ppm), corresponds to an effective temperature of 2,830$^{+50}_{-30}$ K. Our results also show clear indications of variability in Kepler-76b's atmospheric emission and reflectivity, with the phase curve amplitude typically $50.5 \pm 1.3$ ppm but varying between 35 and 70 ppm over tens of days. As is common for hot Jupiters, Kepler-76b's phase curve shows a discernible offset of $\left( 9 \pm 1.3 \right)^\circ$ eastward of the sub-stellar point and varying in concert with the amplitude. These variations may arise from the advance and retreat of thermal structures and cloud formations in Kepler-76b's atmosphere; the resulting thermal perturbations may couple with the super-rotation expected to transport aerosols, giving rise to a feedback loop. Looking forward, the TESS Mission can provide new insight into planetary atmospheres, with good prospects to observe both secondary eclipses and phase curves among targets from the mission. TESS's increased sensitivity in red wavelengths as compared to Kepler means that it will probably probe different aspects of planetary atmospheres.

Tables

Table 1. Table 1: Transit, Eclipse, and Other Photometric Variation Model Parameters

Transit
$T_0$ (BJD)	$2454965.00396\pm 5\times {10}^{-5}$	Time of Primary Transit
Period (days)	$1.54492871\pm 9\times {10}^{-8}$	Orbital Period
$R_p/R\star$	$0.085\pm 0.001$	Planetary Radius/Stellar Radius
$a/R_{\star }$	${5.103}_{-0.060}^{+0.058}$	Orbital Semi-Major Axis/Stellar Radius
$b$	$0.908\pm 0.003$	Impact Parameter
$i$ (degrees)	$79.7\pm 0.1$	Orbital Inclination
$({\gamma }_1,{\gamma }_2)$	(0.313, 0.304)	Limb-Darkening Coefficients (Faigler et al., 2013)
Eclipse/Photometric Variations
$D$ (ppm)	$87\pm 6$	Secondary Eclipse Depth
$A_{\mathrm{planet}}$ (ppm)	$50.5\pm 1.3$	Phase Curve Amplitude
$\delta$ (degrees)	$-9.0\pm 1.3$	Phase Curve Offset
$A_{\mathrm{ellip}}$ (ppm)	$13\pm 1$	Ellipsoidal Variation (EV) Amplitude
$A_{\mathrm{beam}}$ (ppm)	$3.8\pm 0.3$	Beaming Effect Amplitude
Other Parameters
$M_{\star }\left(M_{\odot }\right)$	$1.2\pm 0.2$	Stellar Mass (Faigler et al., 2013)
$g$	0.07	Gravity-Darkening Parameter (Claret & Bloemen, 2011)
${\alpha }_{\mathrm{beam}}$	0.92	Beaming Coefficient (Faigler et al., 2013)
${\alpha }_{\mathrm{ellip}}$	1.02	EV Coefficient
$K_z$ (km/s)	$0.306\pm 0.020$	Radial Velocity (Faigler et al., 2013)
$q_{\mathrm{ellip}}$	$\left(1.9\pm 0.2\right)\times {10}^{-3}$	Planet-Star EV Mass Ratio
$M_{p,\mathrm{ellip}}\left(M_{\mathrm{Jup}}\right)$	$2.4\pm 0.3$	Planetary EV Mass
$q_{\mathrm{beam}}$	$\left(1.7\pm 0.1\right)\times {10}^{-3}$	Planet-Star Beaming Mass Ratio
$M_{p,\mathrm{beam}}\left(M_{\mathrm{Jup}}\right)$	$2.1\pm 0.2$	Planetary Beaming Mass

Equations

$$\Delta F=F_0-A_{\mathrm{ellip}}\cos \left(2\times 2\pi \varphi \right)+A_{\sin }\sin \left(2\pi \varphi \right)-A_{\cos }\cos \left(2\pi \varphi \right),$$

$$\ln L=-\frac{1}{2}\sum_{i}\bigg{(}\frac{d(t_{\rm i})-m(t_{\rm i})}{\sigma}\bigg{)}^{2}+\ln\left(2\pi\sigma^{2}\right),$$

Full text

planets and satellites: atmospheres – planets and satellites: gaseous planets – planets and satellites: individual (Kepler-76b)

††journal: ApJ, 2019 Jan 31††software: numpy (van der Walt et al., 2011), astropy (Astropy Collaboration et al., 2013, 2018), statsmodels (Seabold & Perktold, 2010), scipy (Jones et al., 2001–), lightkurve, PyAstronomy, emcee (Foreman-Mackey et al., 2013), get_lds (Espinoza & Jordán, 2015)

1 Introduction

The age of exoplanet meteorology has arrived – from the first detection of atmospheric thermal emission (Deming et al., 2005; Charbonneau et al., 2005) to the recent detection of the secondary eclipses of an Earth-sized planet (Demory et al., 2016), observation, characterization, and modeling of exoplanet atmospheres has blossomed into a rich field approaching that of solar system meteorology in scope and complexity. Indeed, the much wider range of dynamical, compositional, and radiative environments across which exoplanets have been discovered means the study of these systems will undoubtedly provide insight into the planetary atmospheres in our own solar system.

Among the first key results in the study of exoplanet meteorology was the mapping of atmospheric emission from the transiting hot Jupiter HD189733b (Knutson et al., 2007). In that study, observations in the Spitzer Space Telescope’s 8 $\mu$m band revealed a deep secondary eclipse (0.3381% in differential flux) corresponding to a dayside brightness temperature of more than 1,200 K. Perhaps of more interest to meteorological studies, the data also showed signs of a planetary emission during the entire orbit, peaking as the tidally locked planet’s dayside rotated into view. However, instead of being centered around the eclipse, the planet’s phase curve peak occurred about 2 hours before, indicating the brightest and hottest region in the atmosphere was shifted about 16∘ eastward from the sub-stellar point. Knutson et al. (2007) attributed this offset to advective transport within the atmosphere, a result anticipated by earlier studies of hot Jupiter atmospheres (e.g., Showman & Guillot, 2002).

The mountain of data made available by the Kepler Mission has allowed analysis of phase curves for numerous planets (e.g., Esteves et al., 2015). However, since the Kepler data involve visible wavelengths (with a bandpass spanning 400 to 900 nm), these phase curves likely involve a combination of thermal emission and reflected light from high-altitude aerosols (Sudarsky et al., 2000). The aerosols obscure the lower atmosphere, complicating the determination of atmospheric composition (Sing et al., 2016).

Although the composition, size distribution, and vertical mixing of aerosols are difficult to determine a priori, their influence on the planetary phase curve, in particular the resulting phase curve offset, may be diagnostic. Parmentier et al. (2016) summarize the distribution of observed phase curve offsets and highlight a clear correlation with planetary effective temperature. Cooler planets show offsets toward the west, consistent with the expectation that clouds dominate these planets’ phase curves. By contrast, some hotter planets, such as HD189733 b, show eastward offsets, symptomatic of either homogeneous or little dayside cloud cover and downwind transport of thermally emitting regions, although there are important exceptions, such as CoRoT-2b (Dang et al., 2018). Parmentier et al. (2016) also suggest that the same super-rotating winds that drive hot regions eastward can draw aerosols that form on the nightside across a hot planet’s western (night-to-day) terminator, forming a dayside cloud front whose eastern-most boundary depends on the aerosol composition – the more refractory the aerosols, the farther east the boundary. These clouds can contribute to the planetary phase curve and, since they may be confined to the western hemisphere, they may compete with eastern-hemisphere emission and drive the observed phase curve offset toward the west.

All these effects take place in a highly dynamic atmosphere, in which winds approach sonic speeds (Koll & Komacek, 2018), and the interplay between transport of aerosols, which act to cool the atmosphere, and potent stellar insolation, which drives the winds, produce variable atmospheric structures. This variability likely shows up in the reflected and emitted light from the planets’ atmospheres. Studying aerosol transport in the atmosphere of HD209458b, Parmentier et al. (2013) predicted variations in the secondary eclipse depths as large as 50% in IR wavelengths over tens of days. Analyzing hundreds of days of Kepler observations of the hot Jupiter HAT-P-7 b, Armstrong et al. (2016), indeeed, found significant variations in the phase curve offset, which wandered back and forth across the sub-stellar point over tens of days, likely driven, at least in part, by variations in aerosol transport. Although Armstrong et al. (2016) were not able to robustly detect them, variations in the phase curve amplitude probably accompanied the offset fluctuations since variations in aerosol transport would also be expected to modify the planet’s total geometric albedo. Hindle et al. (2019) explored the effects of magnetic fields on hot Jupiter phase curves and found that planetary-scale equatorial magneto-Kelvin waves could drive westward tilting eddies and produce the observed periodic westward offset in HAT-P-7b’s phase curve reported in Armstrong et al. (2016).

In principle, it should be possible to connect variations in a planet’s atmospheric reflection and emission to phase curve variability, but this requires a large enough population of planets that spans a range of conditions to piece together a cogent story. Following this thread, we analyzed Kepler observations of the hot Jupiter Kepler-76b, a two Jupiter-mass (${\rm M_{Jup}}$) planet with an orbital period of 1.5 days around a hot (6,300 K) star (Faigler et al., 2013). In addition to having a clear atmospheric signal, Kepler-76b also induces observable ellipsoidal variations on its host star (the photometric signature of tidal bulges on the star rotating in and out of view – Welsh et al., 2010) and a beaming signal due to the radial velocity of the star (Loeb & Gaudi, 2003). With nearly 1,400 days of observation, we find significant variability in the planet’s phase curve amplitude, $A_{\rm planet}$, and offset, $\delta$, qualitatively consistent with the variations seen for HAT-P-7 b. We have also updated the transit parameters for Kepler-76b.

This article is organized as follows: in Section 2, we describe our process for conditioning the Kepler dataset and then our methodology for modeling Kepler-76b’s transit, eclipse, and other photometric signals. In Section 3, we discuss our results, and in Section 4, we discuss their implications and future prospects for using TESS data to conduct similar analyses.

2 Data Analysis

2.1 Data Conditioning

We applied the following steps to condition all fourteen available quarters (Quarters 1-5, 7-9, 11-13, 15-17) of Kepler long-cadence (i.e., 30-min integration times) data:

To reduce variations in each quarter’s data to the same scale, we subtracted the median value from each quarter’s PDCSAP_FLUX time-series before dividing through by that same median value. 2. 2.

We next applied a median boxcar filter with a window size equal to four orbital periods. We experimented with window sizes from 1-15 orbital periods (holding fixed the transit parameters to those reported in Faigler et al. 2013), then fitting the out-of-transit photometric signals described below using the Levenberg-Marquardt algorithm (LM) (Newville et al., 2014), including the planet’s eclipse and phase curve; four orbital periods maximized the eclipse depth and minimized the distortion to the phase curve using all available data. To mitigate edge-effect distortions from our boxcar filter, we extended the time-series out a full window length beyond both ends by reflecting the original time-series across its boundary. 3. 3.

Finally, we stitched each quarter’s conditioned data into one long time-series.

Figure 1 illustrates a portion of the raw and conditioned time-series for Kepler-76b.

2.2 Analyzing All Data Folded Together

After applying the above conditioning procedure, we fit a series of photometric models to all the available data. We first analyzed the out-of-transit portion of Kepler-76b’s light curve since those signals (the ellipsoidal variations, the planet’s phase curve, and the Doppler beaming signal) are independent of the transit and can, in principle, contribute to the transit portion (although we found their effects are negligible).

We folded all data on the best-fit period reported by Faigler et al. (2013) and masked out the transit and eclipse portions of the light curve and fit the following model:

[TABLE]

where $F_{0}$ represents a constant baseline and $\phi$ the orbital phase ($\phi=0$ at mid-transit). The second term represents the ellipsoidal variations induced by the planet’s tidal gravity (Welsh et al., 2010), with $A_{\rm ellip}$ its amplitude. The third and fourth terms are a combination of the Doppler beaming signal (Loeb & Gaudi, 2003) and the planet’s reflected/thermally emitted phase curve, with allowance for an offset in the curve’s maximum from superior conjunction at $\phi=0.5$ (Faigler et al., 2013). We also allowed the per-point uncertainty (noise) to be a free parameter and estimated the model likelihood $L$ as

[TABLE]

where $d(t_{\rm i})$ is the datum at time $t_{\rm i}$, $m(t_{\rm i})$ the model for that same time, and $\sigma$ the per-point uncertainty. With this likelihood, we used a Markov Chain Monte Carlo (MCMC) analysis (Foreman-Mackey et al., 2013) with 100 walker chains, each with 5,000 links and a burn-in phase of 2,500 links. We aimed for convergence by requiring small auto-correlation times (e.g., Geyer, 1992) for the mean chain of each fit parameter (none exceeded 80 links) and all Geweke Z-scores (Geweke, 1992) of about 3 or less.

Traditionally, the beaming and planetary phase curve signals are treated separately through the combination A_{\rm beam}\sin\left(2\pi\phi\right)-A_{\rm planet}\cos\big{(}2\pi\left(\phi-\delta\right)\big{)}. This expression can be re-cast using $A_{\rm sin}=A_{\rm beam}-A_{\rm planet}\sin\left(2\pi\delta\right)$ and $A_{\rm cos}=A_{\rm planet}\cos\left(2\pi\delta\right)$. Since the beaming and planetary phase curve signals have the same frequency (once per orbit) and the phase offset $\delta$ allows them to phase up (at least in principle), they are degenerate; completely independent fits for $A_{\rm beam}$, $A_{\rm planet}$, and $\delta$ are not possible.

Figure 2 shows the fit and residuals for Equation 1 to the dataset. Figure 3 shows the resulting best-fit parameters, while Figure 4 shows the corresponding range of values for $A_{\rm planet}$ and $\delta$. To generate the latter figure, we calculated $A_{\rm beam}$ from the radial-velocity semi-amplitude $K_{\rm z}=0.306\pm 0.020$ km/s reported in Faigler et al. (2013) using $A_{\rm beam}=4\ \alpha_{\rm beam}$ and $K_{\rm z}/c=3.8\pm 0.3$ ppm, where $c$ is the speed of light (Loeb & Gaudi, 2003) – we did not use the $A_{\rm beam}$-value reported in Faigler et al. (2013). Then, we Monte-Carlo-sampled 250,000 times from a Gaussian distribution for $A_{\rm beam}$ (with the mean and width given above), as well as the MCMC chains for $A_{\rm sin}$ and $A_{\rm cos}$. Figure 4 shows the clear correlation between $A_{\rm beam}$ and $\delta$. As a final check, we fit the out-of-transit/out-of-eclipse data with an LM approach, holding $A_{\rm beam}$ fixed at 13.5 ppm, as reported in Faigler et al. (2013), and were able to recover the reported $A_{\rm planet}$ but not the $\delta$-value.

Next, we analyzed the in-transit portion of the light curve. Importantly, we accounted for Kepler’s finite integration time of 30-min-per-point by super-sampling each point, i.e. we modeled the light curve at a cadence of 3-min-per-point and then downsampled to 30-min-per-point (cf. Kipping, 2010) (the finite integration had no significant effect on the Equation 1 model). By first subtracting the Equation 1 model from the in-transit portion of the light curve, we removed its (very small) influence on the transit. Next, we used LM to fit a standard transit curve with quadratic limb-darkening (Mandel & Agol, 2002)111As implemented in PyAstronomy – https://github.com/sczesla/PyAstronomy. to that portion of the time-series within two transit durations of the reported mid-transit time, allowing the out-of-transit baseline to vary. This analysis netted initial best-fit transit parameters: $a/R_{\star}$, the ratio of the planet’s semi-major axis to the stellar radius; $b$, the impact parameter; and $R_{\rm p}/R_{\star}$, the planet-to-star radius ratio. We assumed the orbital eccentricity is equal to zero based on the photometric and radial velocity analyses in Faigler et al. (2013) and held fixed the limb-darkening coefficients $\gamma_{1}=0.313$ and $\gamma_{2}=0.304$. We explored two other sets of quadratic limb-darkening coefficients, estimated assuming the stellar parameters from Faigler et al. (2013) and using the model from Espinoza & Jordán (2015) but found different values had negligible impact on our transit analysis.

To check for transit-timing variations, we fit each individual transit (with the full orbital light curve this time and not just the in-transit portion) using LM again and held fixed all transit parameters except the mid-transit times and out-of-transit baseline. For uncertainties, we took the square root of the diagonal elements of the resulting covariance matrix (Press et al., 2007, p. 790) and confirmed the accuracy of the mid-transit times and uncertainties by fitting transit curves with an MCMC analysis (Foreman-Mackey et al., 2013) to several individual transits – we achieved convergence with 100 walkers, each a chain of 1,000 links and a burn-in phase of 300 links. Although most single transits are detectable, they usually comprise one or two data points and so provide no useful constraints on parameters other than the mid-transit times. Finally, we fit the collection of mid-transit times with both linear and quadratic ephmerides using an MCMC analysis (100 walker chains, 5,000 links in each with a burn-in of 1,500 links) and found no significant departure from the latter. We also found no significant periodicities using a Lomb-Scargle periodogram (Lomb, 1976; Scargle, 1982). Figure 5 shows the resulting transit ephemeris, which is consistent to within uncertainties with that reported in Faigler et al. (2013).

With this updated ephemeris, we folded together all available data and fit the resulting transit portion of the light curve, with free parameters $a/R_{\star}$, $b$, $R_{\rm p}/R_{\star}$, and an estimate of the per-point noise (again using Equation 2 for the likelihood). We used the same MCMC analysis as above, with 100 walker chains, each with 32,000 links and a burn-in phase of 16,000 links. We aimed for convergence by requiring small auto-correlation times (e.g., Geyer, 1992) for the mean chain of each fit parameter (none exceeded 250 links) and all Geweke Z-scores (Geweke, 1992) of about 3 or less. Figure 6 shows the resulting parameter distributions, with the best-fit values and uncertainties given in Table 1, and Figure 7 shows the resulting transit curve and residuals.

2.3 Searching for Eclipse and Planetary Phase Curve Variability

We turned to a search for orbit-to-orbit variability in the planet’s phase curve and secondary eclipse. First, we folded together all 944 orbits-worth of data on our best-fit ephemeris and used the LM algorithm to fit the out-of-transit and eclipse portions (accounting for Kepler’s finite integration time) to establish the average fit-values – the resulting eclipse depth $D=87\pm 6$ ppm, a 14.5-$\sigma$ detection. In our search for variability, ideally, we’d analyze the data for each orbit individually to maximize our time resolution. However, with a 14.5-$\sigma$ detection over 944 orbits, we expect that a 3-$\sigma$ detection of the eclipse in the presence of Gaussian noise requires folding together data from about 40 ($=\left(\frac{3}{14.5}\right)^{2}\times 944$) orbits, giving a time resolution of about 60 days. We expect estimates of the parameters $A_{\rm sin}$ and $A_{\rm cos}$ to be more robust, however, since those signals are present throughout most of the orbit (the eclipse only occupies one or two points each orbit).

For this analysis, we marched a window spanning 10, 20, or 40 orbits from the beginning to the end of the dataset, moving one orbital period at a time. We folded all the out-of-transit points in the window together and fit both Equation 1 and an eclipse, which is represented by the same transit light curve model as above but with a variable depth $D$ and limb-darkening parameters set to zero (i.e., a uniform disk). We used the LM algorithm and took the square root of the diagonal of the covariance matrix scaled by the square root of the resulting reduced $\chi^{2}$ as the parameter uncertainties – we compared these uncertainties to those derived from an MCMC analysis for several stacked orbits and found good agreement (to within 1 ppm).

To avoid having our window span large gaps in the data, we required a window to include a large enough number of data points that at least half the desired number of orbits were represented. For example, for the window spanning 10 orbits, we required at least 5 orbital periods worth of data (i.e., about $360=5\times 1.5\ {\rm days}/30\ {\rm min}$ points). We took the mid-eclipse time to be the mid-transit time plus half an orbital period, again assuming zero orbital eccentricity. Figure 8 shows the resulting collection of best-fit parameters for the different windows, where the x-coordinate of each point represents the median observational time for all the points in the window. Since the windows span several orbits and the points are only spaced by one orbital period, the parameter fits are not all statistically independent.

The plots seem to show considerable variability, but comparing the variations to the best-fit average values (top right panel of Figure 8) shows that the eclipse depth departs from its average value, 87 ppm, by less than 3-$\sigma$, where $\sigma$ is calculated by adding uncertainties in quadrature (Taylor, 1997, p. 58). By contrast, $A_{\rm sin}$ and $A_{\rm cos}$ frequently depart by more than 3- and sometimes 5-$\sigma$ from their average values. As expected, though, the more orbits we stacked together, the more the variations are averaged out, but even the points with 40 orbits stacked together (in green) show statistically significant variation. These results, of course, rely on the accuracy of our uncertainty estimates, so we conducted several numerical tests to check them.

First, we generated and fit several synthetic datasets incorporating Gaussian noise and with the same observational times as the real dataset but holding the phase curve parameters and eclipse depth constant. These fits showed variations visually similar to those in Figure 8, but always the variations in fit parameters were within 3-$\sigma$ of the assumed constant values.

We next gathered and phased up about 7,000 datasets involving 10, 20, or 40 orbits worth of out-of-transit data (10 times the number of points in the left panels in Figure 8) but randomly scattered throughout the whole Kepler time-series, i.e., orbits that were not necessarily adjacent in time. We expect that parameter fits to these scattered datasets will show only random scatter, with orbit-to-orbit coherent variations averaged out. We used the Kolmogorov-Smirnov (KS) test (Press et al., 2007, p. 730) to compare the distributions of best-fit parameters (scaled by their respective uncertainties) for these scattered data to the same scaled distributions for the unscattered data (right panels in Figure 8). To judge the range of KS scores we should expect, we randomly sub-sampled the distribution of best-fit parameters for the scattered data 10,000 times and compared these sub-samples to one another. The KS test comparisons indicated that the distribution of eclipse depths for the unscattered data shown in Figure 8 resembled closely the distribution for the scattered data, indicating that the variations in the top left panel of Figure 8 are consistent with statistical noise. The KS test comparisons for $A_{\rm sin}$ suggests the distribution for the unscattered data is probably inconsistent with random scatter (with a KS probability of about 0.1%), and the KS test comparisons for $A_{\rm cos}$ are totally inconsistent with random scatter (KS probabilities $<10^{-11}$ when 40 orbits are stacked together). This trend is not surprising. We typically detected the $A_{\rm cos}$ signal at somewhat greater significance than the $A_{\rm sin}$, and both were detected at much greater significance than the eclipse.

We also employed the same orbit stacking and phasing and fit transit radii (holding other transit parameters fixed) to look for correlations with the phase curve parameters. If the variations in the phase curve arise from some stellar or instrumental effect, we might expect the transit depth to vary in concert; however we found no robust linear correlation.

As a final check, we performed an injection-recovery test for which we generated a synthetic dataset designed to closely mimic the original dataset and included the variations in phase curve parameters shown in Figure 8. For this test, we divided up the full time-series into individual orbits, calculating the mid-transit time for each orbit. We linearly interpolated to that mid-transit time from among the phase curve parameters for 10 orbits stacked together (blue points in Figure 8) to determine the phase curve parameters for the orbit. We generated a phase curve model for that orbit and subtracted it from the original, unconditioned time-series (i.e., from the data shown in the top panel of Figure 1). Then, we shifted all the mid-transit times forward in time by about 800 days, i.e., half the maximum observational time in the Kepler dataset – mid-transit times that moved out beyond this maximum time, we wrapped back around to the beginning of the dataset. We next injected synthetic phase curve signals back into these shifted data but this time using the phase curve parameters from the original mid-transit times. This approach has the effect of retaining the time-ordered structure of the phase curve variations while planting the phase curves in a different noise environment. Presumably, if we can recover the phase curve parameters, they are robust against noise. Indeed, after fully re-conditioning these shifted, synthetic data, we were able to recover the vast majority of parameters to within uncertainties. Not surprisingly, for the points not recovered successfully, most were associated with the eclipse depth, which was detected with the least signal-to-noise originally.

3 Results

Our analysis provides an updated set of transit parameters for the Kepler-76 planetary system that was originally announced in Faigler et al. (2013); our ephemeris agrees well with that previous study, while our transit parameters are different. For instance, Faigler et al. (2013) estimated $R_{\rm p}/R_{\star}=0.0968\pm 0.0003$, while we estimate $R_{\rm p}/R_{\star}=0.085\pm 0.001$, an 11-$\sigma$ discrepancy. We conducted a variety of tests (ignoring the finite-time Kepler integration, applying an alternative data conditioning approach, etc.) and were unable to reproduce the previously reported value. Heller (2019) provided equations for calculating the transit depth for a limb-darkened star. Using these equations and the transit parameters given in Faigler et al. (2013), we find that the depth expected (0.66%) does not seem to agree with the depth of the transit light curve from our study or that shown in Figure 6 of Faigler et al. (2013) (about 0.56%).

In any case, we corroborate detection of planet-induced ellipsoidal variations and the beaming signal, albeit with a different approach and results. Our amplitude for the ellipsoidal variations $A_{\rm ellip}$ is 4-$\sigma$ smaller than that of Faigler et al. (2013). For the beaming signal, we retrieve an amplitude $A_{\rm beam}$ almost 5-$\sigma$ discrepant, although this disagreement may simply arise from our different approach (fitting the planetary phase curve and beaming signals in a non-degenerate way). Using our new transit parameters, we can estimate a new planet-star mass ratio $q$ through the relation $A_{\rm ellip}=\alpha_{\rm ellip}\left(q\sin^{2}i\right)\left(R_{\star}/a\right)^{3}$. In this equation, $\alpha_{\rm ellip}$ depends on the gravity-darkening parameter (0.07 – Claret & Bloemen, 2011) and linear limb-darkening parameter (0.552 – Espinoza & Jordán, 2015) as described in Morris (1985). With these values, $\alpha_{\rm ellip}=1.02$. By re-casting the $A_{\rm ellip}$ equation above in terms of mass ratio $q$ and sampling the distributions of fit parameters from our analysis, we can convert the ellipsoidal signal into $q_{\rm ellip}=\left(1.9\pm 0.2\right)\times 10^{-3}$. In a similar way, we can also estimate a mass ratio for the beaming signal, whose amplitude $A_{\rm beam}=\alpha_{\rm beam}4\left(K_{\rm Z}/c\right)$. Taking $\alpha_{\rm beam}=0.92$, as in Faigler et al. (2013), we find $q_{\rm beam}=\left(1.7\pm 0.1\right)\times 10^{-3}$, which is in good agreement with our estimate for $q_{\rm ellip}$. With the stellar mass $M_{\star}=1.2\pm 0.2$ solar masses $M_{\odot}$ from Faigler et al. (2013), these ratios give planetary masses $M_{\rm p,ellip}=2.4\pm 0.3$ Jupiter masses ${\rm M_{Jup}}$ and $M_{\rm p,beam}=2.1\pm 0.2\ {\rm M_{Jup}}$, respectively. These values agree with the mass estimates in Faigler et al. (2013).

Turning to the results from our phase curve analysis, a priori, we expect that the dayside temperature of Kepler-76b, in radiative equilibrium, lies between 2,000 and 2,300 K (between the extremes of complete atmospheric redistribution and a uniform temperature dayside with no redistribution to the nightside), either of which puts a significant fraction of its blackbody curve within the Kepler bandpass. If the planet’s secondary eclipse represented only thermally emitted light, these temperatures correspond to a depth between 4 and 20 ppm. However, Kepler-76b’s average eclipse depth $D=87\pm 6$ ppm (within 1.3-$\sigma$ of the result in Faigler et al. 2013). With a nightside flux equal to zero, $D$ should be equal to $2\times A_{\rm planet}\cos\left(\phi-\delta\right)$, the planet’s phase curve at mid-eclipse. Converting our average $A_{\rm cos}$ and $A_{\rm sin}$-values, we find that these values agree to within 2.3-$\sigma$, consistent with no nightside flux. This eclipse depth corresponds to a brightness temperature of 2,830${}^{+50}_{-30}$ K and suggests there is a significant reflected component, $\geq 67$ ppm. Following Faigler & Mazeh (2011), we estimate a Kepler bandpass-integrated geometric albedo $p_{\rm geo}\geq 0.25$, also consistent with Faigler et al. (2013).

Finally, we turn to the apparent variations in the planet’s phase curve. From Figure 8, variations in $A_{\rm cos}$ and $A_{\rm sin}$ amount to variations of about 40% in $A_{\rm planet}$ around a median value of $50.5\pm 1.3$ ppm. Meanwhile, $\delta$ has a median value of $-9^{\circ}\pm 1.3^{\circ}$ (an offset east of the sub-stellar point) and swings between about $-40^{\circ}$ (an eastward offset) and $+20^{\circ}$ (a westward offset). Figure 9 shows the range of $A_{\rm planet}$- and $\delta$-values, along with a linear fit which incorporated the per-point uncertainties in $A_{\rm planet}$ and $\delta$ (Boggs & Rogers, 1990).

The Kepler-76 system closely resembles the HAT-P-7 system – both have hot ($>$ 6,000 K) stars hosting large ($>$ 1 ${\rm M_{Jup}}$) short-period planets with observed Kepler secondary eclipses indicating dayside brightness temperatures $>$ 2,700 K. As discussed in Section 1, Armstrong et al. (2016) applied a model similar to ours, searching for changes in the planet’s phase curve, finding statistically significant fluctuations in $\delta$ between about $6^{\circ}$ eastward and $9^{\circ}$ westward of the substellar point over timescales of tens to hundreds of days. That study then applied a model involving atmospheric transport of aerosols and found good qualitative agreement between the model and their observational results. Our results are qualitatively similar, although with much larger apparent swings in $\delta$. In addition, we were able to detect statistically significant variations in the planet’s phase curve amplitude, while Armstrong et al. (2016) did not.

4 Discussion and Conclusions

As discussed in Parmentier et al. (2016), the phase curve offset for an eclipsing planet arises, in part, from the competition between thermal emission from the hottest regions and reflection from highly refractive cloud regions on the dayside. Although the daysides of hot Jupiters like Kepler-76b are thought to be too hot for aerosol condensation, super-rotation in their atmospheres can transport nightside aersols across the western terminator onto the dayside before they evaporate. The same super-rotation is thought to shift the hottest region on the dayside eastward of the sub-stellar point. Thus, either with homogeneous clouds or without clouds at all, a hot Jupiter’s phase curve would exhibit an eastward (toward $\delta<$ 0) shift. On the other hand, reflective clouds swept over the western terminator would drive the phase curve peak west (toward $\delta>0$). According to Parmentier et al. (2016), the degree of phase shift depends, at least in part, on the cloud composition, with some of the more refractory aerosols able to travel farther east than the less refractory ones.

Variations in the phase curve may help clear up this cloudy story. For planets where reflection from clouds contributes significantly to the phase curve, we might expect a correlation between the phase curve amplitude and offset – more clouds would produce more signal with a larger westward shift. Figure 9 shows a robust (17-$\sigma$) positive correlation between $A_{\rm planet}$ and $\delta$. Since the contributions from clouds and from the dayside hotspot are thought to drive the phase curve peak in opposite directions, their relative influences are encoded in this correlation, and future work with more advanced models for the planetary disk (e.g., Louden & Kreidberg, 2017) could help to tease them out.

By making some simplifying assumptions about the phase curve, however, we can draw a few tentative conclusions regarding Kepler-76b’s atmosphere. As discussed above, a plausible maximum equilibrium temperature expected for Kepler-76b is about 2,300 K, and therefore the maximum contribution to the eclipse is probably about 20 ppm. The other 67 ppm must arise from reflected and/or scattered light. Including Rayleigh scattering and thermal emission, the cloudless atmosphere models from Parmentier et al. (2016) suggest an effective geometric albedo in the Kepler bandpass for Kepler-76b $p_{\rm geo}\approx 0.15$, considerably less than we infer. Thus, some additional component seems to be required. Parmentier et al. (2016) also indicate that aerosols swept from the nightside and evaporated on the dayside are likely confined to near the western (night-to-day) terminator. If reflection from these aerosols made up the remaining 67 ppm signal for the eclipse, however, we might expect the phase curve to be offset to the west (i.e., $\delta>0$). This line of reasoning suggests Kepler-76b’s dayside is enshrouded by a reflective and homogeneous cloud deck, and the increase in phase curve amplitude and westward shift in offset in Figure 8 may arise from periodic arrival of additional aerosols from the nightside.

As discussed in Armstrong et al. (2016), the transport of reflective clouds to the dayside would likely perturb the dayside temperature profiles, potentially reducing the day-night temperature contrast. Since the strength of the super-rotation depends, in part, on this contrast (Showman & Polvani, 2011), a reduction in contrast might reduce both windspeeds and the degree of dayside cloud cover, setting up a feedback. Such a feedback could align with the scenario above, and detailed atmospheric models that include this feedback may help determine the timescales for this feedback for Kepler-76b. In addition to the planet’s day-night temperature contrast and rotation rates (which are relatively well constrained), the timescales presumably depend on the optical properties and distribution of the aerosols, which determine their influence on the temperature structure. Thus, comparing the expected to the observed fluctuation timescales could constrain the aerosol properties and mixing in the planet’s atmosphere.

Although the full promise of Kepler observations for studying phase curve variations remains unfulfilled, the advent of the TESS Mission portends even more insight into planetary atmospheres. Already, TESS  is detecting phase curves and secondary eclipses from hot Jupiters (Shporer et al., 2018), and more are expected. For example, Figure 10 shows estimates of signal-to-noise ratios (SNR) expected for secondary eclipses observed by TESS based on the synthetic population of planets in the TESS yield calculations from Barclay et al. (2018). For each synthetic planet from that study, we estimated the eclipse depth $D$ for (1) the case in which the planet reflects all of the light it receives from the star (shown in blue) and (2) the case in which the planet emits (and absorbs) light as a perfect blackbody (shown in orange). For the former case, the eclipse depth is given as $\frac{1}{4}\left(\frac{R_{\rm p}}{a}\right)^{2}$. We then estimated the corresponding SNR by multiplying the transit SNRs given for each planet by planet-star radius ratio squared and then by $D_{\rm reflected}$. For the reflected light eclipse depths, we used the stellar insolation given for each planet and assumed dayside thermal equilibrium with zero albedo for each planet to estimate a brightness temperature and convolved the resulting blackbody curve against the TESS spectral response function222https://tessgi.github.io/TessGiWebsite/the-tess-space-telescope.html. We combined the insolation and planet’s orbital distance to estimate the stellar effective temperature and, likewise, convolved the corresponding blackbody curve against the response function. Finally, we divided the planet’s convolved brightness by the star’s and used the re-scaled transit SNR to calculate the SNR for the thermally emitted eclipse. Figure 10 only shows those SNR-values $>3$, among which the vast majority (54 of 62 total) have $R_{\rm p}\geq 10\ R_{\rm Earth}$.

The synthetic population from Barclay et al. (2018) included 4,553 planets, and Figure 10 suggests only about 1% of these may have secondary eclipses detectable by TESS. This result differs from a similar analysis in Sullivan et al. (2015), which suggested only 0.01% of TESS planets (about 2) will exhibit observable eclipses, with the difference due to the larger population of short-period, giant planets included in the population we used (Barclay et al., 2018). Yet, among those with observable eclipses, our calculation suggests that the eclipse signals will be considerably more sensitive to emitted rather than reflected light, in contrast to the eclipses observed by Kepler; this is largely due to the fact that the TESS spectral response function is more sensitive at longer wavelengths than Kepler. Presumably, this different sensitivity means secondary eclipses detectable by TESS will reveal atmospheric thermal structures in new ways. However, TESS will observe many targets for fewer than 30 days, meaning the kind of variability study presented here would be difficult for all but the most pronounced planetary phase curves or those near the continuous viewing area. The planets with the more easily detectable phase curves, hot Jupiters, also produce the deepest transit signals and are most easily studied from the ground. Barclay et al. (2018) estimated that a significant fraction of TESS discoveries will be these largest planets, many of which will lie within TESS’s continuous viewing zone, providing nearly a year of monitoring. For a planet like Kepler-76b, with an orbital period of 1.5 days, these observations would span hundreds of orbits. Thus, in addition to an exciting menagerie of Earths and super-Earths, TESS is poised to reveal a large population of hot Jupiters, potentially ripe for atmospheric characterization.

This study is based on work supported by NASA under Grant no. NNX15AB78G issued through the Astrophysical Data Analysis Program by the Science Mission Directorate. The authors gratefully acknowledge assistance from Loretta Cannon in crafting this manuscript and helpful input from Nick Cowan, Lisa Dang, Tamara Rogers, Alexander Hindle, and an anonymous referee. Kepler.