© ESO 2021

1 Introduction

A brown dwarf that transits its host star once per orbital period offers a unique opportunity to characterize its atmosphere. The signal of the object can be separated from the light of the host star by time-differential techniques (Charbonneau & Deming 2007), that is, by carrying out a comparison of the brightness of the star-companion system based on either: out-of transit with in-transit or out-of secondary eclipse with in-eclipse. The wavelength-resolved application of these techniques is known as transmission spectroscopy for the transit event and emission spectroscopy for the secondary eclipse observation (see reviews by Kreidberg 2018; Deming et al. 2019; Madhusudhan 2019). In recent years, these methods have been applied very successfully for close-in planets (e.g., Sing et al. 2016; Evans et al. 2016; Tsiaras et al. 2019) because of the more frequent occultations and also because hotter planets have larger atmospheric scale heights and greater thermal emission.

Close-in companions are generally expected to be tidally locked as the tidal synchronization timescale is typically much shorter than the estimated stellar age (Guillot et al. 1996). Thus, they possess permanent day and night sides. The two successful techniques – transmission and emission spectroscopy – probe different regions: the day-night terminator and the day side, respectively, which are separated in longitude and exhibit different physical conditions. One way to observe the dependence of atmosphericparameters (such as temperature) on the planetary longitude is the observation of a full-orbit phase curve (see recent reviews by Shporer 2017; Parmentier & Crossfield 2018). At infrared wavelengths, the temperature difference between the day side and night side causes a brightness variation modulated with the planetary rotation due to varying thermal emission. At optical wavelengths, the brightness modulation is a combination of thermally emitted light and reflected stellar light, both having a maximum around secondary eclipse and a minimum around the transit event. However, typical amplitudes of phase curve modulations are small, amounting to 100 ppm of the stellar light at optical wavelengths (Esteves et al. 2015; Shporer 2017; Shporer et al. 2019) and 1000 ppm in the near-infrared (Beatty et al. 2014; Zhang et al. 2018; Kreidberg et al. 2018), with the variations being very gradual over timescales of days. Thus, exoplanet phase curve observationsare currently out of reach with regard to ground-based telescopes and can only be performed with space telescopes.

The Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) has measured the full-orbit phase curves of multiple close-in gas giants thus far, revealing the day- and night-side brightness temperatures at optical wavelengths as well as the phase curve offset between the brightest measured longitude and the most strongly irradiated longitude at the substellar point (e.g., Shporer et al. 2019; Wong et al. 2020a,b,c; Bourrier et al. 2020; Daylan et al. 2019; Jansen & Kipping 2020; Nielsen et al. 2020; von Essen et al. 2020a). These values provide information on the Bond albedo and the energy transport from the day to the night side (Cowan & Agol 2011). Additionally, the TESS phase curves allowed for constraints on the geometric albedo of individual objects (Shporer et al. 2019; Wong et al. 2020b,c; von Essen et al. 2020a) and the detection of gravitational interactions between planet and host star (Shporer et al. 2019; Wong et al. 2020a,c; Nielsen et al. 2020; von Essen et al. 2020a), allowing for independent mass constrains. When compared with near-infrared phase curve measurements of the same objects, it is possible to study the wavelength dependence of the phase curve parameters and distinguish between several models for the pressure-temperature profile of the planet atmosphere (Parmentier & Crossfield 2018).

In this work, we present our analysis of the TESS phase curve of KELT-1b (Siverd et al. 2012). KELT-1b was discovered with the Kilodegree Extremely Little Telescope (KELT, Pepper et al. 2012) and was found to have a mass in the brown-dwarf regime, via subsequent radial-velocity measurements. The first measurements of the thermal emission of the day side were provided by Beatty et al. (2014) with secondary eclipse observations with the Spitzer space telescope at 3.6 and 4.5 μm and with ground-based observations in the Sloan z’ band. A spectrally-resolved secondary eclipse was presented by Beatty et al. (2017) at the near-infrared H band, which showed a wavelength dependence of the brightness temperature of KELT-1b when compared to the previousSpitzer and z’ band measurements. The day-side emission spectrum favored a model with monotonically decreasing temperaturewith altitude, similar to an isolated brown dwarf of the same temperature. Ultra-hot Jupiters of lower mass showed indications of isothermal pressure-temperature profiles (Parmentier & Crossfield 2018), implying that the higher surface gravity of KELT-1b might play a major role in governing its atmospheric structure. Spitzer full-orbit phase curves at 3.6 and 4.5 μm were published by Beatty et al. (2019), with values for heat redistribution efficiency similar to those of Jupiter-mass ultra-hot planets.

This work is structured as follows. In Sect. 2, we present a detailed description of the observations. In Sect. 3, we present our analysis of the third light contamination, limb darkening coefficients, and primary transit parameters. We conclude that section by providing updated parameters for the host star and companion, which we use in our modeling. In Sect. 4, we detail the six different model components we take into consideration in order to characterize the atmospheric parameters of KELT-1b. In Sect. 5, we present the physical parameters derived from TESS observations and our analysis. In Sect. 6 we set our observations into context. We present our concluding remarks in Sect. 7.

2 Observations and data preparation

KELT-1 (TIC identifier 432549364) was observed by TESS between 8 and 31 of October 2019, during cycle 2 and sector 17, using camera 2. The photometric time series has a cadence of 120 seconds and it was constructed using the Science Processing Operations Center (SPOC) pipeline, based on the NASA Kepler mission pipeline (Jenkins et al. 2016; Jenkins 2017). The light curve analyzed in this work is the one constructed by the Presearch Data Conditioning (PDC) msMAP correction method. The timestamps are provided by SPOC in Barycentric Julian dates (TDB) and are, consequently, left unchanged.

It is well known that the PDC data might show systematic features in the photometry that have not, in turn, been observed in the Aperture Photometry (SAP) data (Bourrier et al. 2020; Shporer et al. 2019; Wong et al. 2020b). The unbinned KELT-1 photometry show variability during secondary eclipse and along the orbital phase. Consequently, to make a thoughtful choice between carrying out an unphysical detrending of the photometry (SAP), for instance, by detrending the data with a spline function or using the more physically-motivated detrending (PDC) for the cost of extra noise, we quantified the amount of correlated noise in the light curves computing and comparing their respective β values (Carter & Winn 2009). As was previously done by other authors (see e.g., Bourrier et al. 2020; Shporer et al. 2019; Wong et al. 2020b), we first removed all data marked by a bad quality flag by the SPOC pipeline and then we ran a 16 point-wide moving median filter to the primary and secondary eclipse-masked light curve. We used this median filter to remove 3σ outliers. As was done by other authors, we removed the first 0.5 days of data from each TESS orbit (Wong et al. 2020b,a). The total number of points removed along the whole process was on the order of 25%. This exercise was repeated over both the SAP and PDC light curves. Finally, we divided the light curve of each orbit by its median value, computed from the primary transit and secondary eclipse masked light curves. We did not carry out any other detrending nor normalization strategy (see e.g., Huang et al. 2018; Wang et al. 2019) to avoid removing phase variations that are of relevance to our analysis. Instead, we compensate for any remaining noise by increasing the individual photometric error bars by the amount of correlated noise.

To estimate the amount of correlated noise in both versions of the light curve SAP and PDC, we employed the method of residual binning described in Winn et al. (2008). From the light curve, we subtracted our joint best-fit model, which is described in detail in Sect. 4. We then divided the light curve residuals into M bins of equal duration, each one holding N data points. If the data are not affected by correlated noise, they should follow the expectation of independent random numbers: $$\widehat{{\sigma }_N}={\sigma }_1{N}^{-1/2}{[M/(M-1)]}^{1/2}.$$(1)

In the equation, ${\sigma }_1^2$ is the sample variance of the unbinned data, and $\widehat{{\sigma }_N^2}$ is the sample variance (or RMS) of the binned data: $${\sigma }_N=\sqrt{\frac{1}{M}{\sum }_{i=1}^M{(\langle \widehat{{\mu }_i}\rangle -\widehat{{\mu }_i})}^2},$$(2)

where $\widehat{{\mu }_i}$ is the mean value of the residuals per bin and $\langle \widehat{{\mu }_i}\rangle$ is the mean value of the means. When correlated noise is present in the data, each σ_N will differ by a factor of β_N from their expectation, $\widehat{{\sigma }_N}$. We may account for correlated noise by enlarging uncertainties by a factor β. We compute the average of the β_N’s derived from different time bins relevant to the timescales involved in the study. For instance, when primary transits are fit, relevant timescales are usually close to the ingress or egress time (see e.g., von Essen et al. 2013). In this work, the study spans the whole orbital phase of KELT-1b. As a consequence, we considered the timescale Δt between a fourth of the transit duration to the whole orbital period, divided into steps of a fourth of the transit duration. This corresponds to a total of 84 Δt’s. We then averaged all the β_N ’s computed in this range to obtain the relevant β. Our results are β_SAP = 8.61 and β_PDC = 3.06. These values do not come as a surprise, as the PDC data have their own detrending, while the SAP data do not. In any case, we adopt the PDC data, but with their individual photometric error bars enlarged by β_PDC to alleviate theeffects caused by correlated noise. We show the photometry of KELT-1 in Fig. 1, for both SAP (red points) and PDC (black points) light curves. TESS stared at the target for about 23 days, with the usual gap due to a data downlink.

Fig. 1 Detrended flux of KELT-1 observed by TESS. SAP photometry is indicated with red points and PDC with black circles. Both light curves are shown after our internal data preparation. The fifteen transits analyzed in this work are indicated in the bottom of the figure with red lines.

3 Analysis and model considerations

3.1 Third light contamination

Companions, whether they are visual or bound by gravity, dilute the depth of transits and must be properly accounted for. Due to TESS’s large pixels projected onthe sky (21 × 21 arcsec each), it is very likely that flux from a neighbouring star will fall into the aperture around our target of interest (Fig. 2). PDC corrects the flux for each star falling into KELT-1’s aperture. The amount of contaminated light found in theCROWDSAP keyword of KELT-1’s fits file is 1.181012%. Nonetheless, we carried out our own investigation to compare the results.

In analyzing the target pixel file, we discovered one relatively bright visual companion included in the TESS aperture (Gaia DR2 2881784379713578496, henceforth, C1), and a second one that is much fainter (Gaia DR2 2881784280929656448, henceforth C2). According to Gaia’s Data Release 2 (Gaia Collaboration 2018), C1 has a magnitude of G_C1 = 15.2428 ± 0.0005 and an effective temperature of ~5300 K, while C2 is about two magnitudes fainter, G_C2 = 17.3132 ± 0.0009, and of unknown temperature. For comparison, KELT-1’s Gaia magnitude is G_KELT-1 = 10.5905 ± 0.0003. Additional targets have been found by Gaia DR2 with magnitudes G > 20 and are, consequently, negligible. Furthermore, Siverd et al. (2012) found a companion ~0.5 arcsec away from KELT-1, with a magnitude difference of ΔH = 5.90 ± 0.10, thus, it is negligible as well at TESS wavelengths. Surprisingly, this target was not detected by Gaia. In addition to these, systematic surveys such as Mugrauer (2019); Piskorz et al. (2015); Belokurov et al. (2020) do not list KELT-1 as having another close-in (<2 arcsec) companion.

To compensate for the difference in distance between KELT-1 and C1, C2, we used Pogson’s law: $$m_{\text{KELT-1}}-m_{\text{C1},\text{C2}}=-2.5\times {\log }_{10}\left(\frac{f_{\text{KELT-1}}}{f_{\text{C1},\text{C2}}\times D}\right),$$(3)

where m corresponds to the Gaia magnitudes, f to the corresponding fluxes, and D is a factor that compensates for the difference in distance. To compute the fluxes, f, we made use of PHOENIX stellar intensities (Göttingen 2018) convolved with Gaia’s filter response. Solving this equation, we determined D and we used it to scale down the spectra within the TESS transmission response. As previously mentioned, for C1 there are reliable measurements of the effective temperature of the star of ~5300 K, so we used PHOENIX spectra matching this value. In the case of C1, the third light contribution ascends to 1.74%. However, there is no information about the effective temperature of C2. As a consequence, we tried several PHOENIX spectra with effective temperatures ranging between 3500 and 8500 K, with a step of 1000 K. In each case, we computed the third light contribution in TESS bandpass. As extreme values, we found that for T_eff = 8500 K the third light contribution was 0.17%, while for T_eff = 3500 K it increased to 0.28%. As a consequence, the third light contribution lies between 1.91 and 2.02%. These values are between 60 and 70% larger than the value reported in CROWDSAP. We did not correct the third light contribution, but we pay special attention to the transit depth, as compared to the literature values. While the difference between 1.18% and 2.0% third light contributionmight be measurable for the transit depth, its effect is negligible for the phase curve amplitude and the secondary eclipse depth.

Fig. 2 Third light contamination in the TESS aperture. Top: target Pixel File (TPF) of KELT-1 showing the aperture mask. Three stars are included and completely blurred. Bottom: ~3 × 3 arcmin field of view centered around KELT-1, the brightest star in the field. The mask and the pixels are schematized with white thick and thin lines, respectively. The field of view has been rotated to be aligned with the ecliptic.

3.2 Limb-darkening coefficients

To represent the stellar center-to-limb intensity variability we used a quadratic limb-darkening law: $$\frac{I(\mu )}{I(1)}=1-u_1(1-\mu )-u_2{(1-\mu )}^2,$$(4)

with linear (u₁) and quadratic (u₂) limb-darkening coefficients (LDCs). In this equation, the intensity is normalized to that of the stellar center, and μ = cos(γ), where γ is the angle between the line of sight and the emergent intensity. To compute the limb-darkening coefficients, we made use of the TESS transmission response and angle-dependent, specific intensity spectra taken from PHOENIX library (Göttingen 2018). For the stellar parameters, we used an effective temperature of T_eff = 6500 K, a surface gravity of log (g) = 4.5, and a metallicity of [Fe/H] = 0.00, all fully consistent with KELT-1 values within errors. Similarly to von Essen et al. (2017) and Claret & Bloemen (2011), we neglect the data points between μ = 0 and μ = 0.075, as the intensity drop given by PHOENIX models is too steep and potentially unrealistic. After integrating the PHOENIX angle-dependent spectra convolved with the TESS response, we fit the derived intensities normalized by the maximum values with Eq. (4) with a Markov-chain Monte Carlo (MCMC) approach. The derived limb-darkening coefficients for KELT-1 are u₁ = 0.319(1) and u₂ = 0.227(3). Errors for the coefficients are derived from the posterior distributions of the MCMC chains after visually inspecting them for convergence. In order to assess the quality of our procedure, we fit the LDCs to TESS primary transit light curves. From their posterior distributions, we obtained consistent results with their PHOENIX counterparts. In order to reduce computational cost, in this work we consider u₁ and u₂ as fixed. To double check our procedure, we found our values comparable to those reported by Claret (2017), equal to u₁ = 0.3079 and u₂ = 0.2295.

3.3 Primary transit parameters derived from TESS light curves

To derive updated transit parameters for KELT-1b for the purposes of our joint model as prior information, we fit TESS transits of KELT-1b. For this, we considered the data around primary transits only, plus and minus ~1.5 h of off-transit data to allow for a proper normalization. As a detrending model, we used a second-order time-dependent polynomial. The degree of the polynomial was determined from a prior fit to the data of a first, second, and third degree polynomials, after which we computed the Bayesian Information Criterion (BIC) from the residuals and analyzed which polynomial systematically minimized the BIC. All the transit light curves were fit simultaneously to a primary transit model (Mandel & Agol 2002), while each primary transit light curve had its own set of detrending coefficients. To evaluate whether the non-zero eccentricity of the system, e = 0.0013 ± 0.0005 (Beatty et al. 2019), has any impact on the derived primary transit parameters, we considered two orbits for KELT-1b, a circular one and an eccentric one, the latter with an eccentricity and argument of periastron fixed to the values specified in Beatty et al. (2019). The transit photometry will not allow us to distinguish between the small-eccentricity and a zero-eccentricity model, but we aim to know whether the derived transit parameters differ. The fitting parameters for the primary transit model in the case of a circular orbit are the semi-major axis, a∕R_S, the inclination in degrees, i, the orbital period, P, the planet-to-star radius ratio, R_P∕R_S, and the mid-transit time of reference, T_c. As described in Sect. 3.2, we considered a quadratic limb-darkening law with fixed limb-darkening coefficients. For a non-circular orbit the primary transit parameters are the semi-major axis, a∕R_S, the inclination in degrees, i, the orbital period, P, the planet-to-star radius ratio, R_P ∕R_S, the time of periastron passage, T_peri, the argument of periastron in degrees, ω, and the orbital eccentricity, e. The latter two parameters e and ω remained fixed throughout the fit, whereas the other were treated as free to fit.

To derive the best-fit values for the model parameters and their corresponding uncertainties we used a Markov-chain Monte Carlo (MCMC) approach, all implemented in routines of PyAstronomy¹ (Patil et al. 2010; Jones et al. 2001). We iterated 10⁶ times, with a conservative burn-in of the first 20% of samples. For all the parameters we considered uniform priors around ±50% their corresponding starting values, which were, in turn, taken from Beatty et al. (2019). These are specified in the last column of Table 1. We computed the mean and standard deviation from the posterior distributions and used these as our best-fit values and 1σ uncertainties. We checked the convergence of the chains by visually inspecting each one of them and by dividing the chains in three sub-chains. In each case, we computed the usual statistics and we considered that a chain converged if the derived parameters were consistent within their uncertainties. The best-fit transit parameters for both the circular and non-circular solution is given in Table 1. The posterior distributions and the corresponding correlations between parameters for the circular orbit can be seen in Appendix A.1. Besides the well known correlation between a∕R_S and i, and P and T_c, the remaining parameters are uncorrelated, with Pearson’s correlation values ranging between − 0.05 and 0.05. The very small eccentricity indicated by Beatty et al. (2019) would not affect the derivation of the parameter values, we find all values to be consistent within 1σ uncertainty. The phase-folded light curves for the circular orbit can be seen in Fig. 3.

From TESS photometry, the ephemeris of the system has the following values: $$\begin{array}{lll}{T}_{c,e=0}\hfill & =\hfill & 2458764.31647\pm 0.00018{\text{BJD}}_{\text{TDB}},\hfill \\ P_{e=0}\hfill & =\hfill & 1.2174928\pm 1.7\times {10}^{-6} \text{days}.\hfill \end{array}$$

The period value is in 1σ agreement to the most recent determination of Maciejewski et al. (2018), but less precise due to our shorter time baseline.

3.4 Updated stellar and companion parameters

The stellar radius and stellar mass are varied as part of our modeling (see Sect. 4.2). As a consequence, in this work we have redetermined the stellar and brown dwarf companion radii and masses, taking advantage of the newly available parallax from Gaia DR2 (Gaia Collaboration 2018), together with the available photometry from all-sky broadband catalogs. We used the semi-empirical approach of measuring the stellar spectral energy distribution (SED) described by Stassun et al. (2017) and Stassun & Torres (2016). We pulled the B_T V_T magnitudes from Tycho-2, the BV gri magnitudes from the APASS catalog, the JHK_S magnitudes from 2MASS, the W1–W3 magnitudes from WISE, the GG_BPG_RP magnitudes from Gaia, and the NUV magnitude from GALEX. Together, the available photometry spans the full stellar SED over the wavelength range 0.2–10 μm (see Fig. 4).

We performed a fit using Kurucz stellar atmosphere models, with the priors on effective temperature (T_eff), surface gravity (log g), and metallicity ([Fe/H]) from the spectroscopically derived values reported in the KELT-1b discovery paper (Siverd et al. 2012). The remaining free parameter is the extinction, A_V, which we limited to the maximum for the line of sight from the dust maps of Schlegel et al. (1998). The resulting fit (Fig. 4) has a reduced χ² of 1.6, and best-fit A_V = 0.20 ± 0.04. Integrating the (unreddened) model SED gives the bolometric flux at Earth of F_bol = 1.616 ± 0.038 × 10⁻⁹ erg s⁻¹ cm⁻². Taking the F_bol and T_eff together with the Gaia DR2 parallax, adjusted by + 0.08 mas to account for the systematic offset reported by Stassun & Torres (2018), we determined the stellar radius. From the empirical relations of Torres et al. (2010) we estimated the stellar mass, which is consistent with that inferred from the spectroscopic log g together with the radius (1.36 ± 0.06 M_⊙). Finally, the estimated mass together with the radius gives the mean stellar density. With the updated value for the stellar radius and the planet-to-star radii ratio derived fitting TESS photometry, we report here the KELT-1b radius in the TESS bandpass. All parameters derived here, and specified in Table 2, are in a 1σ agreement to the previous determination of Siverd et al. (2012) and yield a slight improvement in precision.

Fig. 3 Primary transit light curve and residuals for KELT-1b for a circular orbit as a function of the orbital phase. Black points denote TESS data. Continuous red line corresponds to the best-fit primary transit model. The residuals have been artificially shifted below the transits.

Fig. 4 Spectral energy distribution (SEDs) of KELT-1. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the bandpass. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (black).

4 Modeling

Our model consists of six components fit to the TESS data simultaneously. The components are presented in order of decreasing impact, which was determined by first fitting the different components to TESS data separately. To reduce the number of fitting parameters and to carry out a more realistic approach, if a given parameter was present in another model component, the pair was treated as equal. For each model component, we detail the model parameters, we highlight those that are being fit by bold-facing them in their respective equations, and we underline those that are being treated as equal by our MCMC algorithm. In all cases, the variable t corresponds to time as measured and provided by TESS.

We carried out four fitting approaches, all of them considering zero eccentricity. First, we considered KELT-1b mass fixed to the value reported by Beatty et al. (2019) (henceforth, M1), and we used uniform priors in all the parameters except the orbital period and the mid-transit time of reference, which always had Gaussian priors using the starting values and uncertainties reported in Table 1. The difference between M1 and the second fitting approach, M2 is that M2 considers Gaussian priors on all the parameters. Then M3 and M4 are similar to M1 and M2, respectively, with the only difference that the brown dwarf mass is added as a fitting parameter. These four model approaches are intended to investigate whether the choice of priors has any impact on our results and to assess whether the TESS data would confidently allow for a proper mass determination through data fitting.

4.1 Primary transit

We modeled KELT-1b’s primary transits, PT(t), following Mandel & Agol (2002). The fitting parameters are, as specified in Sect. 3.3, the semi-major axis, a∕R_S, the inclination in degrees, i, the orbital period, P, the planet-to-star radius ratio, R_P∕R_S, and the mid-transit time of reference, T_c. Thus, we have: $$\text{PT}(\text{t})=\text{PT}(\underset{\_}{a/R_S},\underset{\_}{i},\underset{\_}{P},\underset{\_}{R_P/R_S},\underset{\_}{T_c}).$$(7)

The amplitude of the variation connected to the primary transit is its own depth, ~5900 ppm. All parameters are fit and considered as equal by other model components.

4.2 Ellipsoidal variation and Doppler beaming

The secondary eclipses presented by Beatty et al. (2019), observed with Spitzer at 3.6 and 4.5 μm, are well placed nearby the maximum of planetary thermal emission. In constrast to this, after phase-folding the TESS KELT-1b data on its orbital period, we found the secondary eclipse to be placed at the minimum level of the phase curve variation. This is because the phase curve variability in TESS photometry is not dominated by the thermal emission of KELT-1b, as it is for the Spitzer light curves of Beatty et al. (2019), but by the ellipsoidal variations in the host star caused by the gravitational pull of the brown dwarf companion. The amplitude of the effect depends on the properties of the stellar surface, the masses of the companion and the star, and the distance between the two bodies (see e.g., Welsh et al. 2010; Mislis & Hodgkin 2012). Even though we considered the eccentricity to be zero, the expression for ellipsoidal variation accounts for eccentric and inclined orbits (Beatty et al. 2019): $$\begin{array}{lll}\Delta F_{\text{EV}}(t)\hfill & =\hfill & A_{\text{EV}}\left(1-\cos \left(\frac{4\pi (t-\underset{\_}{T_c})}{\underset{\_}{P}}\right)\right),\hfill \\ A_{\text{EV}}\hfill & =\hfill & \beta \frac{\underset{\_}{M_P}}{M_{\text{S}}}{\underset{\_}{\left(\frac{R_S}{a}\right)}}^3{\left(\frac{1+e\cos (\nu )}{1-e}\right)}^3\sin {(\underset{\_}{i})}^3,\hfill \\ \beta \hfill & =\hfill & 0.12\frac{(15+u)(1+g)}{(3-u)}.\hfill \end{array}$$

In this equation, M_P and M_S correspond to the companion and stellar masses, respectively, a∕R_S is the relative distance between companion and star scaled to the stellar radius (semi-major axis for a circular orbit), and β (Morris 1985) depends on the linear limb darkening coefficient u and the gravity darkening coefficient g, which were taken from Claret (2017), interpolated to the stellar parameters of KELT-1. Specifically, for KELT-1b, g = 0.2150, and u = 0.4337. The mass of the star was adopted from this work (Sect. 3.4). The fitting parameters are the mid-transit time, T_c, the orbital period, P, and the semi-major axis, a∕R_S. The mass of the KELT-1b was fixed to the value of Beatty et al. (2019) in the modeling approaches M1 and M2, while in M3 and M4 we employed this mass estimation as prior. From literature values, the amplitude of the ellipsoidal variation is A_EV ~ 400 ppm.

As KELT-1b and star orbit around their common barycenter, the star periodically moves towards and away from the observer with a frequency equal to the orbital period of the companion. This creates a variation of the host star brightness ΔF_DB due to Doppler beaming. Its expression (Beatty et al. 2019) is given by: $$\begin{array}{lll}\Delta F_{\text{DB}}(t)\hfill & =\hfill & (3-\alpha )\frac{K_{\text{amp}}(t)}{c},\hfill \end{array}$$(11)

where $$\begin{array}{lll}{K}_{\text{amp}}(t)\hfill & =\hfill & \frac{2\pi a_{*}\sin (\underset{\_}{i})}{\underset{\_}{P}\sqrt{1-e^2}}\left[\cos (\omega +{\nu }_{\text{S}}(\underset{\_}{T_c}))+e\cos (\omega )\right],\hfill \\ a_{*}\hfill & =\hfill & \left(\frac{\underset{\_}{M_P}}{M_{\text{S}}}\right)\underset{\_}{\left(\frac{a}{R_S}\right)}{R}_{\text{S}}.\hfill \end{array}$$

Here, ω is the argument of the periastron, ν_S is the true anomaly of KELT-1, measured with respect to the mid-transit time, T₀, and R_S is specified in Sect. 3.4. The speed of light is represented by c and α depends on the central wavelength of TESS bandpass and the effective temperature of the star (Loeb & Gaudi 2003). The fitting parameters are a∕R_S, T_c, P, M_P (only in M3 and M4), and i. The amplitude of the Doppler beaming is K_amp ~ 90 ppm, based on the mass of KELT-1b derived in Beatty et al. (2019). It differs from the ellipsoidal variation not only in amplitude, but also in frequency (one repetition per orbital period for the Doppler beaming, two repetitions per orbital period for the ellipsoidal variation). Like the ellipsoidal variations, Doppler beaming also strongly depends on the mass of the companion, hence we vary them together in this description.

4.3 Reflected light, thermal emission and secondary eclipse

Following Cowan & Agol (2008), we model the reflection of starlight and the thermal emission of KELT-1b as a series of first order expansions in sines and cosines: $$\text{PPV}(\text{t})=c_0+c_1\times \sin (\text{2}\pi \text{t/}\underset{\_}{P})+c_2\times \cos (\text{2}\pi \text{t/}\underset{\_}{P}).$$(14)

PPV(t) stands for planetary phase variability. The fitting parameters of this model component are the vertical offset, c₀, the amplitudes of the sine and cosine, c₁ and c₂, respectively,and the orbital period, P. This linear combination of sines and cosines allow for a potential offset between theregion of maximum brightness and the sub-stellar point. The amplitude of the phase curve is around ~100 ppm. During secondary eclipse KELT-1b is hidden behind the star. As a consequence, the thermal emission and reflected light components should not show up during secondary eclipse.

To model the secondary eclipse we make use of a step function (Sackett 1999), in sync with the levels of PPV(t) right before and after secondary eclipse begins and ends, respectively, to allow for continuity. The shape of the step function is limited by the first to fourth contacts, which are in turn determined from the center of the eclipse, CE: $$\text{CE}=\frac{\underset{\_}{P}}{2}\times (1+4 \text{e} \cos (\omega )),$$(15)

with the transit duration between first and fourth contacts: $$T_{\text{dur},\text{1-4}}=\frac{\underset{\_}{P}}{\pi }{\sin }^{-1}\left(\frac{\sqrt{{(1+\underset{\_}{R_P/R_s})}^2-{(\underset{\_}{a/R_S} \cos (\underset{\_}{i}))}^2}}{\underset{\_}{a/R_S}}\right),$$(16)

and the transit duration between second and third contacts: $$T_{\text{dur},\text{2-3}}=\frac{\underset{\_}{P}}{\pi }{\text{sin}}^{-1}\left(\frac{\sqrt{{(1-\underset{\_}{R_P/R_S})}^2-{(\underset{\_}{a/R_S} \cos (\underset{\_}{i}))}^2}}{\underset{\_}{a/R_S}}\right).$$(17)

As the secondary eclipse is not necessarily symmetrically centered with the phase curve, the ingress and egress parts of the step function canhave different levels. As a consequence, the eclipse depth ED is not measured as the difference of the out-of-eclipse flux to the in-eclipse flux, but as the difference between the point in the phase curve that corresponds to the center of the eclipse (as if the eclipse was not there), and the in-eclipse flux (see e.g., Bell et al. 2019, their Fig. 4). The secondary eclipse model depends on the following parameters: $$\text{SE}(t)=\text{SE}(\underset{\_}{a/R_S},\underset{\_}{i},\underset{\_}{P},\underset{\_}{R_P/R_S},\text{e},\omega ,ED).$$(18)

Assuming pure thermal emission and thus no reflected light, using the brightness temperature determined by Beatty et al. (2019) at Spitzer bandpasses in the near IR, the expected secondary eclipse depth in the TESS bandpass is ~300 ppm.

4.4 Stellar activity

After subtracting the primary transit, ellipsoidal variation and Doppler beaming models from the data using literature values, we hoped we would only be left with the reflected light, the thermal emission, and the secondary eclipse of KELT-1b. Visually inspecting the individual orbits, we noticed the shape of the phase curve and the eclipse depth were evidently changing from orbit to orbit. Before concluding that KELT-1b may be exhibiting weather patterns in its atmosphere (see e.g., Armstrong et al. 2016), we addressed other possible sources that could give rise to such variability. One of the sources investigated was stellar activity. As specified in Siverd et al. (2012), the host is an F-type star on the main sequence with a rather shallow convective envelope. As a result, the star is expected to have a low level of activity or none at all.

In order to investigate the presence of spot modulation in the TESS data, we computed a Lomb-Scargle periodogram Lomb (1976); Scargle (1982); Zechmeister & Kürster (2009) of TESS residuals, after the previously mentioned model components of Sects. 4.1–4.3 were fit initially and removed. The resulting periodogram can be found in the bottom panel of Fig. 5. The figure shows unresolved peaks around 0.7 cycles per day (c/d), which could be due to rotational modulation induced by spots. Such a signal usually results in a peak of power around the rotational frequency of the star and will only show resolved peaks if the stellar spots are long-lived and do not migrate in latitude. As a consequence, we could not simply remove the peaks using traditional frequency fitting.

To determine the mean frequency of what appears to be rotational modulation of the host star, we used the phase dispersion minimization technique instead (PDM, Stellingwerf 1978). We divided the frequency space between 0.4 and 0.9 c/d, where most of the power lies, in steps of 1/2Δt, where Δt corresponds to the total time of TESS observations. Then, we converted these frequencies into periods (so-called trial periods), and with each one of these trial periods we phase-folded TESS residuals. We divided the orbital phase into 50 bins, and we computed the overall variance of all the samples per phase-folded dataset (Eq. (2) in Stellingwerf 1978). We compared this value to the variance of the residuals. If the trial period minimizes the dispersion of the phase-folded data, then the PDM should be minimized. The top panel of Fig. 5 shows the overall variances for the different trial periods, converted back to frequencies. The minimum corresponds to ν_min,spots = 0.659 cd⁻¹, slightly rightward of the highest peak. This results was consistent even after changing the frequency range where the PDM was evaluated, and the number of bins into which we divided the phase-folded data. To estimate an uncertainty for this value we carried out a fit of the spectral power with a Gaussian function. To derive the best-fit values, we carried out a least-squared minimization between data and model. For the Gaussian function, we considered the mean value as fixed, and equal to ν_min,spots. The fitting parameters are the amplitude, the offset, and the standard deviation. From the latter, we computed the error of ν_min,spots to be equal to 0.107 cd⁻¹. The fit is shown in Fig. 5 with a blue continuous line. If this periodicity is interpreted as the stellar rotation period with P = 1.52 ± 0.29 d, we find an agreement within 1σ to the stellar period derived from the measured projected rotation velocity of the star of P = 1.33 ± 0.06 d (Siverd et al. 2012).

To represent the observed stellar activity (SA) of KELT-1 we modeled it as follows: $$\text{SA}(t)=A_{spots}\times \sin [2\pi (t/{\nu }_{\text{min},\text{spots}}+{\varphi }_{spots})],$$(19)

where A_spots corresponds to the amplitude of the photometric variability due to stellar spots, ν_min,spots corresponds to the frequency minimizing the PDM, and ϕ_spots is the phase, in units of 2π. As the periodogram shows, the amplitude of this effect is A ~ 210 ppm. The fitting parameters are the amplitude and the phase. We do not fit for the frequency, as we are limited by the total amount of time of TESS on the target.

Fig. 5 Periodicity of the residual light curve modulation. Top: phase dispersion minimization evaluated in the frequency range between 0.4 and 0.9 c/d. Bottom: Lomb-Scargle periodogram (amplitude in ppm versus frequency in cd−1) for the residuals of KELT-1, showing only the frequency range where most of the power lies. Overplotted in blue is the fit of a Gaussian function over the power to estimate an error for νmin. The vertical orange line indicates the measured projected rotation velocity of the star as determined by Siverd et al. (2012). For comparison, the frequency of the orbital motion of KELT-1b is 0.821 c/d.

4.5 Combined model and fitting approach

Our combined model is the addition of the six components previously described, PT(t) + ΔF_EV(t) + ΔF_DB(t) + SA(t) + PPV(t) + SE(t). For M1 and M2, the total number of fitting parameters are 11, while for M3 and M4 they are 12, when the mass of the brown dwarf is included. We note that the vertical offset is arbitrarily included in the PPV(t). As our combined model is a simple linear combination of the six components, it is irrelevant where the offset is. This combined model was fit to the unbinned TESS photometry. For each one of the model approaches, we determined the BIC to determine which fitting approach was the best at reproducing the data.

As in the primary transit fitting approach, to derive the best-fit values for the parameters, we made use of the MCMC. In each one of the previously mentioned cases we iterated 10 000 times, with a burn-in of the first 2000 samples. This burn-in was determined from a posterior visual inspection of the chains, where we also investigated their convergence. As usual, the best-fit values for the parameters, along with their corresponding 1σ uncertainties, are derived from the mean and standard deviation of the posterior distributions.

5 Results

5.1 Parameters derived for KELT-1b

Table 3 shows the parameters obtained from our combined fit, along with those determined from these values. The corresponding best-fit values for the derived parameters, along with their uncertainties, were computed in the usual way from their posterior distributions. TESS phase-folded data, along with our best-fit combined model, can be seen in Fig. 6 for M2, which is the model choice with the lowest BIC value. The figure shows TESS photometry of KELT-1bin black points, phase-folded with the best-fit orbital period. Primary transits are close to phases ϕ = 0 and 1, while the secondary eclipses lie around ϕ = 0.5. The first three panels show TESS unbinned data, while the last two are binned each Δϕ = 0.02 (equivalently, Δt ~ 30 minutes). At this cadence, the photometric precision is 80 ppm.

From the coefficients connected to the PPV(t) component we determined the amplitude of the phase curve: $$\text{A}({\text{c}}_{\text{1}},{\text{c}}_{\text{2}})\text{=}\sqrt{{\text{c}}_{\text{1}}^{\text{2}}{\text{+c}}_{\text{2}}^{\text{2}}}.$$(20)

In addition, we computed the offset between the region of maximum brightness and the substellar point: $${\varphi }_{\text{off}}={\tan }^{-1}(c_1/c_2).$$(21)

The derived values for the different fitting approaches can be seen in Table 3.

5.2 Physical parameters derived from these observations

We determined the equilibrium temperature at the sub-stellar point on KELT-1b (i.e., the irradiation temperature), $$T_0=T_{\text{eff}}\sqrt{R_S/a},$$(22)

where T_eff is the effective temperatureof the star and a∕R_S is the semi-major axis scaled by the stellar radius obtained from the primary transit fitting component of our joint model (Mandel & Agol 2002). From this parameter, we can derive the effective temperature of the day side in the no-albedo, no-circulation limit, $T_{ϵ=0}={(2/3)}^{1/4}{T}_0$ (Burrows et al. 2008; Cowan & Agol 2011), the maximum possible day-side temperature of KELT-1b.

To compute the brightness temperature, T_b(λ), on the day and the night side of KELT-1b, we follow Cowan & Agol (2011), $$T_b(\lambda )=\frac{hc}{\lambda k}{\left[\log \left(1+\frac{\exp (hc/\lambda k{T}_{\text{*b}}(\lambda ))-1}{\psi (\lambda )}\right)\right]}^{-1}.$$(23)

Here, T_*b(λ) is the star’s brightness temperaturein the TESS bandpass, and h, c and k are the Planck constant, the speed of light, and the Boltzmann constant, respectively. For KELT-1, the brightness temperature is T_*b = 6456 K. In the equation, ψ(λ) corresponds to the ratio of the companion’s day- or night-side intensity to the stellar intensity at that wavelength. In particular, for the day-side temperature we made use of the secondary eclipse depth divided by the primary transit depth, as this is a direct measure of the ratio of the companion’s day-side intensity to the stellar intensity, $$\psi {(\lambda )}_{\text{day}}=ED/{(R_P/R_S)}^2.$$(24)

To compute the night-side temperature we used the difference between the secondary eclipse depth and the phase variation amplitude, combined with the offset between the peak of the phase curve and secondary eclipse: $$\psi {(\lambda )}_{\text{night}}=\frac{ED-2A\cos {\varphi }_{\text{off}}}{{(R_P/R_S)}^2}.$$(25)

We then computed the brightness temperatures of the day and the night side of KELT-1b by substituting these two values into Eq. (23). The lowest point of the phase curve in Fig. 7 shows a slightly higher flux level than the bottom of the secondary eclipse, which suggests a night-side emission from KELT-1b in the TESS bandpass. However, this estimated amount of night-side flux is significant to less than 2σ, therefore, the derived night-side temperature T_night of about 1400 K is not well constrained either.

We obtained uncertainties on all these parameters using the posterior probability distributions for the values c₁, c₂, ED and a∕R_S. In particular, for each one of the 8000 MCMC iterations, we computed each one of the previously specified equations. In this way, the values reported in Table 3 have their best-fit values and uncertainties determined from their mean and standard deviations.

Table 3 summarizes the parameters derived in this work and compares them to those computed by Beatty et al. (2019). To begin with, the offset between the peak of the phase curve and secondary eclipse is consistent within 1σ to the value reported by Beatty et al. (2019), regardless of the fitting approach we used. This suggests that the TESS phase curve is sensitive to thermal emission of KELT-1b, as is Spitzer. However, the day-side brightness temperature inferred from TESS data is higher than this reported by Beatty et al. (2019). The difference may be attributed to stronger molecular absorption in the Spitzer IRAC bands compared to that in the TESS band, as discussed in Sect. 6.4.

The day-side brightness temperature in the TESS bandpass is marginally larger than T_ϵ=0², even though it still consistent at the 1σ level. This might suggest that either: (a) the two measurements are indeed consistent, and we can derive reliable Bond albedo and heat redistribution efficiency for KELT-1b (Sect. 6.1); (b) reflected light contributes to the large TESS eclipse depth (Sect. 6.2); or (c) an emission feature in the TESS bandpass increases the brightness temperature above the effective temperature. Here, the TESS data could be probing the deeper hotter layers of the atmosphere compared to the Spitzer bands, which could have strong opacity and, hence, probe cooler upper regions of the atmosphere (Sect. 6.4). We discuss these options in the following sections.

Fig. 6 Data (black points) versus best-fit model (red continuous line) as a function of the orbital phase of KELT-1b. From top to bottom: full data set; same data, but with the best-fit primary transit, PT(t), removed; ellipsoidal variation, EV(t), removed; Doppler beaming and stellar activity, DB(t) and SA(t), respectively, removed, at which point the secondary eclipse and the phase curve are visible; residuals, once the secondary eclipse and the phase curve are removed. The last two panels are shown at a 0.02 orbital phase binning. The last four panels are given in ppm.

Fig. 7 Relative flux of TESS data binned as specified before, showing only the phase curve and secondary eclipse components. In red and blue bars, we indicate the day- and night-side fluxes.

6 Discussion

6.1 Bond albedo and heat redistribution efficiency through a self-consistent analysis combining Spitzer & TESS

To derive the Bond albedo and the heat redistribution efficiency for KELT-1b, we assume that the flux of KELT-1b measured by TESS is solely composed of thermal emission. Furthermore, we make the simplifying assumption that the atmosphere of KELT-1b behaves entirely as a blackbody. In other words, we can infer a Bond albedo that neglects the fact that opacities might be caused by molecular absorption. Potential deviations from this blackbody assumption are considered in Sect. 6.4. We estimate the day-side temperature using the error-weighted mean of the Spitzer and TESS brightness temperatures, that is, by combining the Beatty et al. (2019) values specified in the last column of Table 3and our values in the TESS bandpass for M2, while for the night-side effective temperature we only average the two Spitzer measurements because of the unconstrained TESS value. The resulting day and night effective temperatures are listed in Table 4. The averaged day-side temperature is close to the value of Beatty et al. (2019), since their more precise Spitzer results dominate the error-weighted mean. It is worth mentioning that T_day is below T_ϵ=0.

We then follow Cowan & Agol (2011) to determine the Bond albedo: $$A_B=1-\frac{5T_{\text{night}}^4+3T_{\text{day}}^4}{2T_0^4},$$(26)

and the heat redistribution efficiency, $$ϵ=\frac{8}{5+3{(T_{\text{day}}/T_{\text{night}})}^4}.$$(27)

The resulting parameters are also listed in Table 4.

6.2 Reflected light to explain the deep TESS eclipse

Since the amount of thermally emitted light by KELT-1b is reduced at the optical TESS bandpass compared to the Spitzer bandpasses, a potential contribution of reflected light to the secondary eclipse depth would have a considerable impact on the calculation of the brightness temperature. Therefore, it is tempting to attribute the high TESS day-side brightness temperature to reflected light, which would not have significantly affected the Spitzer measurements. To investigate this possibility, we then model the TESS eclipse depth following Alonso (2018): $$\frac{F_P}{F_S}=A_g{\left(\frac{R_P}{a}\right)}^2\varphi (\alpha )+\frac{B(\lambda ,T_{d,p})}{B(\lambda ,T_s)}{\left(\frac{R_P}{R_S}\right)}^2,$$(28)

where the first term of the right-hand side of the equation corresponds to reflected light and the second term is thermal emission. For the reflected light term, A_g is the geometric albedo, R_P and R_S are the planet and stellar radius, respectively, a is the semi-major axis, and ϕ(α) is the phase function that can be approximated to 1 during secondary eclipse. For the thermal emission term, B(λ, T_d,p) and B(λ, T_s) correspond to the blackbody emission of the companion and star, respectively, at brightness temperatures of T_d,p and T_s. If we assume the brightness temperature in the TESS bandpass to be the same as measured in the Spitzer bandpasses, T_d,p = 2988 ± 60 K (Beatty et al. 2019), then using our derived values for the TESS eclipse depth, the planet-to-star radii ratio, and the semi-major axis scaled to the stellar radius, we estimate a geometric albedo of 25 ± 11%. Hot Jupiter exoplanets have been observed with a range of geometric albedos from low values close to zero (Rowe et al. 2008; Esteves et al. 2015; Mallonn et al. 2019) to values of up to 0.3 and more (Demory et al. 2011; Keating & Cowan 2017; Wong et al. 2020c). A large geometric albedo at optical wavelengths has for example been observed for Kepler-7b with A_g = 32 ± 3% (Demory et al. 2011). The authors interpreted this large geometric albedo as caused by Rayleigh scattering.

Our day and night temperatures satisfy $T_{\text{day}}^4+T_{\text{night}}^4<\frac{2}{3}{T}_0^4$. This means that KELT-1b has a non-zero Bond albedo, if the day-side and night-side atmospheres are treated as blackbodies (Cowan & Agol 2011). In order to evaluate whether the Bond and geometric albedos are self-consistent, we estimate the fraction of incident stellar light in the TESS bandpass: we convolved and afterwards integrated the TESS transmission response with PHOENIX stellar intensity models for T_eff= 6500 K, log (g) = 4.5 and [Fe/H] = 0, matching the values of KELT-1 within uncertainties. The fraction of stellar flux in the TESS bandpass is F_TESS = 31%. The minimum Bond albedo, A_B,min should be A_g × F_TESS = 7.6% when considering the geometric albedo specified in Sect. 6.2. This value is compatible with the Bond albedo determined in Sect. 6.1 (Table 4), but requires the geometric albedo to be nearly zero outside of the TESS bandpass.

Even though this analysis of the Bond and geometric albedo of the day side is self-consistent, the assumption of pure blackbody radiation with the day-side brightness temperature being the same at TESS and Spitzer bandpasses might be overly simplified. The depth of the secondary eclipse at optical wavelength could also be explained by the TESS band probing hotter layers of the atmosphere compared to those probed by Spitzer, due to a non-isothermal temperature profile and wavelength-dependent opacity. This last scenario is first motivated (Sect. 6.3) and investigated (Sect. 6.4) in the following subsections.

6.3 Initial estimation of the effective temperature of KELT-1b

To determine the effective temperature of KELT-1b, we assume at this point that KELT-1b’s emission can be approximated by a blackbody. Under this assumption, we fit all the literature values reported by Beatty et al. (2014, 2017, 2019) and Croll et al. (2015) to our new measurement to synthetic eclipse depths. These values were determined by first computing the flux ratio between KELT-1b, represented by a blackbody function of different temperatures, and PHOENIX intensities, following the stellar parameters specified before, to represent KELT-1. This flux ratio was then integrated within all the transmission functions for which literature data were obtained. For the blackbody we created a grid of temperatures between 2000 and 4000 K, with a step of 1 K. For each one of these temperatures we computed the previously mentioned ratio, after which we determined the χ² between the observations and our model. From a χ² minimization, we obtained T_eff = 3010 ± 78 K. The error on the temperature was computed considering Δχ² = 1.

At first glance, Fig. 8 reveals that both the Beatty et al. (2017) and Croll et al. (2015) eclipse depth measurements are significantly higher compared to the derived effective temperature. Both measurements, along with the z-band measurement reported by Beatty et al. (2014), share the same observing technique: they were derived from incomplete orbital phase coverage, solely from secondary eclipse observations. As the mass of KELT-1b is particularly high, the ellipsoidal variation of the star creates a convex secondary eclipse. If the ellipsoidal variation is not removed, then the measured secondary eclipse, due to this convexity, will be deeper (see Sect. 4.4 of Bell et al. 2019). To quantify how much deeper, we carried out the following exercise, focusing on the H-band data reported by Beatty et al. (2017). First, we created synthetic secondary eclipse data. We used the parameters, cadence, and noise amplitude from Beatty et al. (2017). However, for this exercise, thenoise was purely white. The synthetic data include the secondary eclipse, some thermal emission, the ellipsoidal variation, and the Doppler beaming. Once the synthetic time stamps, data and errors were generated, we fit them. As a model, in this case, we used a second-order time-dependent polynomial and a step function with a flat bottom as the secondary eclipse. The detrending should take into account every effect introduced in generating the synthetic data. We fit the data and model using least square minimization. The fitting parameters featured the detrending coefficients and the secondary eclipse depth. We repeated the same process 1000 times. Within these repetitions, we simply counted how many times was the derived secondary eclipse depth larger than the one we usedto create the synthetic data, equal to the value reported in Beatty et al. (2017) for the H band. Regardless the number of repetitions we carried out, we always obtained the same results. About 80% of the time, the fitted secondary eclipse depth was larger than the input one. As a result, the fit is clearly opting for a deeper eclipse to compensate for the convexity. In addition, the percentage at which the secondary eclipse depth was overestimated is around 10%. In other words, by not taking into account the gravitational pull of KELT-1b induced over the star, and the effect on the secondary eclipse light curve, the secondary eclipse depth is enlarged by about 10%.

When determining the effective temperature of KELT-1b without considering the H- and K- band measurements, the effective temperature decreases to T_eff = 2946 ± 71 K, and ${\chi }_{\text{red}}^2$ decreases from 6.8 to 1.6. As we believe the eclipse depth measurements obtained from incomplete phase curve analysis are overestimated, in order to further carry out a more detailed analysis of the atmospheric properties of KELT-1b, we go on to focus on the TESS and Spitzer measurements, both of which benefit from full-orbit monitoring.

Fig. 8 Eclipse depth, in ppm, as a function of wavelength. Literature values are taken from Croll et al. (2015); Beatty et al. (2014, 2017, 2019). Our new data point is plotted in red. The gray area shows a 1σ contour of the equilibrium temperature of KELT-1b, determined from all data. Transmission responses are plotted in black continuous lines, with the exception of TESS, highlighted in red.

6.4 Self-consistent atmospheric models

In this section, we self-consistently model the day- and night-side spectra of KELT-1 b using the GENESIS atmospheric model (Gandhi & Madhusudhan 2017). GENESIS calculates full, line-by-line radiative transfer under radiative-convective, local thermodynamic, hydrostatic, and thermochemical equilibrium. The model self-consistently calculates the atmospheric temperature profile, thermal emission spectrum, and chemical profiles given the incident irradiation, internal flux, and elemental abundances (via metallicity and C/O ratio). Equilibrium chemical abundances are calculated using the HSC CHEMISTRY (version 8) software package, as in Piette et al. (2020). HSC CHEMISTRY uses the GIBBS solver (White et al. 1958) to calculate chemical abundances by minimising the total Gibbs free energy of the system and includes effects due to thermal dissociation or ionisation, which are relevant with regard to to the highly-irradiated day side of KELT-1 b.

While the equilibrium chemistry calculation includes > 150 molecular, atomic, and ionic species, the species which typically dominate irradiated giant planet atmospheres are H₂O, CH₄, CO, CO₂, NH₃, HCN, C₂ H₂, Na, K, TiO, VO, and H⁻, besides H₂ and He (Burrows & Sharp 1999; Madhusudhan et al. 2016; Arcangeli et al. 2018). We therefore include opacity due to these species in our atmospheric models, using cross-sections calculated as in Gandhi & Madhusudhan (2017) and Gandhi et al. (2020) as well as line transition data from a variety of sources (see also Piette et al. 2020), which include the HITEMP, HITRAN, and ExoMol databases (H₂O, CO and CO₂: Rothman et al. 2010, CH₄: Yurchenko et al. 2013; Yurchenko & Tennyson 2014, C₂ H₂: Rothman et al. 2013; Gordon et al. 2017, NH₃: Yurchenko et al. 2011, HCN: Harris et al. 2006; Barber et al. 2014, TiO: McKemmish et al. 2019, VO: McKemmish et al. 2016, H₂ -H₂ and H₂-He Collision-Induced Absorption: Richard et al. 2012). The Na and K opacities are calculated as in Burrows & Volobuyev (2003) and Gandhi & Madhusudhan (2017), while the H⁻ bound-free and free-free cross sections are calculated using the prescriptions of Bell & Berrington (1987) and John (1988) (see also Arcangeli et al. 2018; Parmentier & Crossfield 2018; Gandhi et al. 2020).

We consider the stellar irradiation with a flux distribution factor on the day side of f =2∕3 (e.g., Burrows et al. 2008). Given the extreme irradiation received by KELT-1 b, a low day-night energy redistribution efficiency is expected (Cowan & Agol 2011). Therefore, we chose to explore models in the limit of no day-night energy redistribution for simplicity. We note, however, that recent studies find that H₂ dissociation and recombination may increase heat transport in highly-irradiated atmospheres (Bell & Cowan 2018; Tan & Komacek 2019; Mansfield et al. 2020). We also account for the internal heat flux emanating from KELT-1b, which is characterized by the internal temperature, T_int. In the absence of incident irradiation, for instance, on the night side, T_int is equivalentto the effective temperature used in stellar contexts. While the day-side temperature structure and spectrum are dominated by the irradiation, the night side is dominated by the internal flux and energy advected from the day side. Given the small number of data points available, we do not perform a full atmospheric retrieval which would have more free parameters than data points.Rather, we report physically plausible self-consistent models that provide a consistent explanation for both theday side and night side data. In order to fit the observed data, we explore models with different T_int, metallicities and C/O ratios, as described below.

We are able to fit the day side data using irradiated models with no redistribution (upper panel of Fig. 9). At the temperatures present in the day-side atmosphere, CO is expected to be a significant opacity source in the IRAC 2 band (Burrows & Sharp 1999; Madhusudhan et al. 2016). However, since the brightness temperature in the IRAC 2band is comparable to the IRAC 1 band and lower than the TESS band, this suggests the absence of a thermal inversion in the infrared photosphere, which is consistent with the findings of previous studies on KELT-1 b (Beatty et al. 2017). A thermal inversion would instead cause a strong CO emission feature in the IRAC 2 band, which is not seen. In turn, a non-inverted temperature profile corresponds to a low optical/infrared opacity ratio (e.g., Hubeny et al. 2003; Guillot 2010), for example, due to a depletion of optical absorbers such as TiO and VO. Such a depletion may be caused by the gravitational settling of TiO and VO due to their high mean-molecular weight (Spiegel et al. 2009) or increased efficiency of cold trap processes (Parmentier et al. 2013; Beatty et al. 2017). We find that a model with solar metallicity and C/O = 0.5 provides a good fit to the data if the abundances of TiO and VO are just 5% of those predicted under thermochemical equilibrium (within 1.4, 0.5, 0.8σ of the TESS, IRAC 1, and IRAC 2 points, respectively, see red spectrum in the top panel of Fig. 9). The continuum of the model spectrum is close to a blackbody spectrum at ~3000 K, with the decreasing temperature gradient resulting in a slight CO absorption feature in the IRAC 2 band.

An alternative scenario that involves depleted TiO and VO is that featuring a high C/O ratio (Madhusudhan et al. 2011). We find that a model with solar metallicity in all elements except carbon and C/O = 1 is able to fit the TESS, IRAC 1, and IRAC 2 data within 1.5, 1.7, and 1.9σ, respectively (blue spectrum in the top panel of Fig. 9). The photosphere in the ~0.6–5 μm range lies between ~0.1 and 10 bar, where the temperature profile is relatively isothermal. As a result, the emergent thermal emission spectrum is also close to a blackbody at ~3000 K. In this model, the atmospheric chemistry is still dominated by CO, which results in some absorption in the IRAC 2 band. We note that given the deep photosphere in KELT-1 b, thermal dissociation is less prevalent in the photospherecompared to lower-gravity hot Jupiters with shallower photospheres. We also note that, given the strong stellar irradiation, the day-side spectrum of KELT-1 b is not sensitive to internal temperatures below ≳ 2000 K. As a result, T_int is constrained only by the night side data and we assume a value of T_int = 1500 K for the day side model.

The day side data can be explained by thermal emission alone, without the need for a high albedo, although the latter cannot be ruled out with current data. Figure 9 also shows the T_int = 1000 K cloudy model from Beatty et al. (2020) (from their Fig. 6) for reference. Our clear model is able to fit the TESS point slightly better than the Beatty et al. (2020) cloudy model, while the fit to the Spitzer data cannot be compared as their model is truncated at 2.5 μm. Since the current data do not rule out a cloud-free atmosphere, further observations are needed to constrain whether clouds are present on the day side of KELT-1 b. Indeed, a plausible solution may be a partially cloudy atmosphere, with contributions from both clear and cloudy regions. Future spectroscopic secondary eclipse or phase curve observations may provide more confident constraints on the presence of clouds in KELT-1 b’s day side.

We find that the Spitzer night-side brightness temperatures are consistent with a blackbody in the range ~1000–1500 K (within 1σ). Within the limit of no day-night redistribution, the data can be explained by non-irradiated self-consistent models with T_int in the range ~1500–2000 K (lower panel of Fig. 9). This places an upper limit of ~2000 K on T_int, as a smaller internal flux is permitted by the data if additional flux is advected from the day side to the night side. In the models shown in Fig. 9, the absorption features in the 3.6 and 4.5 μm bands are due to CH₄ and CO, respectively. We assume a C/O ratio of 0.5, consistent with the day side model shown in red in the top panel of Fig. 9. However, we note that the night side data are not able to constrain chemical abundances and a range of C/O ratios and metallicities are consistent with the night side observations. Furthermore, we note that since these models do not include energy advected from the day side, solutions with lower T_int and greater day-night energy redistribution are also possible.

Fig. 9 Self-consistent modeling of the atmosphere of KELT-1b. Top: model day-side spectra (left) and pressure-temperature (P-T) profiles (right), assuming solar metallicity and an internal temperature of 1500 K. The red spectrum and P-T profile correspond to a C/O=0.5 and depleted TiO/VO (e.g., due to gravitational settling), while the blue model corresponds to a C/O = 1. The black dashed lines correspond to blackbody spectra at 3000 and 3200 K. The cloudy, Tint = 1000 K model from Beatty et al. (2020) is shown in grey for reference. Bottom: model night-side spectra (left) and pressure-temperature profiles (right) for KELT-1 b assuming internal temperatures of 1500 K (red lines) and 2000 K (purple lines). Dashed black lines correspond to blackbody spectra at 750, 1000 and 1500 K. Flux is shown in SI units (W sr−1 m−3). In both panels, data are shown by black points and error bars. Coloured circles show the binned model spectra. Vertical gray bands show the TESS and Spitzer bandpasses.

6.5 Tidal spin-up of the host star

Siverd et al. (2012) formulated the suggestion of a complete synchronization of the KELT-1 brown dwarf-star system, including the stellar rotation spin-up by tides to a value equal to the orbital period of KELT-1b. This suggestion was backed by the theoretical calculation of a short synchronization timescale (Siverd et al. 2012; Matsumura et al. 2010) and the large value of measured stellar rotation vsini of 56 km s⁻¹ (Siverd et al. 2012). In this work, we derived a value for the stellar variability period of P_spots = 1.52 ± 0.29 d (Sect. 4.4). When interpreted as the stellar rotation period, it confirms a significant tidal spin-up of the host star, since field stars of mid F spectral type are expected to rotate with a period between 10 and 15 days at the age of KELT-1 of 1.5–2 Gyr (Barnes 2010; Barnes et al. 2016). Despite evidence of being spun-up, our measurement of stellar variability indicates a slightly longer stellar spin period than the brown dwarf orbital period; thus the process of stellar rotation spin-up might still be ongoing until it reaches full synchronization. However, the uncertainty on the stellar rotational period is only 1σ away from the period of KELT-1b, so they might actually be synchronized at present. We can only speculate if the slight discrepancy between stellar variability period and orbital period might also be caused either by systematic features in the PDC photometry, by difficulties in measuring stellar rotation periods from irregular and low-amplitude photometric variability (Shapiro et al. 2020), or by stellar differential rotation. We note that the independent study of Beatty et al. (2020), in analyzing the same TESS data set modeled the photometric variation (which we interpreted as rotational modulation) differently. The authors considered this variation as correlated noise in the TESS data and approximated it with time-dependent polynomials on the order up to seven in different segments of the data. Because of the low amplitude of the signal of 0.02% and potential deviations from strict periodicity due to a time evolution of magnetic surface elements, it is not a trivial task to distinguish whether the source of the flux variation is astrophysical or instrumental in nature. However, the derived planet phase curve parameters are in decent agreement between Beatty et al. (2020) and our work; thus, they appear to be robust against different versions of modeling of the additional flux variation present in the data.

6.6 Agreement of transit depth between optical and near-IR wavelengths

The limb darkening corrected transit depth, ${(R_p/R_s)}^2$ at optical wavelengths of this work agrees with the joint fit of the Spitzer 3.6 and 4.5 μm data of Beatty et al. (2019) to 0.2 ± 1.0%. The investigation of a potential wavelength dependence of ${(R_p/R_s)}^2$ is known as transmission spectroscopy in exoplanet science and is used to characterize the atmospheres of transiting planets. For favorable targets, the transit depth has been found to vary over wavelength by 5% and more (von Essen et al. 2019; Sing et al. 2016). However, KELT-1b is not suited to the atmospheric characterization by transmission spectroscopy because of its large surface gravity and subsequent small value of atmospheric scale height. A typical signal in the transmission spectrum of five scale heights amplitude amounts to a variation in ${(R_p/R_s)}^2$ of only 0.2%, which therefore would not be detectable for KELT-1b with the current instruments.

However, the sub-percent agreement of ${(R_p/R_s)}^2$ between the optical and near-IR is helpful for this study because it indicates that our analysis does not miss a significant source of third light in the large TESS aperture, since the modification of the parameters of KELT-1b by third light shows usually a strong wavelength dependence (Mallonn & Strassmeier 2016; von Essen et al. 2020b). The uncertainty of 1% does not allow to distinguish between the slightly different values of third light contamination applied in the PDC TESS photometry and which we derived independently (see Sect. 3.1). Their sub-percent difference turns out to be not significant for the currently achieved measurement precision. Also, the agreement of ${(R_p/R_s)}^2$ over wavelength rules out a significant color-dependent modification of the transit light curve by temperature inhomogeneities on the stellar surface (e.g., Oshagh et al. 2014; Mallonn et al. 2018) and confirms the rather low activity level of the star found by the low amplitude of photometric stellar variability in Sect. 4.4. In a recent study on hot Jupiter transmission spectroscopy, Alexoudi et al. (2018) described the emergence of an apparent wavelength-variation of ${(R_p/R_s)}^2$ when the impact factor b applied in the fit deviates from the correct one. The value of b derived hereis nearly identical to the value derived in Beatty et al. (2019). If both values were systematically inaccurate, we would expect a deviation in ${(R_p/R_s)}^2$. Thus, the precise agreement of our optical value with the near-IR value of Beatty et al. (2019) also indicates a good accuracy of the derived b value of KELT-1b.

7 Conclusion

In this work, we present the detection and analysis of the phase curve of KELT-1b at optical wavelengths, analyzing data taken by the Transiting Exoplanet Survey Satellite (TESS) during cycle 2 and sector 17. Due to the high mass of KELT-1b, there are visible distortions in the TESS light curve caused by the gravitational influence that the brown dwarf companion exerts over the star. First, we modeled the data by combining five components, namely: the primary transit, the secondary eclipse, the phase curve that assumes both reflected light and thermal emission, the Doppler beaming, and the ellipsoidal variation models. From our initial fit, we found evidence of variability in the shape of the phase curve from orbit to orbit. Before jumping to the conclusion that we are observing weather patterns in the atmosphere of KELT-1b, we analyzedTESS photometry in more detail and detect what appears to be stellar activity. After adding this component to our model, we re-fit the data and determined the secondary eclipse depth in the TESS bandpass to be 304 ± 75 parts-per-million (ppm). Furthermore, we determined the amplitude of the phase curve to be 129 ± 27 ppm, with a corresponding eastward offset between the region of maximum brightness and the substellar point of 19.2 ± 9.6 degrees; the latter showing good agreement with Spitzer measurements. We determine the day-side brightness temperature in the TESS bandpass of 3201 ± 147 K, slightly higher than what is inferred from Spitzer 3.6 and 4.5 μm data.

Using the error-weighted mean of the Spitzer and TESS brightness temperatures, we computed the Bond albedo of KELT-1b to be A_B = 7.1 ± 7.7%, consistent with zero at 1σ, and a heat redistribution efficiency of ϵ = 0.052 ± 0.006. If the TESS eclipse depth is due to both reflected light and thermal emission, then we derive a geometric albedo of A_g = 25 ± 11%, greater than the geometric albedos of most other hot Jupiters. For this exercise, we are effectively assuming that the day-side brightness temperature in the TESS band is the error weighted mean of the Spitzer and TESS brightness temperatures, so the high geometric albedo inferred might, in fact, be caused by a non-blackbody spectrum, rather than reflected light.

Taking the Bond and geometric albedos at face value, we then asked whether they are self-consistent. We estimated the fraction of incident stellar light in the TESS bandpass: 31%. The minimum Bond albedo should be the product of the geometric albedo and the fraction of incident stellar light in the TESS bandbass, or 7.6%. This is consistent with the Bond albedo inferred from day and night temperatures, but requires that the geometric albedo is nearly zero outside of the TESS bandpass.

The published H and K-band eclipse measurement suggest a significantly higher day-side temperature than either the TESS of Spitzer measurements. We attribute this to unmodeled ellipsoidal variation, which would appear as a deeper eclipse depth in observations with insufficient phase coverage. We therefore opt to ignore eclipse depths based on incomplete phase coverage, and obtain a day-side effective temperature of T_eff = 2946 ± 71 K, with an improvement of ${\chi }_{\text{red}}^2$ from 6.8 to 1.6.

A self-consistent 1D radiative transfer model can explain the day-side fluxes from both TESS and Spitzer with thermal emission alone without invoking reflected light. The somewhat lower brightness temperature in theSpitzer IRAC 2 band compared to that in the TESS band can be explained by CO absorption due to a non-inverted temperature profile.

The Spitzer night-side brightness temperatures place an upper limit of ~2000 K on the internal temperature of KELT-1 b. In the limit of no day-night redistribution, the night side data can be explained using an atmospheric model with an internal temperature of ~1500–2000 K, somewhat higher than the ~850 K internal temperature expected for the brown dwarf (see Beatty et al. 2019) if it were unaffected by the stellar irradiation and if the age of the host star is accurate. However, we note that some of the night-side flux which is being attributed to internal flux could potentially be contributed by day-night energy redistribution, which is not included explicitly in the night side model.

The internal temperature of KELT-1 b may be linked to its inflated radius, considering that the radius for such a brown dwarf would be expected to be below 1 R_J if it were isolated and ≈2 Gyr old (Burrows et al. 1997; Baraffe et al. 2003; Siverd et al. 2012). This inflated radius may be a result of an elevated entropy in its interior and a high T_int. Such a situation is possible if the brown dwarf is young (≈10^8.5 yr) and/or if the irradiation has delayed its cooling significantly; in light of this, it would be worth revisiting the age of the host star, carefully accounting for tidal effects and utilizing the rotational activity-age relationship (Mamajek & Hillenbrand 2008). As such, KELT-1b could provide a valuable testing ground for theories of inflated radii of highly irradiated giant planets.

During the review process of the manuscript of this work, Beatty et al. (2020) published an independent investigation of the same TESS photometry data of KELT-1b. Their derived model parameter values for the secondary eclipse and phase curve variations are broadly consistent with the results presented here.

Acknowledgements

We thank the anonymous referee for the helpful comments, and constructive remarks on this manuscript. C.v.E. and G.T. acknowledge support from the European Social Fund (project No. 09.3.3-LMT-K-712-01-0103) under grant agreement with the Lithuanian Science Council (LMTLT). Funding for the Stellar Astrophysics Centre is provided by The Danish National Research Foundation (Grant agreement No.: DNRF106). This work was supported by a research grant (00028173) from VILLUM FONDEN. G.T. acknowledges Europlanet 2024 RI project funded by the European Union’s Horizon 2020 Research and Innovation Programme (Grant agreement No. 871149). This work made use of PyAstronomy, the SIMBAD data base and VizieR catalog access tool, operated at CDS, Strasbourg, France, and of the NASA Astrophysics Data System (ADS).

Appendix A Posterior distributions

Fig. A.1 Posterior distributions for the primary transit parameters assuming a circular orbit for KELT-1b. In the figure, the red points indicate the best-fit values, and the colors of the contours highlight the 1, 2 and 3σ intervals. To allow for a better visual inspection of the uncertainties, the posteriors were shifted to the best-fit values specified in Table 1.

Fig. A.2 Posterior distributions for the eleven parameters determined from our simultaneous fit, M1.

Fig. A.3 Posterior distributions for the twelve parameters determined from our simultaneous fit, M3.

References

All Tables

All Figures

	Fig. 1 Detrended flux of KELT-1 observed by TESS. SAP photometry is indicated with red points and PDC with black circles. Both light curves are shown after our internal data preparation. The fifteen transits analyzed in this work are indicated in the bottom of the figure with red lines.
In the text

Fig. 2 Third light contamination in the TESS aperture. Top: target Pixel File (TPF) of KELT-1 showing the aperture mask. Three stars are included and completely blurred. Bottom: ~3 × 3 arcmin field of view centered around KELT-1, the brightest star in the field. The mask and the pixels are schematized with white thick and thin lines, respectively. The field of view has been rotated to be aligned with the ecliptic.

In the text

	Fig. 3 Primary transit light curve and residuals for KELT-1b for a circular orbit as a function of the orbital phase. Black points denote TESS data. Continuous red line corresponds to the best-fit primary transit model. The residuals have been artificially shifted below the transits.
In the text

	Fig. 4 Spectral energy distribution (SEDs) of KELT-1. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the bandpass. Blue symbols are the model fluxes from the best-fit Kurucz atmosphere model (black).
In the text

Fig. 5 Periodicity of the residual light curve modulation. Top: phase dispersion minimization evaluated in the frequency range between 0.4 and 0.9 c/d. Bottom: Lomb-Scargle periodogram (amplitude in ppm versus frequency in cd−1) for the residuals of KELT-1, showing only the frequency range where most of the power lies. Overplotted in blue is the fit of a Gaussian function over the power to estimate an error for νmin. The vertical orange line indicates the measured projected rotation velocity of the star as determined by Siverd et al. (2012). For comparison, the frequency of the orbital motion of KELT-1b is 0.821 c/d.

In the text

Fig. 6 Data (black points) versus best-fit model (red continuous line) as a function of the orbital phase of KELT-1b. From top to bottom: full data set; same data, but with the best-fit primary transit, PT(t), removed; ellipsoidal variation, EV(t), removed; Doppler beaming and stellar activity, DB(t) and SA(t), respectively, removed, at which point the secondary eclipse and the phase curve are visible; residuals, once the secondary eclipse and the phase curve are removed. The last two panels are shown at a 0.02 orbital phase binning. The last four panels are given in ppm.

In the text

	Fig. 7 Relative flux of TESS data binned as specified before, showing only the phase curve and secondary eclipse components. In red and blue bars, we indicate the day- and night-side fluxes.
In the text

	Fig. 8 Eclipse depth, in ppm, as a function of wavelength. Literature values are taken from Croll et al. (2015); Beatty et al. (2014, 2017, 2019). Our new data point is plotted in red. The gray area shows a 1σ contour of the equilibrium temperature of KELT-1b, determined from all data. Transmission responses are plotted in black continuous lines, with the exception of TESS, highlighted in red.
In the text

Fig. 9 Self-consistent modeling of the atmosphere of KELT-1b. Top: model day-side spectra (left) and pressure-temperature (P-T) profiles (right), assuming solar metallicity and an internal temperature of 1500 K. The red spectrum and P-T profile correspond to a C/O=0.5 and depleted TiO/VO (e.g., due to gravitational settling), while the blue model corresponds to a C/O = 1. The black dashed lines correspond to blackbody spectra at 3000 and 3200 K. The cloudy, Tint = 1000 K model from Beatty et al. (2020) is shown in grey for reference. Bottom: model night-side spectra (left) and pressure-temperature profiles (right) for KELT-1 b assuming internal temperatures of 1500 K (red lines) and 2000 K (purple lines). Dashed black lines correspond to blackbody spectra at 750, 1000 and 1500 K. Flux is shown in SI units (W sr−1 m−3). In both panels, data are shown by black points and error bars. Coloured circles show the binned model spectra. Vertical gray bands show the TESS and Spitzer bandpasses.

In the text

	Fig. A.1 Posterior distributions for the primary transit parameters assuming a circular orbit for KELT-1b. In the figure, the red points indicate the best-fit values, and the colors of the contours highlight the 1, 2 and 3σ intervals. To allow for a better visual inspection of the uncertainties, the posteriors were shifted to the best-fit values specified in Table 1.
In the text

	Fig. A.2 Posterior distributions for the eleven parameters determined from our simultaneous fit, M1.
In the text

	Fig. A.3 Posterior distributions for the twelve parameters determined from our simultaneous fit, M3.
In the text