Evidence for Atmospheric Cold-trap Processes in the Noninverted Emission Spectrum of Kepler-13Ab Using HST/WFC3

We began our image calibration from the flt images provided by STScI. We first flat-fielded the images using the procedure outline in Section 6.2 of the aXe handbook (Kümmel et al. 2011). This necessitated an initial wavelength solution for our images, which we determined using the method described in Section 2.2, but here using the flt images. We then used this wavelength solution and the flat-field coefficients given in the G141 flat-field data cube provided by STScI to determine, and apply, a flat-field correction over all of our grism images.

We next defined a bad-pixel mask by manually selecting bad pixels on the flat-field image, and we corrected these bad pixels on the flat-fielded images by replacing their values with the median values of the eight pixels immediately around the bad pixels. Ultimately, none of the bad pixels lay within our spectroscopic extraction aperture, rendering this step largely unnecessary.

To identify and remove cosmic-ray hits, we divided each image into two sections: the area around the stellar spectra, and the surrounding area dominated by the sky background. Within the area around the stellar spectra, which we defined to be 140 pixels wide and 50 pixels tall and centered on the two spectral traces in the first grism image, we began our cosmic-ray rejection by creating a model image by median-combining the area in all of our flat-fielded and bad-pixel-corrected images. We then divided each individual image by this model and identified cosmic-ray hits as any time a pixel was more than 2.5 times higher than the median-combined model image. We replaced the pixels in the affected images with the value of that pixel in the model image.

For the lower-signal-to-noise background-dominated section of the image, we again created a median-combined model image from all of our flat-fielded and bad-pixel-corrected images. Instead of dividing each exposure by this model image, here we subtracted it—after fitting for a scaling factor to account for the changing pedestal value of the background. Taking the standard deviation of each median-subtracted exposure, we identified cosmic-ray hits as pixels that lay more than four standard deviations above zero. Again, we replaced the pixels in the affected images with the value of that pixel in the model image.

Finally, we calculated and subtracted the background from each of our grism exposures. To do so, we defined two background regions across the bottom (y-pixels 5–45) and top (y-pixels 240–255) of each of the flat-fielded, bad-pixel- and cosmic-ray-corrected images. Using these two background stripes, we fit for the background in the central portion of each image assuming that the background varied as a 2D plane, and we subtracted this plane fit from each image. We also experimented with estimating the background using the master sky images provided by the Space Telescope Science Institute and found that our 2D plane background estimate was within 5% of the best-fit background estimated using the master images. Since the median background level was 0.52 electrons per second per pixel for the April visit and 0.40 electrons per second per pixel for the October visit, the difference in the background estimate is negligible relative to the average target counts of approximately 1500 electrons per second per pixel.

This process left us with flat-fielded, bad-pixel- and cosmic-ray-corrected, background-subtracted images with which we performed our spectral extraction.

We used the direct image taken at the beginning of each of the visits to establish an initial wavelength solution for that visit. To do so, we used the daofind routine implemented in the PhotUtils Python package to determine the x- and y-pixel locations of Kepler-13A and Kepler-13BC on the detector subarray. We then used the wavelength calibration method described by Kuntschner et al. (2009), with the adjusted wavelength coefficients determined by Wilkins et al. (2014), to calculate a wavelength solution for each star. In the spectra of both Kepler-13A and Kepler-13BC the Paschen–β line at 1.282 μm was clearly visible, and we verified the accuracy of our initial wavelength solution using this feature in the first grism exposure of each visit.

We assumed that this initial wavelength solution was accurate for the first grism exposures in each of the visits, which were taken immediately after the direct image. For each subsequent grism exposure in a visit, we used the spectral alignment method described by Wilkins et al. (2014) to estimate a wavelength shift in the spectrum along the dispersion direction. Briefly, the Wilkins et al. (2014) alignment method median-combines all of the unshifted out-of-eclipse 1D spectra from a visit and uses this as a reference spectrum against which one measures the wavelength shift of an individual exposure's spectrum. We did an initial extraction of all our 1D spectra using the extraction method described in Section 3 and thus computed a master spectrum for each of the two visits. We cross-correlated this master spectrum against each exposure in the visit using a range of shifts from −10 pixels to +10 pixels. We linearly interpolated the master spectrum and fit for a normalization constant at each of the cross-correlation steps. We assumed that the dispersion of the wavelength solution did not change over the course of a visit.

This procedure gave us wavelength shifts for each individual exposure relative to the median-combined master spectrum. Based on our assumption that the wavelength solution derived from the direct image is correct for the first grism exposure in a visit, we subtracted the wavelength shift we measured for this first exposure from our cross-correlated master shifts to determine the shift of each individual exposure relative to our initial wavelength solution (Figure 1).

Zoom In Zoom Out Reset image size

As described in the Introduction, both Kepler-13A and Kepler-13B are mid- to late A dwarfs of nearly equal brightness. Santerne et al. (2012) identified a third star in the Kepler-13 system, which they determined to be on a 65.831-day orbit about Kepler-13B and which they estimated contributes about 1% of the combined system light.

Initial spectroscopic characterization of Kepler-13A and Kepler-13B by Szabó et al. (2011) gave estimated effective temperatures for the two stars of 8200–8500 K, masses around 2 ${M}_{\odot }$, and radii around 2.5 ${R}_{\odot }$. Shporer et al. (2014) subsequently collected an independent set of spectral observations, and their analysis gave effective temperatures of around 7500 K, masses near 1.7 ${M}_{\odot }$, and radii of about 1.7 ${R}_{\odot }$. Concurrently with the measurements of Shporer et al. (2014), Huber et al. (2014) performed a bulk reanalysis of all the Kepler target stars and found that Kepler-13A had a temperature of 9100 K, a mass of 2.5 ${M}_{\odot }$, and a radius of 3 ${R}_{\odot }$.

The lack of agreement on the temperature, mass, and radius of Kepler-13A has complicated previous studies of the planetary system (e.g., Esteves et al. 2015). To distinguish between these three sets of stellar parameters, we compared the stellar density implied by their masses and radii for Kepler-13A to the stellar density measured for Kepler-13A using the Kepler transit light curve (Seager & Mallén-Ornelas 2003). The Szabó et al. (2011) properties give a density of 0.17 g cm⁻³ (they gave no uncertainties on their mass and radius estimates), the Huber et al. (2014) properties give 0.12 ± 0.08 g cm⁻³, and the Shporer et al. (2014) properties give a density of 0.49 ± 0.07 g cm⁻³. The Kepler light curve, as measured by Shporer et al. (2014), gives a density of 0.52 ± 0.03 g cm⁻³.

Since the Shporer et al. (2014) stellar properties are the only set that correctly reproduce the stellar density measured via the transit light curve, we adopted the Shporer et al. (2014) stellar properties as correct and do not consider either the Szabó et al. (2011) or Huber et al. (2014) estimates in our reduction or analysis. We list the aggregate stellar properties for all three stars in Table 1.

Table 1.  Properties of Kepler-13A, Kepler-13B, and Kepler-13C

Parameter	Value	References
Kepler-13A
Teff (K)	7650 ± 250	Sh14
$\mathrm{log}(g)$	4.2 ± 0.5	Sh14
$[\mathrm{Fe}/{\rm{H}}]$	0.2 ± 0.2	Sh14
v $\sin (i)$ (km s−1)	76.96 ± 0.61	Jo14
M* (${M}_{\odot }$)	1.72 ± 0.10	Sh14
R* (${R}_{\odot }$)	1.71 ± 0.04	Sh14
Kepler-13B
Teff (K)	7530 ± 250	Sh14
$\mathrm{log}(g)$	4.2 ± 0.5	Sh14
$[\mathrm{Fe}/{\rm{H}}]$	0.2 ± 0.2	Sh14
v $\sin (i)$ (km s−1)	63.21 ± 1.00	Jo14
M* (${M}_{\odot }$)	1.68 ± 0.10	Sh14
R* (${R}_{\odot }$)	1.68 ± 0.04	Sh14
Kepler-13C (Orbiting Kepler-13B)
P (days)	65.831 ± 0.029	Sa12
e	0.52 ± 0.02	Sa12
M* (${M}_{\odot }$)	0.40 M⊙–0.75 M⊙	Jo14

References. Sa12 = Santerne et al. 2012; Sh14 = Shporer et al. 2014; Jo14 = Johnson et al. 2014.

Download table as:  ASCIITypeset image

As mentioned previously, Kepler-13A is separated from Kepler-13BC by $1\buildrel{\prime\prime}\over{.} 15$ (Shporer et al. 2014), or approximately 8 pixels on the WFC3/IR detector. Although the similar spectral types and brightness of Kepler-13A and Kepler-13B led us to expect that there would be relatively little wavelength-dependent dilution in our eclipse observations, we still wished to rigorously account for the effect of the added light from the wings of Kepler-13BC's spectrum.

To do so, we used the wayne simulator (Varley et al. 2017) to generate artificial 2D spectra of both stellar components, using the stellar properties listed in Table 1. We then subtracted the simulation for one member of the stellar system (e.g., Kepler-13BC) from our observed images to create an undiluted 2D spectrum for the remaining member (e.g., Kepler-13A).

wayne is a simulator for HST/WFC3 spectroscopic images, observed with any of the two infrared grisms (G102 and G141) in both observing modes (staring and spatial scanning). It is able to simulate a number of detector characteristics, such as the read noise, the nonlinearity effect, the dark current, and the wavelength-dependent flat field, as well as positioning issues, such as horizontal and vertical shifts and scan speed variations (Varley et al. 2017, and references within). A key feature in wayne is the field-dependent structure of the spectrum, as it takes into account the changes in the spectrum trace and the wavelength solution when the spectrum moves on the detector (Kuntschner et al. 2009). In addition, to simulate the wavelength-dependent point-spread function (PSF), wayne uses a linear combination of a pair of 2D Gaussian distributions, which results in point sources consistent with the PSF ensquared energy fraction given in the Wide Field Camera 3 Instrument Handbook (Dressel 2016).

For our simulation of the Kepler-13 system, we created two independent simulations, one for each component, and then combined them to a final data set. The position of the simulated spectra was based on the direct (nondispersed) image of the target and the horizontal and vertical shifts that occurred during the observations (Figure 1). In this simulation we included the photon noise, the read noise, the wavelength-dependent flat field, and the sky background, to simulate the flt images. At the stage of the flt images, the remaining reduction steps have already been applied, and for this reason we did not include them. In another set of simulations we did not combine the two companions, and also we did not include any source of noise or detector characteristics, in order to use them as models for subtracting each companion from the original frames.

As inputs for the simulations we used two synthetic spectra for Kepler-13A and Kepler-13BC, rotationally broadened to the measured $v\sin i$ of the stars (Table 1). To create each exposure, the simulator multiplies them by the sensitivity curve of the G141 grism to calculate the expected electron rate as a function of wavelength (Kuntschner et al. 2011). Furthermore, the spectra are scaled based on a combined model for the transit and the systematics (r_a, r_b1, r_b2 in Varley et al. 2017, Equation (8)). Finally, the electrons per wavelength channel are distributed on the detector by randomly sampling the wavelength-dependent PSF described above. However, since this scaling is only approximate, in order to subtract the simulated model spectra from the real ones, we fit for a scaling factor between them. The parameters used for the models are a combination of the stellar parameters listed in Table 1 and the planetary parameters listed in Table 2.

Table 2.  Prior Values for Kepler-13Ab's Properties from Shporer et al. (2014)

Parameter	Units	Value
TSa	Predicted TS (${\mathrm{BJD}}_{\mathrm{TDB}}$)	2,456,776.23411 ± 0.00008
P	Orbital period (days)	1.76358799 ± 3.7 × 10−7
$e\cos \omega$ b	⋯	−0.00015 ± 0.00004
$e\sin \omega$ b	⋯	0.0 ± 0.00005
$\cos i$	Cosine of inclination	0.0714 ± 0.008
${R}_{P}/{R}_{* }$	Radius ratio	0.0845 ± 0.0012
$a/{R}_{* }$	Scaled semimajor axis	4.4 ± 0.16
MPc	Planet mass ($\,{M}_{{\rm{J}}}$)	6.52 ± 1.58
RPc	Planet radius (RJ)	1.406 ± 0.038

Notes.

^aFor the April visit. During fitting, we calculate the eclipse time for the October visit using this T_S and advancing by 95 times the orbital period, as described in Section 3.2. ^bEstimated using the eclipse time and duration given in Shporer et al. (2014) and the first-order approximations for $e\cos \omega$ and $e\sin \omega$ given in Winn (2010). ^cNot a fitting parameter, but provided for reference.

Download table as:  ASCIITypeset image

As usual with HST/WFC3 grism observations, our raw extracted photometry displayed clear correlated noise (top panel of Figure 2). As delineated by Wakeford et al. (2016), these are typically a long-term temporal trend over the course of several orbits, an "L"-shaped hook during the course of an individual orbit, and quick variations within an orbit due to the spacecraft's thermal breathing. Only the first (visit-long slope) and third (thermal breathing) effects were clearly present in our broadband and spectrally resolved photometry of Kepler-13A and Kepler-13BC.

Zoom In Zoom Out Reset image size

The raw photometry did show a downward trend during the course of each orbit's observations (top panel, Figure 2). This is not usually seen in WFC3 observations, but this observation is one of the very few where the SPARS10 sampling sequence at NSAMP = 3 was used. In the similar case of WASP-12b (Mandell et al. 2013; Swain et al. 2013) a similar behavior can be seen. This may be caused when the detector is continuously flushed to prevent charge buildup after an exposure (Deustua 2016, p. 158), in combination with the large number of exposures per orbit. Another possibility is that this trend was introduced by the subtraction process we just described. To test for this, we extracted photometry for Kepler-13A without subtracting off our model for Kepler-13B. As shown in Figure 3, while the no-subtraction photometry shows higher scatter than the subtracted photometry in Figure 2, the orbit-long downward trend is still present, and thus not a result of our subtraction process.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At this point the presence of Kepler-13BC greatly simplified our measurement of Kepler-13Ab's eclipse, by providing a nearly equal brightness and equal color comparison light curve with which to perform differential photometry. Since Kepler-13C contributes a negligible amount of light to the Kepler-13BC light curve, the broadband response of Kepler-13BC to the HST/WFC3 systematics closely matched the response of Kepler-13A.

We initially attempted to create a differential light curve for Kepler-13A by the straight division of the Kepler-13BC light curve, similar to how we would treat this in ground-based photometry. To our initial surprise, this gave a differential light curve for Kepler-13A that displayed a much higher scatter than the raw light curve, about 700 ppm versus about 450 ppm, even though the systematic trends seemed to be mostly gone. Upon further consideration, we realized that this was because our raw HST/WFC3 photometry—unlike ground-based photometry—was dominated by photon, rather than systematic, noise. Thus, while dividing the Kepler-13A by the Kepler-13BC light curve removed the systematic trends, it also added the photon noise scatter of the Kepler-13BC in quadrature to the photon noise scatter in Kepler-13A.

We therefore applied a low-pass filter to the Kepler-13BC light curve, to capture the systematic trends in that system's photometry while removing most of the photon noise. Specifically, we used a second-order Butterworth filter with a cutoff at 0.02 days. As shown by the red line in the middle panel of Figure 2, the thus-filtered Kepler-13BC light curve appeared to correctly capture the trends that we were interested in. We then generated our differential light curve for the white-light Kepler-13A data by dividing the raw Kepler-13A photometry by this filtered Kepler-13BC light curve. As shown in the bottom panel of Figure 2, this removed most of the systematics and made the broadband eclipse clearly visible during both visits.

Another advantage of using Kepler-13BC to make a differential light curve for Kepler-13A is that we were able to make effective use of the first orbit of observations in both of our visits. Typically the entire first orbit of a visit is discarded for precise eclipse or transit observations, due to the increased ramp and thermal breathing affects. In our case, it is a relatively simple matter to remove most of these affects using our differential comparison. The use of the first orbit's data was particularly helpful in the April visit, as without the first orbit we would not have had any pre-eclipse baseline observations.

For the spectrally resolved data we also created differential light curves, but instead of differencing against Kepler-13BC, we differenced against Kepler-13A itself. We were motivated to do so out of a concern that the different placement of the two stars on the detector, plus their slightly different intrinsic spectra, would generate spurious spectral signatures if we differenced against Kepler-13BC. We applied the same low-pass filter as we used in the broadband data to the broadband Kepler-13A photometry, and we divided each of our spectrally resolved light curves by this filtered broadband photometry.

In our spectrally resolved data we therefore did not measure the absolute eclipse of Kepler-13Ab, but rather the change in eclipse depth as a function of wavelength.

In both the broadband and spectrally resolved differential light curves there was residual systematic noise. We chose to fit these residual correlations using a GP regression model, as the traditional parametric methods for dealing with HST/WFC3 systematics were not applicable to our differential light curves.

A GP regression models the observed data as random draws from a multivariate Gaussian distribution about some mean function. As a result, GPs are able to directly model the possible covariances between data points by populating the nondiagonal elements of the covariance matrix that defines the multivariate Gaussian distribution. This is in contrast, for example, to a ${\chi }^{2}$ fitting process, which models the data as random draws from a univariate Gaussian distribution and assumes no covariance between data points—aside from what is inserted via detrending functions. For more detail, Rasmussen & Williams (2006) provide a thorough mathematical overview of GP methods. Gibson et al. (2012) introduced them in the context of astronomical time series observations using archival NICMOS observations of HD 189733b, though they have a longer history in the general astronomical community (e.g., Way & Srivastava 2006). Recently, Cartier et al. (2017) used a GP regression to model HST/WFC3 eclipse observations of WASP-103b, and we follow a similar approach.

We defined our GP model using the notation from Gibson et al. (2012). For each visit, we have a vector of N observed fluxes, $f=({f}_{1},\ldots ,\,{f}_{N})$, and times, $t=({t}_{1},\ldots ,\,{t}_{N})$. Additionally, we recorded K state parameters at each time t with the state vector $x={({x}_{t,1},\ldots ,{x}_{t,K})}^{T}$. We combined these state parameter vectors for each of our N observations in the N × K matrix, X. The multivariate Gaussian distribution underlying our GP model was defined by a combination of a mean function, which in our case is the physical eclipse model $E(t,\phi )$, and a covariance matrix ${\rm{\Sigma }}(X,\theta )$. We used ϕ to denote the set of physical parameters describing the eclipse and θ for the set of "hyperparameters" used to generate the covariance matrix from the X state parameters. The joint probability distribution of our observed data f was then

$$\begin{eqnarray}&&p(f| X,\theta ,\phi )={ \mathcal N }[E(t,\phi ),{\rm{\Sigma }}(X,\theta )],\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal N }$ represents the multivariate Gaussian distribution. Our GP model thus depended on on an eclipse model $E(t,\phi )$ and a generating function—referred to as the covariance kernel—for the covariance matrix ${\rm{\Sigma }}(X,\theta )$. We generated our GP covariance matrices and calculated the GP likelihoods using the George python package (Ambikasaran et al. 2014).

While we used the same eclipse model to fit both the broadband and spectrally resolved data sets, we chose slightly different covariance kernels for each set of observations.

3.1.1. The Eclipse Model, Parameters, and Priors

We used a Mandel & Agol (2002) eclipse model as the mean function in our GP regression, specifically the implementation in the BATMAN python package (Kreidberg 2015). For our broadband data, we fit for the time of the secondary eclipse (T_S), the orbital period (as $\mathrm{log}[P]$), $e\cos \omega$, $e\sin \omega$, the cosine of the orbital inclination ($\cos i$), the radius of the planet in stellar radii (${R}_{P}/{R}_{* }$), the semimajor axis in units of the stellar radii (as $\mathrm{log}[a/{R}_{* }]$), and the depth of the secondary eclipse (δ). This gave our broadband eclipse model eight parameters:

$$\begin{eqnarray}{\phi }_{\mathrm{white}} & = & ({T}_{S},\mathrm{log}P,e\cos \omega ,e\sin \omega ,\cos i,\\ & & {R}_{P}/{R}_{* },\mathrm{log}a/{R}_{* },\delta ).\end{eqnarray} \tag{ 2 }$$

Based on Shporer et al. (2014)s measurements of the system, we have strong prior expectations for all of these parameters, except for the eclipse depth δ (Table 2). In our fitting, we implement each of these as Gaussian priors on these seven parameters. We do not impose any prior on the eclipse depth, as we lack any observations of Kepler-13Ab's eclipse at these wavelengths, which implicitly imposes a uniform prior.

For the spectrally resolved observations the implementation of our eclipse model remains the same, but we only fit for the time of the secondary eclipse and the eclipse depth. All of the other parameters we fix to the values we determine from the broadband fit. This gives our spectrally resolved eclipse model only two parameters:

$$\begin{eqnarray}&&{\phi }_{\mathrm{spec}}=({T}_{S},\delta ).\end{eqnarray} \tag{ 3 }$$

Again, we impose no prior on the eclipse depths, and we use the measured time of eclipse and its associated uncertainty from the broadband fit as a prior on T_S for the spectrally resolved fits.

3.1.2. GP Covariance Kernel and Hyperparameter Priors

To model our broadband observations, we used a linear combination of a squared-exponential kernel and a periodic kernel—both as a function of time. A squared-exponential kernel is usually regarded as a generic choice for generating a GP covariance matrix, as it leads to smooth variations as a function of the generating variable. Our intent was to use it to model the residual background temporal trends. We added on a periodic kernel, which is effectively a squared exponential of a sine function, to model the repeatable covariances between points within an individual orbit.

The point-wise covariances between the observations at times t_i and t_j, which collectively made up our broadband N × N covariance matrix ${{\rm{\Sigma }}}_{\mathrm{white}}$, were then

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{white}}(t,{\theta }_{\mathrm{white}})={{\rm{\Sigma }}}_{\mathrm{SqExp}}(t,\theta )+{{\rm{\Sigma }}}_{\mathrm{Per}}(t,\theta ),\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{SqExp}}(t,{\theta }_{\mathrm{white}})={A}_{t}\,\exp \left[-\displaystyle \frac{{({t}_{i}-{t}_{j})}^{2}}{{L}_{t}^{2}}\right]\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{Per}}(t,{\theta }_{\mathrm{white}})={A}_{p}\,\exp \left[\displaystyle \frac{{\sin }^{2}(\pi [{t}_{i}-{t}_{j}]/p)}{{L}_{p}^{2}}\right].\end{eqnarray} \tag{ 6 }$$

We used ${\theta }_{\mathrm{white}}=\{{A}_{t},{L}_{t},{A}_{p},{L}_{p},p\}$ to denote the hyperparameters used to compute the covariances. These are covariance amplitudes A_t and A_p and covariance length scales L_t and L_p. We set the period of the periodic kernel to be $p\equiv 95.664$ minutes. This was the measured orbital period of HST on the dates of our observations, based on archival two-line elements provided to us by the United States Strategic Command's Joint Space Operations Center.

For our spectrally resolved data, we used a slightly simplified covariance kernel with only a periodic component,

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{spec}}(t,{\theta }_{\mathrm{spec}})={A}_{p}\,\exp \left[\displaystyle \frac{{\sin }^{2}(\pi [{t}_{i}-{t}_{j}]/p)}{{L}_{p}^{2}}\right],\end{eqnarray} \tag{ 7 }$$

where we have ${\theta }_{\mathrm{spec}}=\{{A}_{p},{L}_{p},p\}$ hyperparameters for the spectrally resolved data, with the same definitions as in the broadband case.

We imposed no priors on the members of either set of hyperparameters, other than to require that both L_t and L_p be greater than the eclipse's ingress and egress time of 0.01279 days. This ensured that the GP regression model did not treat the eclipse itself as correlated noise.

To determine the best fits to the broadband and spectrally resolved data, we maximized a log-likelihood function that consisted of the GP model likelihood and an additional term based on the priors in Table 2. For the GP model the log-likelihood for a given set of parameters ϕ and θ was

$$\begin{eqnarray}&&\mathrm{log}{p}_{\mathrm{GP}}(r| X,\theta ,\phi )=-\displaystyle \frac{1}{2}{{\boldsymbol{r}}}^{T}{{\rm{\Sigma }}}^{-1}{\boldsymbol{r}}-\displaystyle \frac{1}{2}\mathrm{log}| {\rm{\Sigma }}| -\displaystyle \frac{N}{2}\mathrm{log}(2\pi ),\end{eqnarray} \tag{ 8 }$$

where ${\boldsymbol{r}}=f-E(t,\phi )$ is the vector of the residuals of our observed data (f) from our eclipse model (E) defined in Section 3.1.1. This log-likelihood follows directly from our definition of the GP model as a multivariate Gaussian in Equation (1).

For all of the physical parameters for which we applied a prior based on previous observations of the system we used Gaussian priors, and we added the log-likelihood of these priors to our GP log-likelihood to compute the total log-likelihood for a given model as

$$\begin{eqnarray}&&\mathrm{log}{p}_{\mathrm{tot}}(\phi ,\theta | f,X)=\mathrm{log}{p}_{\mathrm{GP}}(r| X,\theta ,\phi )+\sum \mathrm{log}{p}_{\mathrm{prior}}.\end{eqnarray} \tag{ 9 }$$

Our fitting process was composed of two stages for both the broadband and the spectrally resolved data: an initial Nelder–Mead maximization of Equation (9), followed by an MCMC exploration of the likelihood space around this maximum to determine parameter uncertainties and to verify that we had identified the global likelihood maximum. To conduct the MCMC fitting, we used the emcee Python package (Foreman-Mackey et al. 2013). Our MCMC runs consisted of a 500-step burn-in, followed by a 5000-step production run using 30 walkers. We initialized the walkers by scattering them about the initial Nelder–Mead maximum using random draws from Gaussian distributions in each parameter with $1\sigma$ widths equal to corresponding prior widths. At the end of the production runs we calculated the Gelman–Rubin statistic and judged the MCMC to have converged if the statistics for all the parameters were less than 1.1.

We fit the two visits simultaneously using the same set of physical parameters and hyperparameters. We experimented with using differing hyperparameters for each of the visits, but we found that this did not change the results above our final $1\sigma$ uncertainties, nor did it substantially alter the uncertainties themselves. For the spectrally resolved data we did allow for different hyperparameters between different spectral channels.

Tables 3 and 4 list the results for the broadband eclipse and the spectrally resolved eclipses, respectively. Our best fit to the broadband data gave a broadband eclipse depth of $\delta =734\pm 27\,\mathrm{ppm}$ (Figure 4). The standard deviation in the residuals to the best-fit model was 436 ppm, which was nearly equal to the median per-point flux uncertainty of 400 ppm from photon statistics. The uncertainty on the broadband measurement was 1.5 times what one would expect based on pure photon noise alone.

Table 3.  Median Values and 68% Confidence Intervals for the Broadband Eclipse Observations

Parameter	Units	Value
GP Hyperparameters
At	Sqr. exp. covariance amplitude	$2.4\times {10}^{-6}{\,}_{-1.6\times {10}^{-6}}^{+6.5\times {10}^{-5}}$
Lt	Sqr. exp. covariance length scale	${1.48}_{-0.32}^{+0.92}$
Ap	Periodic covariance amplitude	$1.4\times {10}^{-7}{\,}_{-6.7\times {10}^{-8}}^{+1.8\times {10}^{-7}}$
Lp	Periodic covariance length scale	${2.38}_{-0.57}^{+0.69}$
Eclipse Model Parameters
TS	Eclipse Time (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2456776.234119}_{-0.000079}^{+0.000073}$
$\mathrm{log}(P)$	Log orbital period (days)	$0.24639714\pm 9.2\times {10}^{-8}$
$e\cos \omega$		−0.00015 ± 0.00004
$e\sin \omega$		0.0 ± 0.00005
$\cos i$	Cosine of inclination	0.069 ± 0.008
${R}_{P}/{R}_{* }$	Radius of planet in stellar radii	0.0844 ± 0.0012
$\mathrm{log}(a/{R}_{* })$	Log semimajor axis in stellar radii	0.6325 ± 0.0075
δ	Eclipse depth (ppm)	734 ± 28
Derived Parameters
P	Orbital period (days)	1.76358803 ± 0.00000037
$a/{R}_{* }$	Semimajor axis in stellar radii	4.29 ± 0.08
i	Inclination (deg)	86.04 ± 0.44
b	Impact parameter	0.296 ± 0.031
TFWHM	FWHM duration (days)	${0.1263}_{-0.0022}^{+0.0025}$
τ	Ingress/egress duration (days)	0.0119 ± 0.0003
T14	Total duration (days)	${0.1382}_{-0.0024}^{+0.0027}$
e	Orbital eccentricity	0.00016 ± 0.00004
ω	Argument of periastron (deg)	${128}_{-300}^{+44}$

Download table as:  ASCIITypeset image

The differential depths we measured in our spectrally resolved light curves (Figures 5 and 6) had a median uncertainty of 64 ppm, which is 1.2 times the photon noise expectation. The median standard deviation of the residuals to the best-fit models for the spectrally resolved data was approximately 1250 ppm. The differential eclipse depths showed statistically significant variation away from a flat line, with ${\rm{\Delta }}{\chi }^{2}=64.3$, or ${\rm{\Delta }}{\chi }^{2}/\mathrm{dof}=4.6$ with 14 degrees of freedom. This corresponds to a $5.6\sigma$ detection of variation with wavelength.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a verification of our results, we also performed the exact same fitting procedures, using the exact same priors, on broadband and spectrally resolved light curves from Kepler-13BC. As expected, the broadband Kepler-13BC light curve gave an eclipse depth of ${\delta }_{\mathrm{BC}}=-10\pm 34\,\mathrm{ppm}$. As shown in Figure 5, the differential spectrally resolved eclipse depths for Kepler-13BC are also all consistent with zero, with ${\rm{\Delta }}{\chi }^{2}/\mathrm{dof}=0.3$ with 14 degrees of freedom.

We modeled the measured eclipse depth of Kepler-13Ab using the hot Jupiter atmospheric model described in Madhusudhan & Seager (2009) and Madhusudhan (2012). This model is composed of a 1D plane-parallel atmosphere in hydrostatic equilibrium and local thermodynamic equilibrium. The emergent spectrum is computed using a 1D line-by-line radiative transfer solver in the planetary atmosphere. The atmospheric pressure–temperature (PT) profile and chemical composition are free parameters of the model, with six parameters for the PT profile and a parameter each for each chemical species included in the model. The range of molecules considered and the sources of opacity are discussed in Madhusudhan (2012). For these models we set the mixing ratios of the major molecular constituents (${{\rm{H}}}_{2}{\rm{O}}$, ${\mathrm{CH}}_{4}$, ${\mathrm{CO}}_{2}$, CO) to the chemical equilibrium solar values. The inverted temperature profile in Figure 7 used solar TiO and VO abundances, while the decreasing temperature profile set the TiO and Vo abundances to zero.

Zoom In Zoom Out Reset image size

By way of comparison, the model spectra we use have the same qualitative features as the Burrows et al. (2008) models shown in Figure 11 of Shporer et al. (2014). The primary modeling difference between our results and those of Shporer et al. (2014) is that our WFC3 measurements give a higher average dayside temperature, leading to higher temperatures in our model atmospheres. Compared to the Fortney et al. (2008) models from Shporer et al. (2014), our WFC3 spectrum is best fit by the "no TiO" case, but we find a much more pronounced water feature at 1.4 μm than predicted by the Fortney models.

We find that our observations are relatively well reproduced by an atmospheric model with a monotonically decreasing temperature profile (Figure 7), with ${\chi }^{2}=29.31$ (${\chi }^{2}/\mathrm{dof}\,=1.72$, or $2.15\sigma$). We exclude the best-fit isothermal model atmosphere, with a temperature of 3000 K, with ${\chi }^{2}=46.98$ (${\chi }^{2}/\mathrm{dof}=2.76$, or $3.84\sigma$), and we exclude an inverted temperature profile ${\chi }^{2}=133.45$ (${\chi }^{2}/\mathrm{dof}=7.85$, or $10.9\sigma$). We also tested a range of monotonically decreasing temperature profiles that stepped toward and away from a purely isothermal profile (Figure 8). The green-blue profile closest to isothermal in Figure 8 is only marginally consistent with the observations ($2.96\sigma$) and is the upper limit allowed by the data.

It is important to note that our WFC3 observations and the observations of Shporer et al. (2014) only probe pressure levels in the planetary atmosphere between about ∼ 10⁻² and ∼ 10⁰ bars. Based solely on these observations, it is therefore possible that a temperature inversion could exist on Kepler-13Ab's dayside higher up and at lower pressure levels.

Due to the relatively high mass of Kepler-13Ab, we were interested to compare its dayside emission spectrum to a brown dwarf of similar effective temperature. Besides investigating the general spectral differences between the dayside of a hot Jupiter and a brown dwarf, we also wished to see whether the measured surface gravity of Kepler-13Ab coincided with the surface gravity one would estimate using spectral matching to brown dwarf templates, or using spectral indicators for surface gravity.

Compared to its isolated brown dwarf equivalent, the dayside of Kepler-13Ab is, needless to say, considerably hotter. From its age and mass we would expect that Kepler-13Ab would have a surface effective temperature of approximately 400–800 K if it did not receive any insolation from its host star (Burrows et al. 2003). But due to the propinquity of Kepler-13A, the dayside of the planet appears to be closer in temperature to a mid- to late M dwarf. We therefore compared our emission spectrum to a series of spectral templates from the SpeX Library covering M4 to T9 spectra for field objects and a set of low-gravity templates covering M5 to T9 from Allers & Liu (2013). To do so, we converted the normalized flux measurements in the templates to a set of simulated eclipse depths by assuming that Kepler-13A was a 7650 K blackbody and then scaling the resulting relative flux based on the estimated J-band surface brightnesses for Kepler-13A and the spectral templates from the BT-Settl models (Allard et al. 2011) using the Caffau et al. (2011) values for solar abundances.

We find that the dayside of Kepler-13Ab has a rough spectral type of M7 ± 2 using both the field and low-gravity templates. Using only our 15 WFC3 measurements, we find ${\chi }^{2}=21.19$ (${\chi }^{2}/\mathrm{dof}=1.51$, or $1.66\sigma$) for the field template and ${\chi }^{2}=20.40$ (${\chi }^{2}/\mathrm{dof}=1.70$, or $1.88\sigma$) for the low-gravity template. As shown in Figure 9, in both cases this result is driven by the presence of the 1.4 μm ${{\rm{H}}}_{2}{\rm{O}}$ absorption feature. We note, though, that the apparent improvement on the goodness of fit of these template spectra relative to the atmosphere models in Section 4.1 is due to the fact that here we are only using the WFC3 in our comparisons, due to a lack of wavelength coverage in the spectral templates.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The similarity of the higher-gravity field templates and the low-gravity templates in this temperature range means that we are not able to postdict Kepler-13Ab's surface gravity using either template matching or the surface gravity indicators suggested by Allers & Liu (2013). Nevertheless, the general agreement of both the field and low-gravity templates indicates that despite the heavy insolation Kepler-13Ab receives, its atmosphere appears similar to isolated brown dwarf atmospheres. Since isolated brown dwarfs are primarily heated from within and have correspondingly monotonically decreasing atmospheric temperature–pressure profiles, this strengthens our conclusion from the atmospheric models that Kepler-13Ab's dayside neither is isothermal nor possesses a stratospheric temperature inversion at any pressure level.

Shporer et al. (2014) measured a center time for Kepler-13Ab's eclipse that was very close to being half of an orbital period away from the transit center time, with a displacement of −2.6 ± 7.5 s. As they noted, due to the light-travel time across the diameter of Kepler-13Ab's orbit, for a perfectly circular orbit one would expect the eclipse time to actually be delayed by $+34\pm 0.7\,{\rm{s}}$. The apparent earliness of the eclipse implies a slightly eccentric orbit and is what drives the nonzero value of $e\cos \omega$ listed in Table 2.

As described in Section 3.1.1, we used the predicted eclipse time from Shporer et al. (2014) as a prior on our broadband fitting process. Since our observations poorly sample the predicted eclipse ingress and egress (i.e., Figure 4), this prior on T_S dominates in our fits, and we simply recover it in the posterior distribution for T_S from our MCMC fitting. In our adopted fit, we therefore find that the eclipse center time is −1.8 ± 7.5 s earlier than the transit center time plus exactly half the orbital period, which, as we expected, is effectively identical to the value measured by Shporer et al. (2014).

To determine how well our HST/WFC3 can independently constrain the eclipse center time, we conducted an additional broadband fit using no prior on the eclipse center. This gave ${T}_{S}=2456776.2310\pm 0.0025$, which is −269 ± 216 s earlier than the transit center time plus exactly half the orbital period. The large uncertainties on this timing offset make it consistent with the expectations from both Shporer et al. (2014) and the perfectly circular case, so our observations are not able to meaningfully determine the eccentricity of Kepler-13Ab's orbit.

One complication in interpreting the optical eclipse depth for Kepler-13Ab is the relatively high amount of reflected light expected from the planet. Specifically, the eclipse depth from reflected light alone should be ${\delta }_{\mathrm{ref}}={A}_{g}{({R}_{p}/a)}^{2}={A}_{g}(316)$ ppm, where A_g is the geometric albedo of the dayside at a particular wavelength. Since the Kepler-band eclipse measured by Shporer et al. (2014) is 173.7 ± 1.8 ppm, it is possible—though unlikely –that the optical eclipse is entirely due to reflected light. For our near-IR (NIR) observations, the theoretical expectation is that hot Jupiters at $\approx 3000$ K should have dayside geometric albedos of effectively zero at these wavelengths (Sudarsky et al. 2000). We therefore did not consider possible signatures of reflected light in our WFC3 spectra and excluded the Kepler eclipse point from our evaluation of our thermal emission models.

Based on those thermal emission models that assume that the WFC3 and further NIR data contain no reflected light, we calculate that Kepler-13Ab should have an eclipse 133 ppm deep in the Kepler bandpass. The difference between this prediction for the thermal emission and the observed depth of 173.7 ± 1.8 ppm implies that Kepler-13Ab has a geometric albedo of A_g = 0.12 in the Kepler bandpass. This is substantially lower than the value found by Shporer et al. (2014), A_s = 0.33. The difference occurs because our atmosphere models place the interior isotherm at a higher temperature than in the 2750 K blackbody used by Shporer et al. (2014) in their albedo calculation, which causes the planetary thermal emission to be higher in the optical.

Our measured geometric albedo and the expectations from Sudarsky et al. (2000) are thus both consistent with our assumption that the WFC3 and other NIR eclipse depths contain no substantial components from reflected light. This is in line with what one would expect from other measurements of hot Jupiters' eclipses, most of which have $0.05\lt {A}_{s}\lt 0.2$ in the Kepler bandpass (Heng & Demory 2013). Note, though, that Heng & Demory (2013) had to assume the amount of thermal emission present in the Kepler bandpass for their analysis and could not measure it as we have done here, which makes this only a general comparison. Additionally, the high stellar insolation that Kepler-13Ab receives makes it possible that Kepler-13Ab may have different reflective properties compared to the cooler planets considered by Heng & Demory (2013). One way to test this would be to observe an optical eclipse spectrum for the planet. In the case that Kepler-13Ab has a substantial optical albedo, the optical eclipse spectrum would be a combination of the reflectance spectrum (expected to increase toward the blue; e.g., Sudarsky et al. 2000) and the thermal emission spectrum (expected to decrease toward the blue), which would make the observed eclipse spectrum relatively flat.

As described in the Introduction, the lack of clear stratospheric temperature inversions for planets cooler than approximately 3000 K led Haynes et al. (2015) to recently suggest that that TiO/VO-driven inversions may only be present in the atmospheres of extremely hot giant planets. The Beatty et al. (2016) spectrally resolved observations of the H-band eclipse of the transiting brown dwarf KELT-1b (3200 K) showed a monotonically decreasing PT profile. Due to the high surface gravity of KELT-1b (22 times that of Jupiter), this led them to suggest that surface gravity also plays a role in the presence of a thermal inversion in hot Jupiters. Specifically, they argued that this was evidence for the cold-trap methods of sequestering TiO/VO described by Spiegel et al. (2009) and Parmentier et al. (2013).

Briefly, both Spiegel et al. (2009) and Parmentier et al. (2013) describe a process of TiO/VO gas particles condensing and gravitationally settling out of the upper atmosphere. Spiegel et al. (2009) envisioned this as a "vertical" cold trap, where TiO/VO gas on the dayside of a hot Jupiter randomly crosses the condensation boundary in the atmosphere, condenses, and falls into the planetary interior. Parmentier et al. (2013) considered a "day–night" cold trap, where TiO/VO molecules condense on the cooler nightside of a hot Jupiter and also settle into the interior. In both cases the efficiency of a cold trap is determined by the interplay of the rate of gravitational settling and the strength of some vertical lofting mechanism to bring TiO/VO condensates back into the upper atmosphere where they can revaporize. Both analyses found that under reasonable assumptions these two cold traps should be able to deplete a hot Jupiter's stratosphere of gas-phase TiO/VO.

In the cases of both Kepler-13Ab and KELT-1b, the dayside temperature profiles both remain too hot at depth to allow for a vertical cold trap to exist. Down to pressures of 10² bars, at no point are either of these profiles expected to cross to cooler than the TiO or VO condensation curves, which vary from approximately 2200 K at 10² bars down to approximately 1600 K at 10⁻⁵ bars. The most recent modeling of such a process was done by Parmentier et al. (2016), who found that a vertical cold trap of TiO/VO was mostly effective in planets with equilibrium temperatures less than 1900 K, which is several hundred degrees cooler than the planets we consider. It is also interesting to note that simulations by Tremblin et al. (2017) show that vertical cold traps may be impossible in hot Jupiters, due to advection in their atmospheres. In either case, if a cold-trap process is occurring in these atmospheres—and in the atmospheres of the other extremely hot giant planets—it is very probably caused by a day–night cold trap.

To first order, the rate of gravitational settling in a hot Jupiter's atmosphere will be given by the freefall timescale within that atmosphere. This will be the scale height of the atmosphere divided by the freefall terminal velocity of condensates,

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}=\displaystyle \frac{H}{{V}_{\mathrm{term}}}=\left(\displaystyle \frac{{k}_{{\rm{B}}}T}{{\mu }_{{\rm{m}}}g}\right)\,{\left(\displaystyle \frac{2{a}^{2}g({\rho }_{{\rm{p}}}-\rho )}{9\eta }\right)}^{-1}.\end{eqnarray} \tag{ 10 }$$

Here T is the atmospheric temperature, ${\mu }_{{\rm{m}}}$ is the mean molecular weight of the atmosphere, g is the gravitational acceleration, a is the particle radius, ${\rho }_{{\rm{p}}}$ is the particle density, ρ is the atmospheric density, and η is the viscosity of the gas. If we make the assumption that mean molecular weight and viscosity of hot Jupiters' atmospheres, as well as the condensate and atmospheric density, are all effectively the same, then the freefall timescale goes as

$$\begin{eqnarray}&&{\tau }_{\mathrm{ff}}\propto {{Ta}}^{-2}{g}^{-2}.\end{eqnarray} \tag{ 11 }$$

Typically, turbulent diffusion and large-scale vertical mixing are treated as the dominant vertical lofting mechanisms for condensates in the upper portion of a hot Jupiter's atmosphere, but the exact efficiency of these processes is poorly understood. By analogy to molecular diffusion, the efficiency of vertical mixing is parameterized by the effective diffusion coefficient, k_zz (Banks & Kockarts 1973). Unlike the molecular diffusion coefficient, however, the value of k_zz is the result of inherently chaotic processes that are difficult to model. While mixing length theory can be used to derive analytic estimates for k_zz (e.g., Gierasch & Conrath 1985) in the convective regions of a hot Jupiter's atmosphere, these estimates are not applicable to the upper, radiative portion of the atmosphere probed by eclipse observations. This make the precise value of k_zz difficult to predict, and typical values of k_zz used in brown dwarf modeling can cover three to four orders of magnitude (e.g., Saumon et al. 2007).

If we make the assumption that the efficiency of vertical mixing is approximately the same in all hot Jupiters' atmospheres, then for fixed values of a, the condensate size, cold traps will therefore be less efficient in hotter atmospheres with higher values of T, but dramatically more efficient as surface gravity increases.

Though more detailed modeling is necessary to precisely assess the role of cold-trap processes in Kepler-13Ab's atmosphere, let us proceed from the conclusion in Parmentier et al. (2013) that for particle sizes greater than "a few microns" a day–night cold trap is capable of depleting gas-phase TiO/VO in HD 189733b's atmosphere. Using the scaling in Equation (11), we may extrapolate from this point to other temperatures and surface gravities and compare against the other hot Jupiters with spectrally resolved NIR emission measurements.

Based on their 1.1–1.7 μm emission and the modeling in their respective papers, we categorized the other planets with spectrally resolved NIR eclipse observations as having an inverted PT profile (WASP-33b, WASP-121b), as having an isothermal profile (WASP-12b, WASP-103b), or as being ambiguous (WASP-19b, Bean et al. 2013; WASP-4b, Ranjan et al. 2014; TrES-3b, Ranjan et al. 2014; CoRoT-2b, Wilkins et al. 2014). Of these four ambiguous planets, TrES-3b and CoRoT-2b both show WFC3 spectra consistent with both an isothermal and decreasing PT profile. The temperatures for the isothermal models of both (1800 and 1780 K, respectively) are very close to the condensation temperature of TiO/VO. Depending on the exact temperature structure of their atmospheres, it is therefore possible that TiO and VO have condensed out of their upper atmospheres independently of any possible cold-trap mechanisms. For WASP-4b, the WFC3 measured spectrum agrees with both a 2000 K isothermal model and a carbon-rich, monotonically decreasing model. WASP-19b was observed from the ground by Bean et al. (2013), and their emission spectrum is not precise enough to distinguish the temperature structure of the planet. We categorized Kepler-13Ab, as well as KELT-1b, as having a monotonically decreasing PT profile. Extrapolating from HD 189733b according to Equation (11), we then plotted these eight planets alongside the limits of where we would expect cold traps to inhibit inversions for particle sizes of a = 1 μm, a = 3 μm, and a = 5 μm.

As shown in Figure 10, this set of classifications and predictions relatively well represents the planets with HST/WFC3 observations of their dayside emission. In the shaded red region, where we would expect the atmospheres to be inverted or isothermal for reasonable particle sizes, are WASP-33b, WASP-12b, and WASP-103b. Similarly, Both Kepler-13Ab and KELT-1b lie in the region where we would expect to see decreasing PT profiles assuming that the condensate particles grow larger than a = 1 μm.

Zoom In Zoom Out Reset image size

We note that Heng & Demory (2013) have also considered the efficiency of cold traps but instead choose to compare vertical terminal velocities to vertical mixing velocities to determine whether a condensate sinks into the planetary interior, without considering the atmospheric scale height. Based on numerical simulations of the vertical mixing velocity in planets cooler than 1750 K, this leads Heng & Demory (2013) to estimate that the efficiency of cold-trap processes should be roughly independent of temperature and should scale inversely with condensate size and surface gravity. This nicely contrasts with our Equations (10) and (11), which propose a different scaling. Which of these is correct—if any—would therefore provide us with good insight into the dominant processes within hot Jupiters' atmospheres, and could indicate whether vertical mixing velocities are roughly constant, as we have assumed, or they vary coherently with planetary properties, as assumed by Heng & Demory (2013).

Currently, the global maximum size for condensate particles in the atmospheres of extremely hot Jupiters with daysides near 3000 K is poorly constrained by atmospheric models. Based on the Rayleigh scattering signatures seen in transmission spectra, the maximum particle size at high altitudes along the planetary limb appears to be on the order of 0.1 μm (Wakeford & Sing 2015). But 3D models of HD 187933b's (${T}_{\mathrm{eq}}\approx 1500$ K) atmosphere by Lee et al. (2016) show a wide range of possible particle sizes as a function of longitude, latitude, and depth within the atmosphere. Parmentier et al. (2016) recently analyzed Kepler-band albedo estimates from eclipse measurements made by Esteves et al. (2015), in an effort to constrain condensate properties on planetary daysides and at higher pressures than those probed by transmission observations. Parmentier et al. (2016) found that an average condensate size of 0.1 μm was able to replicate the reflective properties of the cooler planets (${T}_{\mathrm{eq}}\approx 1600$ K), but that the two hotter planets in the Esteves et al. (2015) data set near ${T}_{\mathrm{eq}}\approx 2200$ K required a larger average condensate size. What the condensate size is in a hot Jupiter with a dayside near 3000 K has not been well modeled.

Based on the toy model diagrammed in Figure 10, with more spectrally resolved observations of hot Jupiters' daysides, it may be possible to observationally constrain the general condensate size in these planets' atmospheres by mapping out where stratospheric temperature inversions occur in extremely hot Jupiters. Though this assumes that cold traps are the dominant mechanism to inhibit inversion, it will also be necessary to determine whether one can even speak of a single "condensate size" for all hot Jupiters. To more adequately assess the role that cold traps play in hot Jupiters' atmospheres, we therefore need more 2D or 3D atmosphere models, such as those in Parmentier et al. (2016) and Lee et al. (2016), that include the variation of temperature structure as a function of longitude and latitude, models of particle and condensate growth processes within hot Jupiters' atmospheres, and spectrally resolved observations of more systems to validate the results of both.