Exploring the Atmospheric Dynamics of the Extreme Ultrahot Jupiter KELT-9b Using TESS Photometry

The KELT-9 system was observed by camera 1 of the TESS spacecraft during Sectors 14 (UT 2019 July 18 to August 14) and 15 (UT 2019 August 15 to September 11). KELT-9 is listed in the TESS input catalog (TIC; Stassun et al. 2018) as TIC 16740101 and is included in the list of preselected target stars that are observed with a 2 minute cadence using an 11 × 11 pixel subarray centered on the target. The photometric data were processed through the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016), hosted at the NASA Ames Research Center, which is largely based on the predecessor Kepler mission pipeline (Jenkins et al. 2010, 2017).

Even after instrumental systematics corrections by the SPOC pipeline, there are residual long-term trends in the photometry. We corrected for these by simultaneously fitting generalized polynomials in time to each data segment alongside the astrophysical phase-curve model. The order of each segment's polynomial detrending function was optimized by minimizing the Bayesian information criterion (BIC); across the 10 segments, the optimal orders ranged from zero to 2.

After downloading the data from the Mikulski Archive for Space Telescopes (MAST), we proceeded to analyze the light curves following methodologies that are largely identical to the ones described in previous TESS phase-curve papers (Shporer et al. 2019; Wong et al. 2020b). In the analysis presented in this paper, we used the presearch data conditioning (PDC; Smith et al. 2012; Stumpe et al. 2014) light curves from the SPOC pipeline, which have been corrected for instrumental systematics and contamination from nearby stars. We also experimented with fitting the uncorrected simple aperture photometry (SAP) light curves. The SAP data show more significant long-term systematic trends. Correcting these trends using linear combinations of the cotrending basis vectors (CBVs) provided by MAST removed a majority of the systematics and yielded phase-curve fit parameters that are consistent with those derived using the PDC light curves. However, some long-term temporal trends persist, and the corrected SAP light curves have larger scatter than the corrected PDC light curves, resulting in less precise fitted astrophysical parameters. Likewise, we constructed light curves from the target pixel files and defined various extraction apertures but found that the quality of the detrending (using either CBVs or polynomials in time) was generally significantly poorer than in the case of the SPOC-generated SAP light curve.

Data obtained during Sectors 14 and 15 were occasionally affected by severe stray-light contamination on the detector. The SPOC pipeline flagged those exposures and set the fluxes as NaN. In addition, data points preceding and following the episodes of stray-light contamination exhibited higher-than-normal background fluxes and were also flagged in the PDC data. Out of 19,337 2 minute frames obtained by TESS in Sector 14, 8603 (44%) were flagged, the vast majority due to stray light; in Sector 15, 37% of the data points (6973 out of 18,757) were flagged. We removed all flagged data points from the light curve and then applied a 16-point moving median filter to the photometric time series to trim 3σ outliers. The outlier trimming removed an additional 1.6% (172 out of 10,734) and 1.7% (195 out of 11,784) of the points from the Sectors 14 and 15 time series, respectively.

The transits of KELT-9b show significant deviations from the typical transit light-curve shape. These asymmetries are attributable to the gravity-darkened photosphere of the rapidly rotating host star and have been reported for a handful of other transiting systems with hot host stars, including Kepler-13A (Barnes et al. 2011), Kepler-462 (Ahlers et al. 2015), HAT-P-70 (Zhou et al. 2019), and MASCARA-4 (Ahlers et al. 2020b). A detailed, dedicated analysis of the gravity-darkened transit of KELT-9b in TESS photometry is described in Ahlers et al. (2020a). In our phase-curve analysis, we trimmed the transits from the light curve prior to fitting.

During normal spacecraft operation, the photometric observations during each Sector are interrupted by scheduled momentum dumps. These events can cause discontinuities in flux and, in some cases, short-term photometric variability before and/or after. These instrumental artifacts are not fully removed by the PDC pipeline and can present difficulties for the systematics detrending during our fits. Following previous work, we split each orbit's time series into segments separated by the momentum dumps. We found that the first segment of each spacecraft orbit exhibits a discernible initial flux ramp. After inspecting the binned residuals from the phase-curve fit to each segment, we removed the transient flux features by trimming the first 0.25 day worth of data from segments 1, 4, and 7 and the first 1.5 days worth of data from segment 9. The outlier-trimmed, median-normalized light curve is shown in the top panel of Figure 1, with the momentum dumps and light-curve segments indicated in blue and the trimmed points marked in red and green. The final light curve contains 17,489 data points.

Zoom In Zoom Out Reset image size

Our phase-curve fitting was done using the Python-based ExoTEP pipeline (e.g., Benneke et al. 2019; Shporer et al. 2019; Wong et al. 2020a). The secondary eclipse light curve λ_e(t) is modeled using batman (Kreidberg 2015) and defined to be zero at mid-eclipse and unity out of eclipse. The phase-curve variation is modeled as an nth-order Fourier series at the orbital phase (ϕ ≡ (t − T₀)/P, where T₀ is the transit time, and P is the orbital period) and is split into two components, one describing the photometric signal from the host star and another containing terms attributed to the planet's brightness variation:

$$\begin{eqnarray}&&{\psi }_{p}(t)=\bar{{f}_{p}}+{B}_{1}\cos (2\pi \phi +\delta ),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\psi }_{* }(t)=1+\displaystyle \sum _{k=1}^{n}{A}_{k}\sin (2\pi k\phi )+\displaystyle \sum _{k=2}^{n}{B}_{k}\cos (2\pi k\phi ).\end{eqnarray} \tag{ 2 }$$

The planet's atmospheric brightness modulation around its mean value $\bar{{f}_{p}}$ has a semiamplitude B₁ and a phase shift δ. It follows that the secondary eclipse depth (i.e., the planet's dayside flux at superior conjunction) is precisely the value of ψ_p(t) at phase π, while the nightside flux of the planet's atmosphere is calculated at phase 0: ${D}_{d}=\bar{{f}_{p}}+{B}_{1}\cos (\pi +\delta )$ and ${D}_{n}=\bar{{f}_{p}}+{B}_{1}\cos (\delta )$.

In general, the brightness modulation at visible wavelengths is a combination of contributions from reflected and emitted light. While the thermal emission component is well described by a simple sinusoidal function (see, for example, the numerous Spitzer phase-curve studies; e.g., Wong et al. 2016; Beatty et al. 2019), the reflected component can take on a more complicated form. In the case of Lambertian scattering, the brightness modulation contains an additional second-order term at the first harmonic of the orbital phase (e.g., Faigler & Mazeh 2015). However, as discussed in Section 4.1, the geometric albedo of KELT-9b is consistent with zero, indicating a negligible reflected component in the TESS phase curve. Therefore, the model in Equation (1) is expected to provided an accurate description of the atmospheric brightness modulation.

All remaining photometric variability is associated with the host star. Out of these terms, the first harmonic of the cosine B₂ is attributable to the leading term of the photometric variation arising from ellipsoidal distortion of the host star (e.g., Morris 1985; Morris & Naftilan 1993; Pfahl et al. 2008). The second-order term of the ellipsoidal distortion component is at the fundamental of the orbital phase; however, the predicted amplitude of this term is more than an order of magnitude smaller than the leading-order term and well below the typical error bars on the sinusoidal amplitudes from our fits.

We included the fundamental of the sine term in the star's flux model. A signal at this frequency and phase can arise from Doppler boosting (e.g., the TESS phase curve of WASP-18; Shporer et al. 2019), which stems from periodic blue- and redshifting of the stellar spectrum, as well as modulations in the photon emission rate in the observer's direction due to the radial velocity of the host star induced by the orbiting planet (e.g., Shakura & Postnov 1987; Loeb & Gaudi 2003; Shporer et al. 2010). In a free, unconstrained fit, the amplitude of this term (A₁) is degenerate with a phase shift in the planet's brightness modulation in the combined phase curve. To address this degeneracy while still including the term in the phase-curve modeling, we computed the predicted strength of the Doppler boosting signal for the KELT-9 system using the formulation in Esteves et al. (2013) and obtained 2.1 ± 0.3 ppm, which we used as a Gaussian prior on A₁ in the joint fits presented in this paper.

The combined normalized astrophysical phase curve and secondary eclipse model is as follows:

$$\begin{eqnarray}&&\psi (t)=\displaystyle \frac{{\psi }_{* }(t)+{\lambda }_{e}(t){\psi }_{p}(t)}{1+\bar{{f}_{p}}}.\end{eqnarray} \tag{ 3 }$$

We modeled any remaining long-term photometric variability in the PDC light curves due to residual uncorrected systematics or stellar variability using generalized polynomials in time:

$$\begin{eqnarray}&&{S}_{N}^{\{i\}}(t)=\displaystyle \sum _{j=0}^{N}{c}_{j}^{\{i\}}{\left(t-{t}_{0}\right)}^{j}.\end{eqnarray} \tag{ 4 }$$

Here t₀ is the time of the first exposure in segment i, and N is the order of the detrending polynomial. We optimized the polynomial order of the detrending polynomial for each segment by carrying out individual segment phase-curve fits and considered both the Akaike information criterion ($\mathrm{AIC}\equiv 2\gamma -2\mathrm{log}L$) and the BIC ($\mathrm{BIC}\equiv \gamma \mathrm{log}m-2\mathrm{log}L$); here γ is the number of free parameters in the fit, m is the number of data points in the segment, and L is the maximum log-likelihood. For the 10 data segments, minimizing the AIC and BIC yielded the same optimal polynomial orders: 2, 1, 1, 0, 0, 1, 4, 4, 1, and 3.

The combined phase-curve and systematics model of KELT-9 is

$$\begin{eqnarray}&&f(t)={S}_{N}^{\{i\}}(t){| }_{i=1-6}\times \psi (t).\end{eqnarray} \tag{ 5 }$$

The systematics-removed light-curve segments are shown in the bottom panel of Figure 1.

We detected a separate periodic variation in the TESS photometry with a frequency unrelated to the orbital period of KELT-9b. The top left panel of Figure 2 shows the Lomb–Scargle periodogram of the residuals from a joint light-curve fit using the model in Equation (5). A strong signal with a characteristic period of roughly 0.316 day is evident. This additional periodicity was also detected in the Spitzer 4.5 μm light curve (Mansfield et al. 2020).

Zoom In Zoom Out Reset image size

We obtained light curves produced by the Quick Look Pipeline (Huang et al. 2018) for all nearby sources brighter than 13th magnitude in the TESS band that lie within 6' of KELT-9. When inspecting the individual periodograms, we did not find any significant periodicity at the frequency seen in the KELT-9 light curve. We therefore attribute this periodic signal to stellar pulsations on KELT-9. In our light-curve analysis, we addressed this contribution to the overall photometric variability by simultaneously fitting for the stellar pulsation signal in our joint phase-curve analysis. The stellar pulsations were modeled by a simple sinusoid function with three additional free parameters,

$$\begin{eqnarray}&&{\rm{\Theta }}(t)=1+\alpha \sin \left(2\pi \xi \right)+\beta \cos \left(2\pi \xi \right),\end{eqnarray} \tag{ 6 }$$

where ξ ≡ (t − T₀)/Π is the phase given a stellar pulsation period Π, and α and β are the semiamplitudes of the Fourier components of the pulsation signal. This stellar pulsation model is multiplied to the host star's photometric modulation described in Equation (2) to produce the corrected astrophysical phase curve and secondary eclipse model:

$$\begin{eqnarray}&&\psi ^{\prime} (t)=\displaystyle \frac{{\psi }_{* }(t){\rm{\Theta }}(t)+{\lambda }_{e}(t){\psi }_{p}(t)}{1+\bar{{f}_{p}}}.\end{eqnarray} \tag{ 7 }$$

The full phase-curve model used in our light-curve fits, incorporating both stellar pulsations and systematics, is as follows:

$$\begin{eqnarray}&&f^{\prime} (t)={S}_{N}^{\{i\}}(t){| }_{i=1-6}\times \psi ^{\prime} (t).\end{eqnarray} \tag{ 8 }$$

We have assumed that the stellar pulsations can be modeled as a simple sinusoid. In general, the shape of the photometric modulation induced by pulsations need not follow such a simple functional form and may consist of additional harmonics and/or other characteristic frequencies. However, we found that including just the fundamental at the stellar pulsation period removed all significant signals in the periodogram of the resultant residuals, as illustrated in the bottom left panel of Figure 2. As an independent check, we experimented with adding higher-order terms at the stellar pulsation period or sinusoidal terms with different periods and found no improvement to the AIC or BIC.

As discussed previously, we removed the transits from the TESS light curve prior to fitting. To constrain the shape of the secondary eclipse, we placed Gaussian priors on the impact parameter b, scaled semimajor axis a/R_*, and planet–star radius ratio R_p/R_* using the results from a dedicated analysis of the gravity-darkened TESS transits (Ahlers et al. 2020a): b = 0.14 ± 0.04, a/R_* = 3.36 ± 0.14, and R_p/R_* = 0.0791 ± 0.0018. We also used the transit ephemeris derived from the same analysis as the priors: T₀ = 2,458,683.4449 ± 0.0003 BJD_TDB and P = 1.4811235 ± 0.0000011 days.

We carried out an ensemble of joint fits of all 10 data segments from Sector 14 and 15 TESS photometry, varying the number of phase-curve terms included in the astrophysical model. Some of these fits are listed and numbered in Table 1, with the variables included as free parameters in each instance indicated by check marks. Phase-curve harmonics up to k = 3 were considered in the astrophysical phase-curve model. In order to examine whether accounting for the stellar pulsations has a notable effect on the phase-curve results, we carried out a joint fit using the phase-curve model in Equation (5), which does not account for stellar pulsations (#1; fit A), as well as one that simultaneously fit for the phase curve and stellar pulsations (#2; fit B) following Equations (6)–(8). We also checked for eccentricity in the system via a joint fit that allowed the orbital eccentricity e and argument of periastron ω to vary freely (#3). The remaining joint fits (#4–#7) contained subsets of the phase-curve terms in fit B, thereby probing the statistical significance of the various detected components independent of the formal uncertainties on the associated amplitudes derived from the fits.

The ExoTEP pipeline utilizes the affine-invariant Markov Chain Monte Carlo (MCMC) ensemble sampler emcee (Foreman-Mackey et al. 2013) to simultaneously calculate the posterior distributions of all free parameters. In addition to the astrophysical and systematics detrending parameters, we defined a uniform per-point uncertainty σ_i for each of the 10 segments to ensure that the reduced χ² value is unity and self-consistently generate realistic uncertainties on the fitted parameters. The number of walkers was set to four times the number of free parameters. The length of each chain was 50,000 steps, and we used only the last 40% of each chain for calculating the posterior distributions. To check for convergence, we applied the Gelman–Rubin test (Gelman & Rubin 1992) and ensured that the diagnostic value $\hat{R}$ was well below 1.1.

At superior conjunction, the planet–star flux ratio (i.e., secondary eclipse depth) can be expressed as a function of various fundamental properties of the system (e.g., Charbonneau et al. 2005):

$$\begin{eqnarray}&&{D}_{d}={\left(\displaystyle \frac{{R}_{p}}{{R}_{* }}\right)}^{2}\displaystyle \frac{\int \tau (\lambda ){F}_{p,\nu }(\lambda ,{T}_{B,p,\mathrm{day}})d\lambda }{\int \tau (\lambda ){F}_{* ,\nu }(\lambda ,{T}_{\mathrm{eff}})d\lambda }+{A}_{g}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{2}.\end{eqnarray} \tag{ 9 }$$

The first term describes the contribution from the planet's thermal emission. Here R_p/R_* is the planet–star radius ratio, T_eff is the host star's effective temperature, and T_B,p,day is the disk-averaged brightness temperature of the planet's observer-facing hemisphere. Here F_p,ν(λ, T_p) and F_*,ν(λ, T_eff) are the photon-weighted planetary and stellar fluxes, which are computed from the corresponding unweighted functions (in wavelength units) via the multiplicative conversion factor λ/hc. These fluxes are integrated across the TESS transmission function τ(λ). The second term represents reflected light off the planet's dayside atmosphere, where A_g is the geometric albedo in the TESS bandpass, and a is the orbital semimajor axis.

To account for the uncertainties in the stellar parameters and propagate them to the brightness temperature estimates, we calculated the integrated stellar flux in the TESS bandpass for a grid of PHOENIX models (Husser et al. 2013) spanning the ranges T_eff = [8000, 12,000] K, [M/H] = [−1.0, +0.5], and $\mathrm{log}g=[3.5,4.5]$. To obtain a simple empirical model of the integrated stellar flux, we fit a linear polynomial in (T_eff, [M/H], $\mathrm{log}g$) to the calculated values; including higher-order terms did not improve the fit. We computed the posterior distribution of the dayside brightness temperature using a Monte Carlo sampling method: in order to account for uncertainties in the system parameters, we allow them to vary in the MCMC fit while constraining them with Gaussian priors derived from their published values and uncertainties. For the stellar parameters, we use the measurements reported by Gaudi et al. (2017; T_eff = 10,170 ± 450 K, $\mathrm{log}g=4.093\pm 0.014$, [M/H] = −0.03 ± 0.020), while for R_p/R_* and a/R_*, we use the values from Gaudi et al. (2017; R_p/R_* = 0.08004 ± 0.00041) and Mansfield et al. (2020; a/R_* = 3.153 ± 0.011), which have significantly smaller uncertainties than the TESS-band measurements published in Ahlers et al. (2020a). We assume a blackbody for the planet's emission spectrum. Both direct measurements and atmospheric models demonstrate that hot Jupiters tend to have very low reflectivities (e.g., Heng & Demory 2013). Furthermore, the extremely hot dayside of KELT-9b is expected to preclude the formation of any condensate clouds that might enhance the geometric albedo (see Section 4.4). Here we allow for a range of geometric albedos spanning zero through 0.2.

For zero geometric albedo, we derive a dayside brightness temperature of T_B,p,day = 4600 ± 100 K. Across the range of geometric albedos we consider, the dayside temperature varies from 4360 to 4600 K. The measured dayside temperature of KELT-9b is comparable to those of mid-K dwarfs. An analogous calculation using the measured nightside flux (and no reflected light) yields a very high nightside temperature of T_B,p,night = 3040 ± 100 K. This means that the nightside of KELT-9b is hotter than the dayside of almost all known exoplanets and comparable to the atmosphere of a mid-M dwarf.

The difference between the dayside and nightside temperatures of a planet reflects the efficiency of heat transport from the dayside to the nightside, e.g., by zonal winds. The more inefficient the heat recirculation is, the larger the expected day–night temperature contrast. Following the simple thermal balance considerations in Cowan & Agol (2011), the relationship between irradiation, heat recirculation, and dayside and nightside temperatures can be expressed as

$$\begin{eqnarray}&&{T}_{B,p,\mathrm{day}}={T}_{* }\sqrt{\displaystyle \frac{{R}_{* }}{a}}{\left(1-{A}_{B}\right)}^{1/4}{\left(\displaystyle \frac{2}{3}-\displaystyle \frac{5}{12}\epsilon \right)}^{1/4},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{T}_{B,p,\mathrm{night}}={T}_{* }\sqrt{\displaystyle \frac{{R}_{* }}{a}}{\left(1-{A}_{B}\right)}^{1/4}{\left(\displaystyle \frac{\epsilon }{4}\right)}^{1/4},\end{eqnarray} \tag{ 11 }$$

where A_B is the Bond albedo (i.e., the fraction of total incident stellar irradiation that is reflected), and is a fiducial parameter that indicates the relative efficiency of heat transport, with zero corresponding to no day–night heat recirculation and 1 corresponding to full recirculation.

By simultaneously fitting to the measured dayside and nightside temperatures in the TESS bandpass using Equations (10) and (11) and constraining T_* and a/R_* with the same Gaussian priors presented above while keeping and A_B unconstrained, we obtain ${A}_{B}={0.19}_{-0.11}^{+0.12}$ and  = 0.39 ± 0.05. These results indicate relatively efficient day–night heat transport. All of the quantities derived in this subsection are listed in Table 2. The measured heat recirculation parameter will be placed into context with other planets in Section 4.4.

4.1.1. Comparison to Mansfield et al. (2020)

In binary systems, the orbiting companion raises a tidal bulge on the surface of the host star, which incurs periodic variations in the star's sky-projected area. The observer-facing profile of the host star is expected to come to maximum at quadrature, producing a photometric modulation with a leading term at the first harmonic of the cosine. The amplitude of this variation can be related to fundamental properties of the system via the theoretical formalism described in Kopal (1959). Using the notation of Morris (1985) and Morris & Naftilan (1993), we can express the dominant amplitude B₂ as

$$\begin{eqnarray}&&{B}_{2}=-{Z}_{\mathrm{ellip}}\displaystyle \frac{{M}_{p}}{{M}_{* }}{\left(\displaystyle \frac{{R}_{* }}{a}\right)}^{3}{\sin }^{2}i,\end{eqnarray} \tag{ 12 }$$

where M_p/M_* is the planet–star mass ratio, and the prefactor Z_ellip depends on the star's linear limb- and gravity-darkening coefficients and is of order unity. The planet also experiences ellipsoidal distortion, but the amplitude is smaller than the stellar signal by a factor of at least ${({R}_{p}/{R}_{* })}^{3}{\bar{f}}_{p}/{({M}_{p}/{M}_{* })}^{2}\sim 0.5 \%$ and can thus be safely ignored.

Using limb- and gravity-darkening coefficient values computed by Claret (2017) in the TESS bandpass for the nearest available combination of stellar parameters ($\mathrm{log}{T}_{\mathrm{eff}}=3.978$, $\mathrm{log}g=4.00$, [Fe/H] = 0.10, u₁ = 0.3494, γ₁ = 0.4355), we calculate a predicted ellipsoidal distortion amplitude of 44 ± 6 ppm. Beyond the leading-order term, there are smaller contributions at the second harmonic of the cosine and the fundamental of the orbital period (i.e., $\cos (6\pi \phi )$ and $\cos (2\pi \phi )$). Using the formulas in Wong et al. (2020c), we compute these expected amplitudes to be 2.7 ± 0.4 and 1.6 ± 0.2 ppm, respectively.

Looking back to the results from our phase-curve fit, we find that while the total amplitude of the photometric variation at the first harmonic (39.6 ± 4.5 ppm) is consistent with the predicted amplitude of 44 ± 6 ppm, there is a very significant phase offset in the modulation relative to the expectation. In Figure 4, we plot the measured first harmonic modulation (solid curve) alongside the predicted signal at the expected phasing (dotted curve).

We also measured a phase-curve variation at the second harmonic with an amplitude of ${15.5}_{-4.3}^{+4.4}$ ppm, and once again, there are significant discrepancies in both amplitude and phasing between the measured signal and the predicted second-order term of the ellipsoidal distortion. Given the relatively low significance of the second harmonic detection, the measured signal may be due to residual short-term instrumental systematics or other noise in the photometry. Looking at Figure 3, we can see that the phase-folded residuals from the joint fit contain some traces of time-correlated noise. When the KELT-9 system is revisited by TESS during the Extended Mission, we will obtain significantly more photometry, from which more robust constraints can be placed on the signal at the second harmonic. We note that unexpected second harmonic phase-curve signals have been detected in other systems, including HAT-P-7 (Armstrong et al. 2016) and Kepler-13A/KOI-13 (Esteves et al. 2013; Shporer et al. 2014).

Mismatches between the measured and expected ellipsoidal distortion signals have been previously reported in other systems with early-type host stars, including KOI-54 (Burkart et al. 2012), KOI-74 (Rowe et al. 2010; van Kerkwijk et al. 2010; Ehrenreich et al. 2011; Bloemen et al. 2012), and KOI-964 (Carter et al. 2011; Wong et al. 2020c). These discrepancies might arise due to an oversimplification of the tidal dynamics of hot stars with largely radiative envelopes, where interactions between the dynamical tide and the gravitational influence of the orbiting companion can yield large steady-state deviations in both amplitude and relative phase from the equilibrium tide predictions (Pfahl et al. 2008). Modeling of the tidal response of KOI-54 (M = 2.32 ± 0.10 M_⊙, R = 2.19 ± 0.03 R_⊙)—an A-type star similar to KELT-9 (M = 2.32 ± 0.16 M_⊙, R = 2.418 ± 0.058 R_⊙; Borsa et al. 2019)—revealed phase shifts in the resultant ellipsoidal distortion modulation of up to a quarter of an orbit, depending on the proximity of the orbital period to resonances and the damping of the stellar oscillation modes (Burkart et al. 2012).

The host star's rapid rotation, combined with the peculiar near-polar orbit of KELT-9b, may also explain the phase shift in the measured first harmonic modulation. Ahlers et al. (2020a) modeled the gravity-darkened transits of KELT-9b from the TESS light curve and simultaneously resolved the three-dimensional stellar spin axis and the orbital misalignment. They measured a sky-projected stellar obliquity of λ = −88° ± 15°, a stellar inclination of ψ = −52°(+8°, −7°), and a spin–orbit misalignment of φ = 87° (+10°, −11°). When combined with the orbital inclination measurement of i = 872 ± 04, this means that KELT-9b passes almost perfectly over the stellar poles and that the spin axis of the host star lies in a plane that is nearly perpendicular to the sky plane. Furthermore, given the stellar inclination, KELT-9b passes over the stellar pole roughly one-sixth of an orbit prior to mid-transit and mid-eclipse.

KELT-9 has a rapid spin, with a measured sky-projected rotational velocity of $v\sin {i}_{* }=114\pm 1.3$ km s⁻¹ (Gaudi et al. 2017) and a stellar rotation period of ${16}_{-4}^{+3}$ hr; the rapid rotation also generates a substantial stellar oblateness of 0.089 ± 0.017 (Ahlers et al. 2020a). Due to gravity darkening, the star is brightest at the poles, and as the planet orbits around the host star, the nonhomogeneous stellar effective temperature distribution alters the irradiation of the planet. Meanwhile, the stellar oblateness changes the size of the projected stellar disk that the planet's dayside sees across the orbit. These two effects combine constructively to produce an irradiance modulation that comes to maximum twice per orbit, when the planet passes near the stellar pole. The corresponding periodic photometric modulation stemming from these processes contributes an additional signal in the planet's flux at the first harmonic of the orbital period.

From the gravity darkening and oblateness modeling in Ahlers et al. (2020a), we produced an effective temperature map of the host star and generated the total irradiance I experienced by KELT-9b as a function of orbital phase, using the model described in Ahlers (2016). The equilibrium dayside temperature of KELT-9b is related to the incident irradiance via T_day ∝ I^1/4, and we computed the relative change in thermal emission from the planet associated with the varying dayside temperature (see the numerator in the first term of Equation (9)). After scaling the resultant modulation by the secondary eclipse depth (i.e., dayside flux relative to the star's brightness), we obtained the predicted phase-curve contribution from the planet's time-varying irradiance, plotted in Figure 5.

Zoom In Zoom Out Reset image size

From the figure, it is evident that the phasing of this irradiance signal is significantly offset from the expected ellipsoidal distortion modulation. Taking into account the uncertainties in the stellar spin axis solution, the pre-eclipse irradiance maximum occurs at an orbital phase of 0.36 ± 0.03. Notably, the alignment of the irradiance signal is consistent at 1.5σ with that of the measured first harmonic phase-curve component from the TESS light curve, which has a pre-eclipse maximum at an orbital phase of ${0.408}_{-0.012}^{+0.011}$.

For additional observational perspectives on the role of periodic atmospheric heating modulations, we turn to the published full-orbit 4.5 μm thermal phase curve. Mansfield et al. (2020) detected a significant phase-curve signal at the first harmonic of the orbital period with a semiamplitude of 89 ± 22 ppm and a pre-eclipse maximum that occurs at an orbital phase of 0.455 ± 0.022. In Figure 5, we plot this measured first harmonic modulation, with the amplitude scaled by the ratio of KELT-9b's measured relative dayside brightnesses in the TESS and Spitzer 4.5 μm bandpasses: 650/3131. The adjusted amplitude is statistically identical to that of the irradiation model, while the phase alignment of the maxima is consistent with the measurement from the TESS phase curve at the 1.9σ level.

Lastly, if we assume that the tidal response of the star occurs as predicted from theory, we can remove the contribution of ellipsoidal distortion to the overall first harmonic signal by subtracting the expected ellipsoidal modulation. The resultant curve is plotted in red in Figure 5. Notably, the phasing of this ellipsoidal distortion-corrected signal is almost identical to the measured first harmonic term from Spitzer. Since the expected ellipsoidal distortion amplitude at thermal wavelengths is negligible compared to the reported signal from Mansfield et al. (2020), any observed modulation at 4.5 μm is most likely due to thermal emission from the planet.

The similarity between the ellipsoidal distortion-corrected TESS and Spitzer phase curves at the first harmonic lends strong support to the hypothesis that the unusual phase shift we observe is caused by the response of the planet's atmosphere to the time-varying irradiation across its orbit. Of course, the assumption that the star behaves as predicted from equilibrium tide theory is subject to many caveats, as explained earlier, and deviations in the ellipsoidal distortion from expectations will, in turn, affect the corresponding contribution from time-varying irradiance. Future multiwavelength phase-curve observations of KELT-9 and detailed atmospheric modeling of KELT-9b's response to time-varying irradiation can help break the degeneracy between the individual phase and amplitude discrepancies stemming from ellipsoidal distortion and irradiation.

The growing number of well-measured phase curves has motivated the search for systematic trends that may elucidate the underlying physical processes driving the atmospheric dynamics in exoplanet atmospheres. KELT-9b lies at the high-temperature extreme of the known exoplanets, and we now place the planet into the wider context of hot gas giants.

The most salient trend that has emerged from phase-curve studies is the relationship between day–night temperature contrast and the level of stellar irradiation. As discussed in the Introduction, phase-curve observations have revealed a systematic increase in the day–night temperature contrast with increasing dayside equilibrium temperature up to roughly 2500 K, after which the trend appears to reverse (e.g., Schwartz et al. 2017; Zhang et al. 2018; Keating et al. 2019). This behavior is largely consistent with the predictions of general circulation models (e.g., Perna et al. 2012; Perez-Becker & Showman 2013; Komacek & Showman 2016). The addition of H₂ dissociation and recombination to the energy balance of the hottest planets has been invoked to explain the turnaround in the observed trend at ∼2500 K (Bell & Cowan 2018; Komacek & Tan 2018; Tan & Komacek 2019), though only a handful of planets in this temperature range have previously been observed across the full orbital phase.

Our analysis of the KELT-9b phase curve allows us to test the predictions of the H₂ recombination hypothesis. Keating et al. (2019) calculated heat recirculation parameter values for 13 hot Jupiters and brown dwarfs from Spitzer phase curves, as defined in Section 4.1. In Figure 6, we plot the published values from their work along with our measurements for KELT-9b. Given the different nightside brightness temperatures and recirculation efficiencies, as well as the possibility that the TESS observations probe systematically different pressure levels than Spitzer measurements, we have included the values for both bandpasses. The equilibrium temperature ${T}_{\mathrm{eq}}={T}_{* }\sqrt{{R}_{* }/2a}$ is used as a proxy for the level of stellar irradiation on the dayside atmosphere.

Zoom In Zoom Out Reset image size

In both bandpasses, the value of KELT-9b is elevated relative to the average value for cooler planets with dayside equilibrium temperatures between 2000 and 3000 K. The inclusion of KELT-9b reinforces the previously suggested turnaround in the trend of day–night temperature contrast with irradiation level predicted by the H₂ recombination hypothesis. Further phase-curve studies of exoplanets are needed, particularly for planets with equilibrium temperatures between 3000 and 3500 K, in order to robustly confirm or refute these apparent trends. No such planets have hitherto been discovered.

Keating et al. (2019) previously reported a nearly constant nightside temperature of about 1100 K for a sample of 12 hot Jupiters, which they attributed to the presence of opaque nightside clouds. In Figure 7, we have added measured nightside temperatures from Mansfield et al. (2020) and our study to the values from Keating et al. (2019). The addition of the KELT-9b data reveals a tentative increasing trend in nightside temperature as a function of insolation across the full data set. Crucially, the significantly higher nightside temperature of KELT-9b supports the prediction that recombination of H₂ contributes significantly to nightside heating for the hottest exoplanets (Bell & Cowan 2018; Komacek & Tan 2018).

Zoom In Zoom Out Reset image size

In Figure 7, we have added estimates of the nightside temperature of KELT-9b given various mixing ratios of atomic and molecular hydrogen in the planet's atmosphere, following the methodology of Kitzmann et al. (2018). Assuming a pure hydrogen gas in chemical equilibrium (Gail & Sedlmayr 2014; Heng et al. 2016) and an approximate photospheric pressure of 0.1 bar consistent with Bell & Cowan (2018) and Kitzmann et al. (2018), the nightside temperature of 3050 K corresponds to an atmosphere with 50% H₂ and 50% H. If we consider the full statistical range of nightside temperatures obtained (3050 ± 110 K), then the percentage of H₂ present is 40%–60%. These estimates place KELT-9b in a very different regime than that considered in Bell & Cowan (2018), where the nightside is assumed to be H₂-dominated.

In Figure 7, we have also overlaid condensation temperatures at 0.1 bar and solar metallicity for enstatite (MgSiO₃), titanium oxide (TiO), and corundum (Al₂O₃). Enstatite and corundum are, respectively, the least and most refractory of the mineral clouds that are expected to form in thermochemical equilibrium (e.g., Burrows et al. 2006), while titanium oxide is believed to provide seed particles for high-temperature clouds (e.g., Helling & Woitke 2006). The nightside temperatures of KELT-9b in the TESS and Spitzer 4.5 μm bandpasses are ∼500–1000 K higher than even the condensation temperature of corundum, implying that the nightside atmosphere is cloud-free at photospheric pressures.

To study the planetary atmosphere in more detail and demonstrate what future observations may reveal, we used the open-source HELIOS radiative transfer code (Malik et al. 2017, 2019) and the HELIOS-O ray-tracing code (Gaidos et al. 2017) to model the dayside, nightside, and transmission spectra. In this model framework, the emission spectra are computed self-consistently in radiative–convective equilibrium such that the thermal structure is an outcome of the computation, rather than an assumption. We fixed the stellar effective temperature and the planet-to-star radius ratio to the previously published values from Gaudi et al. (2017) and Mansfield et al. (2020) that we used when fitting for the blackbody brightness temperatures in Section 4.1: T_eff = 10,170 K and R_p/R_* = 0.08004. For spectral line lists, we used the ExoMol spectroscopic databases for H₂O (Barber et al. 2006; Polyansky et al. 2018), CH₄ (Yurchenko & Tennyson 2014), VO (McKemmish et al. 2016), and TiO (McKemmish et al. 2019); for CO, we used Li et al. (2015); for CO₂, we used the HITEMP database (Rothman et al. 2010); and spectral lines for Fe, Fe⁺, Ti, Ti⁺, Ca, Ca⁺, Na, K, Y, Y⁺, AlO, and SH were taken from Kurucz & Bell (1995). In addition, we included the continuum opacities from H⁻ (John 1988) and collision-induced absorption (H₂–H₂, H₂–He, H–He; Abel et al. 2011; Karman et al. 2019), as calculated using the open-source HELIOS-K opacity calculator (Grimm & Heng 2015).¹⁹ Equilibrium chemistry was computed using the open-source FastChem code, which includes almost 600 chemical species (Stock et al. 2018). Solar metallicity was assumed. For the dayside emission spectra, we computed an interpolated PHOENIX stellar model and used it to irradiate the model atmosphere while varying the heat recirculation parameter (as defined in Section 4.1). For the nonirradiated nightside, the only free parameter is the interior temperature. In the case of the transmission spectra, we assumed isothermal transit chords.

Figure 8 shows a sampling of models for the dayside emission, nightside emission, and transmission spectra. To facilitate comparison with the measured relative planetary fluxes, the planetary emission spectra are shown as a ratio with the stellar spectrum: F_p/F_*. We have overplotted the measured secondary eclipse depths and nightside fluxes from this paper, as well as the values from Gaudi et al. (2017) and Mansfield et al. (2020). For the dayside, we find that the model predicts planetary emission spectra that are close to a blackbody; we note that the spectral features apparent in the modeled F_p/F_* curves stem from absorption features in the star's spectrum, not the planet's. The three data points are well matched by the model with  = 0.4, in excellent agreement with the value of  = 0.39 ± 0.05 we deduced from thermal balance calculations based on the measured dayside and nightside blackbody brightness temperatures (Section 4.1).

Zoom In Zoom Out Reset image size

The model nightside emission spectra are shown in the middle panel of Figure 8. Our measured TESS-band flux is in good agreement with the 3000 K model, while the Spitzer 4.5 μm flux lies intermediate between the 2500 and 3000 K models. All of the model spectra display large deviations from a blackbody spectrum primarily due to absorption features of H₂O and CO, as well as the wavelength-dependent H⁻ opacity. In particular, the broad CO absorption in the Spitzer 4.5 μm bandpass yields a predicted flux that is significantly lower than the corresponding flux when assuming a blackbody. Therefore, this modeling helps explain the discrepant nightside blackbody brightness temperatures measured in the TESS and Spitzer 4.5 μm bandpasses (see Section 4.1).

When examining the model–data comparisons, we recall that the forward models shown here assume fixed values for T_eff and R_p/R_*. However, both of these input parameters have relatively large measurement uncertainties, particularly in the case of stellar temperature (450 K). Changes to one or both of these parameters can incur wavelength-dependent variations in F_p/F_*. When experimenting with altering the values of these inputs within their respective 1σ confidence regions, we found that the relative model brightness at Spitzer wavelengths compared to the TESS-band brightness can vary by 10% or more. In other words, the band-integrated model points shown in the top two panels of Figure 8 have intrinsic error bars that are significant in relation to the distance between adjacent models.

The bottom panel of Figure 8 shows predictions for the transmission spectra of KELT-9b across the wavelength range 0.5–5 μm. We have chosen to normalize the model spectra to match the Spitzer 4.5 μm transit depth from Mansfield et al. (2020). We have also included the TESS-band transit depth from Ahlers et al. (2020a)—6260 ± 290 ppm—which has large error bars but is consistent with the full range of models shown. Meanwhile, we do not include the much larger transit depth from Gaudi et al. (2017), (R_p/R_*)² = 6770 ± 70 ppm; that measurement was derived from transit fits that did not account for the rotational bulge and gravity-darkened limb profile of the rapidly rotating host star (as was done in Ahlers et al. 2020a), despite being in a wavelength range where those effects are significant.

The most salient trend in the model spectra is the disappearance of the familiar water feature at 1.4 μm with increasing temperature. As the temperature increases from 2000 to 4000 K, the water abundance decreases by 3 orders of magnitude, while the H⁻ abundance increases by 2 orders of magnitude (see Figure 2 of Kitzmann et al. 2018). Given the measured temperatures of KELT-9b, we expect the transmission spectrum at low spectral resolution to be a largely featureless spectral slope across the optical and near-infrared through ∼2 μm, with the exception of the prominent alkali absorption peaks. The gradient of the spectral slope in this wavelength range is determined by the wavelength-dependent behavior of the H⁻ cross section (John 1988). Meanwhile, at longer wavelengths, our models predict large absorption features, consistent with previous modeling of UHJs (e.g., Lothringer & Barman 2019).

The uncertainties in all published broadband transit depths are significantly larger than the range spanned by the modeled transmission spectra. Future spectroscopic observations with the Hubble Space Telescope and/or James Webb Space Telescope, coupled with careful transit shape modeling accounting for the rotationally distorted and gravity-darkened star, may achieve the ∼10 ppm precision needed to adequately characterize the transmission spectrum, specifically the spectral gradient in the visible and near-infrared, and allow for a definitive observational test of atmospheric models. Likewise, spectroscopic secondary eclipse observations in the near-infrared and the measurement of absorption features in the dayside emission spectrum will provide strong constraints on the atmospheric chemistry.