TOI-2109: An Ultrahot Gas Giant on a 16 hr Orbit

TOI-2109, listed in the TESS input catalog (TIC; Stassun et al. 2018) as TIC 392476080, was observed by the TESS spacecraft from UT 2020 May 13 to 2020 June 8 during the Sector 25 campaign. Photometry of the target was obtained from the full frame images (FFIs), which have a cadence of 30 minutes. Initial light-curve extraction and analysis were carried out using the MIT Quick Look Pipeline (QLP; Huang et al. 2020a, 2020b). Subsequent vetting identified a transit-like signal in the target's light curve with a depth of roughly 7000 ppm and a duration of ∼1.8 hr that occurred every ∼0.67 days. The planet candidate was added to the list of TESS objects of interest (TOIs) as TOI-2109.01. Table 1 lists astrometric and photometric information about the target provided in TIC Version 8.1 (Stassun et al. 2019).

Table 1. Target Information

Parameter	Value	Source
TIC	392476080	TIC V8 a
R.A.	16h52m45s	Gaia DR2 b
Decl.	$+16^\circ 34^{\prime} 48^{\prime\prime}$	Gaia DR2 b
μra (mas yr−1)	−8.449 ± 0.043	Gaia DR2 b
μdec (mas yr−1)	−9.257 ± 0.051	Gaia DR2 b
Parallax (mas)	3.788 ± 0.039	Gaia DR2 b
Distance (pc)	262.04 ± 2.73	Gaia DR2 b
Epoch	2015.5	Gaia DR2 b
BT (mag)	10.731 ± 0.032	Tycho-2 c
VT (mag)	10.268 ± 0.029	Tycho-2 c
GBP (mag)	10.3638 ± 0.0011	Gaia DR2 b
G (mag)	10.11376 ± 0.00034	Gaia DR2 b
GRP (mag)	9.73916 ± 0.00094	Gaia DR2 b
TESS (mag)	9.7857 ± 0.0061	TIC V8 a
J (mag)	9.382 ± 0.024	2MASS d
H (mag)	9.129 ± 0.026	2MASS d
K (mag)	9.070 ± 0.021	2MASS d
W1 (mag)	9.059 ± 0.023	WISE e
W2 (mag)	9.093 ± 0.020	WISE e
W3 (mag)	9.062 ± 0.030	WISE e

Notes.

^a Stassun et al. (2018). ^b Gaia Collaboration et al. (2018). ^c Høg et al. (2000). ^d Cutri et al. (2003). ^e Cutri et al. (2013).

Download table as:  ASCIITypeset image

Two different extractions of the TOI-2109 photometry are available on the Mikulski Archive for Space Telescopes ⁴⁹ (MAST). The first extraction was carried out using the Science Processing and Operations Center (SPOC) pipeline based at NASA Ames Research Center (Jenkins et al. 2016) as part of the TESS Light Curves from Full Frame Images (TESS-SPOC) High Level Science Products project (Caldwell et al. 2020). The corresponding datafile contains two versions of the light curve: (1) the simple aperture photometry (SAP) light curve, which consists of the flux extracted from the optimized aperture on the FFI, and (2) the pre-search data conditioning (PDC) light curve, which was produced by the SPOC pipeline's detrending routine that utilizes the empirically determined co-trending basis vectors (CBVs) to remove instrumental systematics trends common to all sources on the detector (Smith et al. 2012; Stumpe et al. 2012, 2014). The second data extraction available on MAST is the QLP-derived photometry; the SAP light curve contained therein was extracted using a different aperture than the SPOC light curves.

We used the ExoTEP pipeline (Benneke et al. 2019; Wong et al. 2020a) to analyze the TESS photometry. Prior to fitting the light curves, we removed all points in the time series with a NaN flux value or a nonzero quality flag (as indicated by the SPOC pipeline). We also applied a 16 point wide moving median filter to the transit-trimmed light curves and removed 5σ outliers. From an initial time series of 1231 points, these preprocessing steps removed 82 points (6.7% of the data). Next, we divided the light curve into four segments that are separated by the scheduled momentum dumps (i.e., when the spacecraft thrusters were engaged to reset the onboard reaction wheels, leading to increased pointing jitter and poorer data quality) and data downlink interruptions, similar to what was done in previous analyses of TESS photometry (e.g., Wong et al. 2020b, 2020d, 2020e, 2021). Figure 1 shows the three versions of the TOI-2109 light curve: SPOC-PDC, SPOC-SAP, and QLP-SAP.

Zoom In Zoom Out Reset image size

Aside from the transits, there are periodic brightness modulations that are commensurate with the orbital period (i.e., a phase-curve signal), which were initially discerned from careful inspection of the phase-folded QLP photometry produced as part of the initial candidate vetting process. In addition, there are long-term flux trends in the data attributable to uncorrected instrumental systematics. The time-correlated noise is particularly severe in the SPOC-PDC light curve, which also displays sharp flux ramps near the beginning and end of several data segments.

While the SPOC detrending routine that produces the PDC photometry typically reduces instrumental systematics and improves time-correlated noise behavior, the presence of stellar variability and/or residual systematics features that are not shared by other sources on the detector can result in poorer data quality, as was previously reported, for example, in the TESS light curve of the active planet-hosting star WASP-19 (Wong et al. 2020b). In our TESS light-curve analysis, we only considered the SPOC-SAP and QLP-SAP light curves.

We acquired ground-based time-series photometry of primary transits of TOI-2109 as part of the TESS Follow-up Observing Program ⁵⁰ (TFOP) Sub Group 1 (SG1; seeing-limited time-series photometry) collaboration. The TFOP SG1 network includes both professional and amateur astronomers at more than a hundred observatories around the world. Observers choose targets to follow up using the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen 2013).

For TOI-2109, 20 full- and partial-transit observations were obtained between 2020 July 31 and 2021 June 26. The photometric data sets contributed by TFOP SG1 observers were uploaded to the ExoFOP-TESS repository. ⁵¹ These observations are summarized in Table 2; detailed descriptions of the instruments and observing methodology are provided in the following subsections.

Table 2. TFOP Primary Transit Observations

Date (UT)	Telescope	Filter	Coverage	Exposure time (s)	Duration (hr)	Aperture radius (arcsec)
2020 Jul 31	FLWO KeplerCam (1.2 m)	$z^{\prime}$	Partial	5	2.8	7.4
2020 Jul 31	ULMT (0.6 m)	$r^{\prime}$	Full	16	2.9	4.3
2020 Jul 31 a	LCOGT McD (0.4 m)	$g^{\prime} ,i^{\prime}$	Full	17, 25	3.8	6.8
2020 Aug 2	MLO (0.36 m)	I	Full	30	3.9	6.7
2020 Aug 4 a	TCS MuSCAT2 (1.5 m)	$g^{\prime} ,r^{\prime} ,i^{\prime} ,{z}_{s}$	Full	3, 3, 3, 3	2.5	13.9
2020 Aug 24	LCOGT SSO (1.0 m)	B, zs	Full	10, 24	3.6	15.6
2021 Apr 7	WBRO (0.24 m)	R	Full	100	4.8	8.7
2021 Apr 7	LCOGT SAAO (1.0 m)	$i^{\prime}$	Full	16	5.1	9.0
2021 Apr 7	GdP (0.4 m)	$i^{\prime}$	Full	20	4.4	6.6
2021 Apr 8	LCOGT SSO (1.0 m)	$i^{\prime}$	Partial	16	3.5	7.4
2021 Apr 9	LCOGT CTIO (1.0 m)	$i^{\prime}$	Partial	16	3.5	7.4
2021 Apr 9	ULMT (0.6 m)	$i^{\prime}$	Full	32	4.4	4.3
2021 Apr 9	FLWO KeplerCam (1.2 m)	$z^{\prime}$	Full	5	5.2	6.7
2021 Apr 11	FLWO KeplerCam (1.2 m)	$z^{\prime}$	Full	5	4.7	6.7
2021 May 14 a	TCS MuSCAT2 (1.5 m)	$g^{\prime} ,r^{\prime} ,i^{\prime} ,{z}_{s}$	Full	25, 25, 50, 50	3.4	10.9
2021 May 22 a	TCS MuSCAT2 (1.5 m)	$g^{\prime} ,r^{\prime} ,i^{\prime} ,{z}_{s}$	Full	10, 5, 10, 5	3.9	10.9
2021 May 24	FTN MuSCAT3 (2.0 m)	$g^{\prime} ,r^{\prime} ,i^{\prime} ,{z}_{s}$	Full	20, 12, 15, 33	5.7	6.1
2021 May 25 a	TCS MuSCAT2 (1.5 m)	$g^{\prime} ,r^{\prime} ,i^{\prime} ,{z}_{s}$	Full	5, 5, 10, 5	3.5	10.9
2021 Jun 12	TCS MuSCAT2 (1.5 m)	$g^{\prime} ,r^{\prime} ,i^{\prime} ,{z}_{s}$	Full	5, 5, 10, 5	3.0	10.9
2021 Jun 26	FTN MuSCAT3 (2.0 m)	$g^{\prime} ,r^{\prime} ,i^{\prime} ,{z}_{s}$	Full	20, 12, 15, 33	6.6	5.3

Note.

^a These four MuSCAT2 observations and the McD light curve were affected by dome issues and clouds, respectively, and were not included in the final set of ground-based light curves analyzed in this paper.

Download table as:  ASCIITypeset image

Unless otherwise noted, the data reduction and photometric extraction of the TFOP SG1 observations were performed using the AstroImageJ (AIJ) software package (Collins et al. 2017). The extraction aperture radii ranged from 43 to 156 across the data sets. The nearest Gaia DR2 star to TOI-2109 is a faint neighbor at a projected separation of 18'', well outside all of the apertures used for the ground-based photometry. In the case of the heavily defocused observations on 2020 August 4 and 24, some light from the neighboring star would have fallen within the aperture; however, given that the neighbor is 7 mag fainter than the target star, the dilution is negligible. These ground-based time-series observations excluded all stars in the vicinity of TOI-2109 (at separations larger than ∼1'') as the source of the transit signal. We describe our high-angular-resolution search for smaller-separation companion stars in Section 2.4.

2.2.1. FLWO KeplerCam

2.2.2. ULMT

2.2.3. LCOGT McD, SSO, SAAO, and CTIO

2.2.4. MLO

2.2.5. TCS MuSCAT2

2.2.6. WBRO

2.2.7. GdP

2.2.8. FTN MuSCAT3

On UT 2021 March 6, we observed a secondary eclipse of TOI-2109b in the K_s band (λ_eff = 2.13 μm) with the Wide-field Infrared Camera (WIRC; Wilson et al. 2003) on the Hale 200'' Telescope at Palomar Observatory, California, USA. The data were taken with a beam-shaping diffuser that increased our observing efficiency and improved guiding stability on this bright target (Stefansson et al. 2017; Vissapragada et al. 2020). In order to mitigate a known detector systematic at short exposure times, we initiated the observations with a four-point dither near the target to construct a background frame. We then collected 239 exposures, each consisting of 15 coadds of 0.92 s for a total per-image exposure time of 13.8 s.

The data were dark-subtracted, flat-fielded, and corrected for bad pixels using the same methodology as in Vissapragada et al. (2020). We scaled and subtracted the aforementioned background frame from each image to remove the full frame background structure. Next, we performed circular aperture photometry with the photutils package (Bradley et al. 2020) on TOI-2109 and two nearby comparison stars of similar brightness located within the FOV. We experimented with a range of photometric extraction apertures from 10 to 25 pixels in 1 pixel steps (with 025 pixel⁻¹), eventually selecting a 13 pixel (325) aperture that minimized the per-point scatter of the fitted photometry.

The raw extracted photometry is shown in the top panel of Figure 2. The light curve begins just before the start of ingress , when the target reached an airmass above 2. While the sky was clear throughout the ∼3.8 hr observation, the sky background in the vicinity of TOI-2109 varied considerably, particularly during the last hour, resulting in more severe residual systematics.

Zoom In Zoom Out Reset image size

A bound or line-of-sight companion in close proximity to the target can create a false positive transit signal if it is an eclipsing binary. The so-called "third-light" flux from a companion star can also yield an underestimated planetary radius if not accounted for in the transit model (Ciardi et al. 2015). Likewise, the photometric contamination can lead to nondetections of small planets residing within the same exoplanetary system (Lester et al. 2021). The discovery of binary systems containing close, bound companions, which exist around nearly half of FGK-type stars (Matson et al. 2018), can help further our understanding of exoplanet formation, dynamics, and evolution (Howell et al. 2021). In order to search for close-separation companions unresolved by TESS and other seeing-limited ground-based follow-up observations, we carried out high-resolution speckle imaging of TOI-2109.

The target was observed on UT 2020 September 17 using the 'Alopeke speckle instrument on Gemini-North. 'Alopeke provides simultaneous speckle imaging in two bands (562 and 832 nm), producing a reconstructed image with robust contrast limits on companion detections (e.g., Howell et al. 2016). Five sets of 1000 × 0.06 s exposures were collected and passed through a standard Fourier analysis in our reduction pipeline (see Howell et al. 2011). Figure 3 shows the resultant 5σ contrast curves and the reconstructed speckle images in both bands. We find TOI-2109 to be a single star with no companion brighter than 5–9 mag below the target star's brightness from the diffraction limit (∼20 mas) out to 1.2''. At the distance of TOI-2109 (d = 262 pc; Table 1), these angular limits correspond to physical separations of 5 au to 314 au, respectively.

Zoom In Zoom Out Reset image size

2.5.1. TRES

2.5.2. FIES

The TESS light curves of TOI-2109 (Figure 1) show photometric variability that is synchronous with the orbital period of the planet. We used a full-orbit phase-curve model to fit the light curves. To address the instrumental systematics present in the photometry, we experimented with several different methods for detrending the light curves. Significant periodic brightness modulations attributed to stellar variability were also detected in the TESS photometry and included in our light-curve model.

4.1.1. Full-orbit Phase-curve Model

Following previous phase-curve analyses of TESS data (e.g., Shporer et al. 2019; Wong et al. 2020b, 2020d, 2020e, 2021), we used a simple sinusoidal phase-curve model that treats the stellar and planetary fluxes separately:

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

$$\begin{eqnarray}&&{F}_{p}(t)=\bar{{f}_{p}}-{A}_{\mathrm{atm}}\cos (\phi +\psi ),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{F}_{* }(t)=1-{A}_{\mathrm{ellip}}\cos (2\phi )+{A}_{\mathrm{Dopp}}\sin (\phi ).\end{eqnarray} \tag{ 3 }$$

The orbital phase ϕ is given by 2π(t − T_c ), where T_c is the mid-transit time. The transit and eclipse light curves are represented by λ_t and λ_e . The planetary flux F_p , which is expected to be dominated by thermal emission from the atmosphere, is defined relative to the average flux level $\bar{{f}_{p}}$, with the semiamplitude of the atmospheric brightness modulation represented by A_atm. The parameter ψ denotes a shift in the planetary phase curve, which may arise from an offset dayside hotspot due to superrotating equatorial winds or inhomogeneous clouds. In this parameterization, the secondary eclipse depth is ${D}_{d}=\bar{{f}_{p}}-{A}_{\mathrm{atm}}\cos (\pi +\psi )$, and the nightside flux is${D}_{n}=\bar{{f}_{p}}-{A}_{\mathrm{atm}}\cos \psi$.

The stellar flux F_* contains terms corresponding to two separate physical processes (e.g., Faigler & Mazeh 2011, 2015; Shporer 2017): (1) the tidal response of the stellar surface to the gravitational pull of the orbiting companion, typically referred to as ellipsoidal distortion (A_ellip), and (2) the modulation in the band-integrated flux due primarily to the periodic Doppler shifting of the stellar spectrum, i.e., Doppler-boosting (A_Dopp). Here the sign convention for the various phase-curve amplitudes was chosen so as to produce positive values under normal circumstances.

In an unconstrained fit, A_Dopp and ψ are degenerate, as a phase shift in the planetary phase-curve signal can be absorbed by the coefficient of the $\sin \phi$ term (i.e., where the Doppler-boosting signal lies). Therefore, when fitting the TESS light curves alone, we did not consider any phase offset and simply fit for the total harmonic power at the sine of the orbital frequency, which we denote as A₁. It follows that the values for A₁ that we obtained from our dedicated TESS phase-curve analysis contain contributions from both Doppler boosting on the host star and a phase offset in the planet's atmospheric brightness modulation. In our final global analysis of all available data sets (Section 4.4), we leveraged the constraint on planet mass provided by the RV measurements to self-consistently model the Doppler-boosting signal and disentangle the phase offset.

We mention in passing that the ellipsoidal distortion component of the star's phase-curve modulation contains additional higher-order terms at other harmonics of the cosine (e.g., Morris 1985; Morris & Naftilan 1993). The model in Equation (3) only includes the first-order term at the first harmonic of the cosine. However, the second-order term, which is at the second harmonic (i.e., $\cos 3\phi$), has a theoretically predicted amplitude that is at least an order of magnitude smaller than the leading-order term. In the context of our light-curve fits, this makes the expected second-order ellipsoidal distortion amplitude smaller than the characteristic uncertainties on the phase-curve amplitudes. We also note that previous analyses of Kepler phase curves revealed anomalously large second-harmonic phase-curve signals on HAT-P-7 and KOI-13, which were attributed to the large spin–orbit misalignments in both of those systems (e.g., Esteves et al. 2013, 2015). In contrast, the TOI-2109 system is well aligned (see Section 4.4.1), and we therefore do not expect an additional contribution to the photometric modulation at the second harmonic. Indeed, when fitting the TESS light curve using only the leading-order term (i.e., $\cos 2\phi$), we did not find any periodicity in the residuals at the second harmonic of the orbital phase (see Figure 4). Therefore, we ignored all higher-order terms of the ellipsoidal distortion in the final analysis.

Zoom In Zoom Out Reset image size

In the ExoTEP pipeline, both the transit and secondary eclipse light curves are modeled using batman (Kreidberg 2015). For all of our TESS light-curve fits, we allowed the mid-transit time T_c , orbital period P, radius ratio R_p /R_*, impact parameter b, and scaled orbital semimajor axis a/R_* to vary freely. Due to the low 30 minute cadence of the TESS data, the ingress and egress are not well resolved in the light curve. As such, we did not allow the quadratic limb-darkening coefficients u₁ and u₂ to vary freely, but instead placed Gaussian priors. The mean values were set to the coefficients tabulated in Claret (2018) for the nearest available combination of stellar parameters (u₁ = 0.33, u₂ = 0.22), and the width of each Gaussian was conservatively set at 0.05.

4.1.2. Systematics Detrending

The SPOC-SAP light curve is not corrected for instrumental systematics. We downloaded the CBVs ⁵² for the specific TESS camera and detector on which the target was located (Camera 1, CCD 4) and carried out a customized detrending procedure. The systematics were modeled as a linear combination of the CBVs ν_k (t):

$$\begin{eqnarray}&&{S}_{\mathrm{CBV}}(t)={c}_{0}+\sum _{k=1}^{8}{c}_{k}{\nu }_{k}(t).\end{eqnarray} \tag{ 4 }$$

Eight CBVs were determined by the SPOC pipeline in Sector 25 and included in the downloaded light-curve files.

We fit the SPOC-SAP light curve to the combined phase-curve and systematics model using two approaches. For the first approach, we included all eight CBVs in the detrending model (CBV-full), while for the second approach we only included the CBVs with significant coefficients in the fit (CBV-opt). Each of the four data segments was fit separately, with the optimal combination of CBVs determined for each segment using the Bayesian information criterion (BIC). The sets of CBVs included in the CBV-opt fits are (3, 7), (2, 3, 7), (2, 3, 6, 7), and (1, 2, 3, 6) for segments 1–4, respectively.

The QLP-SAP light curve was generated using a different aperture than the SPOC pipeline. Therefore, the SPOC-generated CBVs are not applicable, and we instead utilized a standard polynomial in time to model the long-term systematics in each data segment:

$$\begin{eqnarray}&&{S}_{{\rm{poly}}}(t)=\displaystyle \sum _{k=0}^{N}{a}_{k}(t-{t}_{0}{)}^{k}.\end{eqnarray} \tag{ 5 }$$

Here, t₀ is the first timestamp of the segment, and N is the order of the detrending polynomial. When determining the optimal polynomial order for each data segment, we considered both the BIC and the Akaike information criterion (AIC; Akaike 1974); the AIC penalizes the addition of free parameters less severely than the BIC, resulting in higher-order polynomials. With simultaneous stellar variability modeling (Section 4.1.3), the optimal orders for the four segments when considering the BIC are 1, 3, 2, and 2, while the AIC prefers 10, 8, 10, and 5, respectively.

For both SPOC-SAP and QLP-SAP light-curve fits, the best-fit systematics model for each segment was removed from the photometry prior to the joint fits of all four segments. The joint fits did not include any additional systematics modeling.

4.1.3. Stellar Variability

Close inspection of the residuals from the SPOC-SAP light-curve fits revealed additional short-term time-correlated brightness variations. Figure 4 shows the Lomb–Scargle periodogram of the residuals as a function of the planet's orbital frequency. While the phase-curve model has fit away all photometric modulation that is synchronous with the orbit, there are two prominent peaks corresponding to periods of roughly Π₁ = 0.61 days and Π₂ = 0.97 days; additionally, a smaller peak is located at twice the frequency of the 0.61 day period, i.e., a variation at the first harmonic. These signals indicate that two distinct stellar variability frequencies are present in the data. The 0.97 day signal lies close to the host star's rotation period as implied by the measured $v\sin {i}_{* }$ and R_* (${P}_{\mathrm{rot}}/\sin {i}_{* }=1.05\pm 0.04$ days; see Section 3).

We checked the Lomb–Scargle periodograms for the QLP-SAP light curves of all targets within $4^{\prime}$ of TOI-2109: TIC 284467967 (T = 13.2 mag), TIC 284467994 (T = 13.5 mag), TIC 392476048 (T = 11.6 mag), and TIC 392476087 (T = 11.6 mag). None of them shows any periodicity near Π₁ or Π₂. Therefore, we assumed that the stellar variability signal belongs to the target host star.

We modeled the two variability signals using generalized sinusoids at the characteristic periods Π₁ and Π₂:

$$\begin{eqnarray}&&\begin{array}{l}{{\rm{\Psi }}}_{1}(t)=1+{\alpha }_{1}\sin ({\phi }_{1})+{\beta }_{1}\cos ({\phi }_{1})\\ \quad \,+\alpha {{\prime} }_{1}\sin (2{\phi }_{1})+\beta {{\prime} }_{1}\cos (2{\phi }_{1}),\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{2}(t)=1+{\alpha }_{2}\sin ({\phi }_{2})+{\beta }_{2}\cos ({\phi }_{2}).\end{eqnarray} \tag{ 7 }$$

The zero-points of the corresponding variability phases are set at the reference time T_ref = 2,458,997 BJD_TDB, which is the integer Julian date closest to the median of the TESS time series: ϕ₁ = 2π(t − T_ref)/Π₁ and ϕ₂ = 2π(t − T_ref)/Π₂. The coefficients marked with primes are the amplitudes of the first harmonic terms at the ∼0.61 day period.

To simultaneously retrieve the parameter values from the full-orbit phase-curve and stellar variability models, we multiplied Ψ₁(t) and Ψ₂(t) by F_*(t) in Equation (1). The amplitudes and variability periods were allowed to vary freely in the fits. For completeness, we also present the results from the initial SPOC-SAP light-curve fits that did not account for stellar variability. The stellar variability was included in the full light-curve model when fitting the QLP-SAP photometry; QLP-SAP light-curve fits without stellar variability modeling necessitated very high orders (>15) in the systematics detrending polynomials for every data segment, making the analysis untenable.

4.1.4. Fitting and Model Selection

ExoTEP uses the affine-invariant Markov Chain Monte Carlo (MCMC) routine emcee (Foreman-Mackey et al. 2013) to compute the posterior distributions of all fit parameters. Each fit consisted of two steps. In the first MCMC run, we included a per-point uncertainty parameter σ, which was allowed to vary freely to ensure that the resultant best-fit model has a reduced χ² of one. The value of σ represents the scatter in the light curve and includes the contributions from both photon noise (i.e., white noise) and time-correlated noise (i.e., red noise) at the 30 minute cadence of the observations.

To address the effect of red noise at longer timescales on the TESS light-curve fit results, we computed the standard deviation of the residuals, binned at various intervals, and compared the resultant values to the expected $1/\sqrt{n}$ scaling for pure white noise (e.g., Pont et al. 2006; see Wong et al. 2020d for details on the specific implementation described here). In all cases, the binned residual scatter showed a positive deviation from the white noise trend (see Figure 5 for an example plot). This indicates the presence of significant time-correlated noise at timescales longer than a few hours, which is particularly relevant for our phase-curve analysis, because the characteristic flux modulations from the atmospheric brightness variation and ellipsoidal distortion occur on those timescales.

Zoom In Zoom Out Reset image size

We accounted for the contribution of additional red noise in our fits by computing the average fractional deviation β between the binned residual scatter and the white noise trend for bin sizes up to 16 (i.e., 8 hr, or roughly half the orbital period) and multiplying this ratio by the fitted per-point uncertainty σ from the first MCMC run. Typical values of β for fits that included stellar variability modeling ranged from 1.15 to 1.25. We then fixed the flux uncertainty to this inflated value σ_r ≡ β σ and ran the MCMC analysis a second time. The resultant posteriors are broader, reflecting the added contribution of red noise to the overall photometric uncertainty.

We carried out fits to the full TESS light curves from both SPOC and QLP using various systematics detrending techniques (SPOC: CBV-full and CBV-opt; QLP: Poly-BIC and Poly-AIC). Table 5 lists the results from these fits, along with the corresponding fitted scatter levels σ and red-noise-inflated per-point uncertainties σ_r . Unsurprisingly, the SPOC-SAP fits that did not account for stellar variability yielded significantly higher scatter at both the native 30 minute cadence and longer timescales. Meanwhile, the noise level in the residuals from the QLP-SAP fits is systematically higher than the residual scatter from the SPOC-SAP fits.

Table 5. Comparison of SPOC and QLP Light-curve Fits with Different Fitting and Detrending Methods

	SPOC-SAP	SPOC-SAP	SPOC-SAP	SPOC-SAP	QLP-SAP	QLP-SAP
	with Variability	with Variability	without Variability	without Variability	with Variability	with Variability
Parameter	CBV-full	CBV-opt	CBV-full	CBV-opt	Poly-BIC	Poly-AIC
Transit and Orbital Parameters
Rp /R*	0.0789 ± 0.0020	0.0790 ± 0.0027	${0.0765}_{-0.0012}^{+0.0030}$	${0.0768}_{-0.0012}^{+0.0025}$	${0.0789}_{-0.0018}^{+0.0028}$	0.0795 ± 0.0024
Tc a	997.16653 ± 0.00016	997.16652 ± 0.00017	997.16653 ± 0.00024	997.16650 ± 0.00024	997.16636 ± 0.00019	997.16642 ± 0.00017
P (days)	0.672469 ± 0.000015	0.672469 ± 0.000015	0.672470 ± 0.000022	0.672468 ± 0.000020	0.672478 ± 0.000018	0.672483 ± 0.000014
b	$0.{61}_{-0.18}^{+0.10}$	${0.62}_{-0.24}^{+0.13}$	${0.29}_{-0.24}^{+0.35}$	${0.37}_{-0.25}^{+0.24}$	${0.53}_{-0.20}^{+0.17}$	${0.59}_{-0.23}^{+0.13}$
a/R*	$2.{60}_{-0.24}^{+0.30}$	${2.57}_{-0.30}^{+0.39}$	${3.04}_{-0.51}^{+0.12}$	${2.96}_{-0.38}^{+0.19}$	${2.73}_{-0.37}^{+0.27}$	${2.62}_{-0.31}^{+0.34}$
u1 b	0.34 ± 0.04	0.34 ± 0.05	0.34 ± 0.05	0.34 ± 0.05	0.34 ± 0.05	0.33 ± 0.05
u2 b	0.24 ± 0.05	0.22 ± 0.04	0.23 ± 0.04	0.23 ± 0.05	0.23 ± 0.05	0.23 ± 0.05
Stellar Variability Parameters
Π1 (day)	0.61393 ± 0.00051	0.61389 ± 0.00048	...	...	0.61450 ± 0.00058	0.61437 ± 0.00051
α1 (ppm)	−11 ± 15	−15 ± 15	...	...	−7 ± 18	−6 ± 17
β1 (ppm)	−149 ± 14	−158 ± 15	...	...	−180 ± 21	−169 ± 17
$\alpha {{\prime} }_{1}$ (ppm)	−71 ± 15	−65 ± 17	...	...	−62 ± 22	−65 ± 18
$\beta {{\prime} }_{1}$ (ppm)	64 ± 15	61 ± 17	...	...	88 ± 18	95 ± 18
Π2 (day)	0.9675 ± 0.0013	0.9676 ± 0.0014	...	...	0.9690 ± 0.0013	0.9688 ± 0.0013
α2 (ppm)	25 ± 17	21 ± 15	...	...	44 ± 24	29 ± 18
β2 (ppm)	−222 ± 14	−213 ± 15	...	...	−262 ± 20	−248 ± 16
Phase-curve Parameters
$\bar{{f}_{p}}$ (ppm)	367 ± 39	363 ± 42	387 ± 61	371 ± 69	384 ± 52	380 ± 44
Aatm (ppm)	360 ± 22	357 ± 21	366 ± 27	367 ± 31	401 ± 29	396 ± 25
Aellip (ppm)	240 ± 19	243 ± 23	223 ± 27	231 ± 28	245 ± 25	249 ± 23
A1 (ppm) c	31 ± 16	40 ± 14	25 ± 24	20 ± 23	25 ± 19	31 ± 16
Derived Parameters
Dd (ppm) d	726 ± 46	720 ± 47	751 ± 71	736 ± 81	785 ± 56	778 ± 47
Dn (ppm) d	8 ± 44	3 ± 46	23 ± 63	9 ± 70	−18 ± 61	−16 ± 54
i (deg)	${\bf{76}}.{{\bf{5}}}_{-4.0}^{+5.1}$	${76.2}_{-5.2}^{+6.6}$	${84.6}_{-9.1}^{+4.5}$	${82.8}_{-6.7}^{+5.0}$	${78.8}_{-6.1}^{+4.9}$	${77.0}_{-5.2}^{+6.1}$
Fit-quality Metrics
σ (ppm) e	295	293	349	351	367	356
σr (ppm) e	359	363	547	549	472	409

Notes.

^a BJD_TDB − 2,458,000. ^b Limb-darkening coefficients were constrained by priors: u₁ = 0.33 ± 0.05, u₂ = 0.22 ± 0.05. ^c A₁ is the total harmonic power in the photometry at the sine of the orbital phase, which includes the Doppler-boosting signal from the host star and any phase shift in the planet's atmospheric brightness modulation. ^d Dayside flux (secondary eclipse depth) and nightside brightness of the planet, derived from the fitted average planetary flux $\bar{{f}_{p}}$ and atmospheric brightness modulation amplitude A_atm. ^e σ: scatter in the residuals from the best-fit phase-curve model; σ_r : per-point uncertainty, inflated to account for red noise.

Download table as:  ASCIITypeset image

When comparing the values from the six listed sets of results, we report a high level of mutual consistency. Most notably, the results from the SPOC-SAP fits that did or did not account for stellar variability agree with each other at much better than the 1σ level, which indicates that our treatment of stellar variability does not have any significant effect on the phase-curve results, aside from improving the time-correlated noise and reducing parameter uncertainties. Likewise, comparisons of fits that utilized full CBV detrending versus optimized CBV detrending show full statistical consistency, as do the QLP-SAP light-curve fits with Poly-BIC versus Poly-AIC detrending.

Looking at the results of the SPOC-SAP and QLP-SAP light-curve fits (with stellar variability modeling) side by side, we find that the transit-shape parameters and orbital ephemeris are mutually consistent to well within 1σ. Similarly, most of the phase-curve parameters and derived quantities such as secondary eclipse depth D_d and nightside flux D_n agree with one another at better than the 1σ level. The exception is the atmospheric brightness modulation amplitude A_atm, for which the QLP-SAP photometry prefers a value that is up to 1.2σ larger than the corresponding measurements derived from the SPOC-SAP light curve. Meanwhile, the stellar variability parameters from the SPOC-SAP and QLP-SAP fits are broadly consistent, with no deviations larger than 2σ. All in all, we find that the astrophysical parameter values of interest are highly robust to the specific choice of photometric extraction, systematics detrending methodology, and stellar variability modeling.

The versions of the SPOC-SAP light curve corrected using the full and optimized CBV detrending methods yielded very similar fit quality with respect to both residual scatter and time-correlated noise. We selected the former for the main results of this paper, on account of the marginally better red noise level.

4.1.5. Results

The phase-folded SPOC-SAP light curve, corrected for systematics using the full CBV detrending model and with the best-fit stellar variability signals removed, is shown in Figure 6. In addition to the high-S/N secondary eclipse with a depth of 726 ± 46 ppm, we retrieved significant (>10σ) phase-curve amplitudes corresponding to the atmospheric brightness modulation of the planet and the ellipsoidal distortion of the star. The nightside flux is consistent with zero. We also obtained a marginal phase-curve amplitude at the sine of the orbital frequency (31 ± 16 ppm).

Zoom In Zoom Out Reset image size

The transit-shape parameters b and a/R_* are not well constrained by the TESS light curve alone, due to the low cadence of the observations. The best-fit values indicate an orbit that is moderately inclined from edge-on. The precision of these values is substantially improved when including the ground-based light curves and spectroscopic transit observations in the final joint fit (Section 4.4). The two fitted stellar variability periods are Π₁ = 0.61393 ± 0.00051 days and Π₂ = 0.9675 ± 0.0013 days. Figure 4 shows the TESS light curve phase-folded on the two different variability periods, with the best-fit full-orbit phase-curve model removed.

The extremely small planet–star separation makes any significant orbital eccentricity highly unlikely. Nevertheless, to probe for possible deviations from a circular orbit, we carried out a fit that included $e\cos \omega$ and $e\sin \omega$ as free parameters. The constraints on eccentricity here are primarily driven by the relative timing of the secondary eclipse. We obtained $e\cos \omega =-0.0005\pm 0.0029$ and $e\sin \omega =0.003\pm 0.045$, corresponding to a formal 2σ upper limit of e < 0.07. Therefore, we conclude that the orbit of TOI-2109b is indeed consistent with circular.

The raw light curves obtained from the various ground-based observations described in Section 2.2 were affected by instrumental systematics and observing conditions, such as airmass and sky background level. To detrend these systematics, we considered all possible linear combinations of relevant quantities, including the measured centroid position of the target (x, y), width of the target's point-spread function Δ_PSF, airmass AM, sky background level in the vicinity of the target F_sky, and the total flux of nearby companion stars F_comp used to derive the differential photometry. For every data set, an additional baseline offset c₀ was used to properly normalize the light curve. We also experimented with modeling a linear trend in time t. In the case of the WIRC K_s -band light curve of the secondary eclipse, the extracted fluxes from the two selected companion stars F₁ and F₂ were included as additional detrending vectors (see Section 2.3). The optimized combinations of detrending vectors were determined through minimizing the BIC for each data set.

In our fits, we only considered ground-based light curves with full transit coverage, due to the possibility of significant biases to the transit timing and transit-shape parameters when modeling partial light curves. Each transit data set was fit using ExoTEP. The B-, $r^{\prime}$-, R-, $i^{\prime}$-, I-, and $z^{\prime}$/z_s -band quadratic limb-darkening coefficients were constrained by priors derived from the tabulated values in Claret et al. (2013), with 1σ Gaussian widths uniformly set to 0.05. The full set of limb-darkening priors used in our analysis is given in Table 6.

Table 6. Priors on Limb-darkening Coefficients

Parameter				Value
u1,TESS				0.33 ± 0.05
u2,TESS				0.22 ± 0.05
u1,B				0.56 ± 0.05
u2,B				0.25 ± 0.05
${u}_{1,g^{\prime} }$				0.54 ± 0.05
${u}_{2,g^{\prime} }$				0.22 ± 0.05
${u}_{1,r^{\prime} }$				0.40 ± 0.05
${u}_{2,r^{\prime} }$				0.23 ± 0.05
u1,R				0.39 ± 0.05
u2,R				0.22 ± 0.05
${u}_{1,i^{\prime} }$				0.33 ± 0.05
${u}_{2,i^{\prime} }$				0.21 ± 0.05
u1,I				0.31 ± 0.05
u2,I				0.20 ± 0.05
${u}_{1,z^{\prime} /{z}_{s}}$				0.26 ± 0.05
${u}_{2,z^{\prime} /{z}_{s}}$				0.21 ± 0.05

Download table as:  ASCIITypeset image

The transit and secondary eclipse depths measured from a light-curve fit can be systematically affected by the unmodeled photometric variability of the system during the eclipsing event. This variability includes the planetary phase curve and modulations in the stellar brightness, which shift the out-of-eclipse baseline. To minimize the possibility of biases in our ground-based transit light-curve fits while simultaneously preserving a sufficient out-of-transit baseline, we excluded all data points that lie more than 0.1 × P from the mid-transit time.

To leverage the high precision of the simultaneous multiband MuSCAT2 and MuSCAT3 photometry, we initiated our ground-based light-curve analysis by jointly fitting each set of transit light curves. The MuSCAT2 transit fit yielded tight constraints on the mid-transit time (T_c = 2,459,378.45955 ± 0.00017 BJD_TDB) and transit-shape parameters (b = 0.743 ± 0.020, a/R_* = 2.274 ±0.057). The first set of MuSCAT3 light curves provided even more precise measurements: T_c = 2,459,358.957820 ± 0.000086 BJD_TDB, b = 0.730 ± 0.012, a/R_* = 2.321 ± 0.031. The second set of MuSCAT3 transit photometry had the best photometric precision of all, yielding T_c = 2,459,391.90867 ± 0.00010 BJD_TDB, b = 0.766 ± 0.011, and a/R_* = 2.220 ± 0.032.

The individual MuSCAT2 and MuSCAT3 R_p /R_* values are listed at the top of Table 7. The detrended transit light curves are plotted in Figure 7. Notably, the MuSCAT3 #1 transit depths show a high level of achromaticity, with all four measurements lying within 1.1σ of each other; meanwhile, the MuSCAT2 and MuSCAT3 #2 transit depths show slightly larger variance (2.4σ and 2.0σ, respectively). The broad consistency in transit depths across the various photometric bands serves as supporting evidence against the false positive scenario of a blended eclipsing binary.

Zoom In Zoom Out Reset image size

Table 7. Individual Ground-based Transit Fit Results

Data Set	UT Date	Detrending Vectors a	Rp /R*	σ (ppm)
MuSCAT2 $g^{\prime}$	2021 Jun 12	...	0.0799 ± 0.0012	2555
MuSCAT2 $r^{\prime}$	2021 Jun 12	...	0.0834 ± 0.0011	2597
MuSCAT2 $i^{\prime}$	2021 Jun 12	...	0.0817 ± 0.0012	2113
MuSCAT2 zs	2021 Jun 12	...	0.0798 ± 0.0010	2914
MuSCAT3 $g^{\prime}$ #1	2021 May 24	...	0.0818 ± 0.0007	810
MuSCAT3 $r^{\prime}$ #1	2021 May 24	t, AM, Fcomp, ΔPSF	0.0821 ± 0.0011	811
MuSCAT3 $i^{\prime}$ #1	2021 May 24	t, AM, Fcomp, ΔPSF	0.0827 ± 0.0012	914
MuSCAT3 zs #1	2021 May 24	Fcomp, ΔPSF	0.0811 ± 0.0008	666
MuSCAT3 $g^{\prime}$ #2	2021 Jun 26	Fsky	0.0815 ± 0.0008	785
MuSCAT3 $r^{\prime}$ #2	2021 Jun 26	Fsky	0.0809 ± 0.0006	808
MuSCAT3 $i^{\prime}$ #2	2021 Jun 26	AM, Fsky	0.0829 ± 0.0008	942
MuSCAT3 zs #2	2021 Jun 26	ΔPSF, Fsky	0.0813 ± 0.0007	635
SSO B	2020 Aug 24	AM, Fcomp, ΔPSF, Fsky	0.0777 ± 0.0075	3711
ULMT $r^{\prime}$	2020 Jul 31	x, t, ΔPSF	0.0794 ± 0.0022	2580
WBRO R	2021 Apr 7	x, t	0.0795 ± 0.0025	1631
SAAO $i^{\prime}$	2021 Apr 7	AM, Fcomp	0.0810 ± 0.0025	2500
GdP $i^{\prime}$	2021 Apr 7	AM, Fcomp, ΔPSF	0.0811 ± 0.0018	2239
ULMT $i^{\prime}$	2021 Apr 9	Fcomp	0.0820 ± 0.0021	2123
MLO I	2020 Aug 2	t, AM, Fcomp, ΔPSF, Fsky	0.0860 ± 0.0066	4260
KeplerCam $z^{\prime}$ #2	2021 Apr 9	AM, Fcomp, ΔPSF, Fsky	0.0837 ± 0.0018	2610
KeplerCam $z^{\prime}$ #3	2021 Apr 11	...	0.0831 ± 0.0018	3291
SSO zs	2020 Aug 24	...	0.0811 ± 0.0033	2257

Note.

^a The optimized sets of detrending vectors used in the light-curve fits. Refer to the text for the definitions of variables.

Download table as:  ASCIITypeset image

The estimates of b and a/R_* from the MuSCAT2 and MuSCAT3 light-curve fits, which mutually agree at the 2.2σ level, have uncertainties that are almost an order of magnitude smaller than the corresponding errors derived from the TESS light curve alone (Table 5), highlighting the power of these high-cadence, high-S/N light curves in resolving the detailed transit geometry of the TOI-2109 system. For the remaining ground-based transit light curves, we used the most precise set of b and a/R_* values (from the MuSCAT3 #2 fit) as Gaussian priors. Priors on T_c were derived by interpolating the TESS, MuSCAT2, and MuSCAT3 timing measurements, assuming a linear orbital ephemeris. The full results of our individual ground-based transit light-curve fits are shown in Table 7. In addition to the measured transit depth, the optimized detrending vector set and best-fit uniform per-point scatter σ are provided for each light curve. The entries from non-MuSCAT facilities are sorted by bandpass. The corresponding systematics-corrected transit light curves are shown in Figure 8.

Zoom In Zoom Out Reset image size

For the K_s -band secondary eclipse fit, we utilized the same priors on orbital ephemeris and transit-shape parameters as in the ground-based transit light-curve fits; additionally, we fixed R_p /R_* to the weighted average of the individual MuSCAT3 #2 depths. Due to the relatively large planet–star flux contrast ratio in the K_s band and the sizable temporal baseline, we expect some variation in the out-of-eclipse flux due to the planet's atmospheric brightness modulation. This phase-curve signal can bias the measured eclipse depth if not accounted for in the analysis (e.g., Bell et al. 2019). We therefore included an additional quadratic function in time to remove any curvature from the planet's phase curve that is present in the out-of-eclipse data.

We measured an eclipse depth of 2027 ± 88 ppm and obtained a best-fit per-point uncertainty of 937 ppm. For the optimal set of detrending vectors, we utilized the two companion star fluxes F₁ and F₂, the airmass AM, and the sky background level F_sky. The bottom panel of Figure 2 shows the detrended light curve.

To obtain a preliminary mass measurement of TOI-2109b, we used the radvel package (Fulton et al. 2018) and modeled the planet's orbital RV signal assuming a circular orbit and no additional planets in the system. The assumption of e = 0 was motivated by the results of our TESS light-curve fits (Section 4.1.5). The full transit duration is roughly 0.15 in orbital phase, and we excluded all RV measurements that were obtained within 0.08 in orbital phase of the mid-transit time; these RVs are denoted in Table 3 by the superscripts. Due to the host star's rapid rotation, the precision of the RV measurements is quite poor, with some RV uncertainties as high as 1.5 km s⁻¹ or more. Nevertheless, the phase-folded RVs show good phase coverage across the two quadratures, and the RVs calculated from our spectroscopic transit observations provide high-S/N data points near the primary transit. For the final results presented in this paper, we removed all RV measurements with uncertainties greater than 1 km s⁻¹.

In this stand-alone analysis of the RV signal, we placed Gaussian priors on the mid-transit time T_c and orbital period P, using the median and 1σ values from the TESS phase-curve fit results (Table 5). In addition to the orbital ephemeris parameters and the RV semiamplitude K_p , we fit for the systemic RV offset and jitter of each instrument, γ_i and σ_jit,i . The parameter space was sampled using the default MCMC routine within radvel.

We obtained a 7.1σ detection of the planetary RV signal, with a semiamplitude of 850 ± 120 m s⁻¹. Including the RVs with uncertainties greater than 1 km s⁻¹ yielded a K_p value that agrees with the previously listed amplitude at better than the 0.1σ level. No significant long-term RV trends were measured when allowing for linear and quadratic temporal terms in the RV model. Using the stellar mass determined from our SED fitting (${M}_{* }={1.447}_{-0.078}^{+0.075}$ M_☉; Table 4) and the orbital parameters from the TESS phase-curve fit, we derived a planet mass of M_p = 4.77 ± 0.70 M_Jup. The phase-folded and offset-corrected RVs are plotted in Figure 9.

Zoom In Zoom Out Reset image size

Visual inspection of the residuals for the RVs taken during transit does not reveal any significant deviations. Such deviations can arise due to the Rossiter–McLaughlin (RM) effect, wherein the planet occults regions of the stellar disk with different rotational velocities and creates systematic aberrations in the stellar line shapes and resultant RVs (Rossiter 1924; McLaughlin 1924). The maximum value of the RV anomaly due to the RM effect is given by (e.g., Gaudi & Winn 2007; Albrecht et al. 2011)

$$\begin{eqnarray}&&{K}_{\mathrm{RM}}\sim {\left(\displaystyle \frac{{R}_{p}}{{R}_{* }}\right)}^{2}v\sin {i}_{* }\left(\sqrt{1-{b}^{2}}+b\tan \lambda \right)\cos \lambda ,\end{eqnarray} \tag{ 8 }$$

where λ is the sky-projected obliquity of the planet's orbit. Using the median values of R_p /R_*, b, and λ from our global analysis of the TESS photometry, ground-based light curves, RVs, and spectroscopic transit observations (Section 4.4; Table 8), we find K_RM = 0.37 km s⁻¹. This value is smaller than the average uncertainty of the RV measurements obtained during the planetary transit. We also note that the estimate in Equation (8) neglects the effects of limb and gravity darkening, which reduce the maximum RM anomaly, particularly for systems that are close to aligned. Therefore, we posit that the precision of the RVs, which is severely affected by the fast stellar rotation, is insufficient for securing the detection of an RM signal.

Table 8. Results from the Joint Fit to the TESS Photometry, Ground-based Light Curves, Radial Velocities, and Spectroscopic Transit Observations

Parameter	Description	Value
Fitted Parameters
Rp /R*	Planet–star radius ratio	0.08155 ± 0.00022
Tc	Mid-transit time (BJDTDB)	2,459,378.459370 ± 0.000059
P	Orbital period (day)	0.67247414 ± 0.00000028
b	Impact parameter	0.7481 ± 0.0073
a/R*	Scaled semimajor axis	2.268 ± 0.021
${\bar{f}}_{p}$	Average TESS-band relative planetary flux (ppm)	370 ± 41
Dd,K	Ks -band secondary eclipse depth (ppm)	2012 ± 80
Aatm	Atmospheric brightness modulation semiamplitude (ppm)	362 ± 19
ψ	Phase offset of the atmospheric brightness modulation (deg)	4.0 ± 2.3
Aellip	Ellipsoidal distortion semiamplitude (ppm)	245 ± 19
Mp	Planet mass (MJup)	5.02 ± 0.75
M* b	Stellar mass (M☉)	1.453 ± 0.074
Teff b	Stellar effective temperature (K)	6540 ± 160
${\gamma }_{\mathrm{TRES}}$	Radial-velocity offset for TRES (km s−1)	−25.64 ± 0.11
γFIES	Radial-velocity offset for FIES (km s−1)	−25.61 ± 0.21
${\sigma }_{\mathrm{jit},\mathrm{TRES}}$	Radial-velocity jitter for TRES (km s−1)	0.22 ± 0.15
σjit,FIES	Radial-velocity jitter for FIES (km s−1)	0.37 ± 0.24
$v\sin {i}_{* }$ b	Sky-projected stellar rotational velocity (km s−1)	81.2 ± 1.6
λ	Sky-projected obliquity (deg)	1.7 ± 1.7
vnonrot	Nonrotational broadening component (km s−1)	7.5 ± 1.6
Π1	First stellar variability period (day)	0.61395 ± 0.00055
α1	Sine semiamplitude at Π1 (ppm)	−9 ± 15
β1	Cosine semiamplitude at Π1 (ppm)	−153 ± 16
$\alpha {{\prime} }_{1}$	Sine semiamplitude at 1/2 × Π1 (ppm)	−68 ± 16
$\beta {{\prime} }_{1}$	Cosine semiamplitude at 1/2 × Π1 (ppm)	62 ± 15
Π2	Second stellar variability period (day)	0.9674 ± 0.0013
α2	Sine semiamplitude at Π2 (ppm)	24 ± 16
β2	Cosine semiamplitude at Π2 (ppm)	−222 ± 15
u1,TESS a	TESS-band quadratic limb-darkening coefficient	0.29 ± 0.03
u2,TESS a	TESS-band quadratic limb-darkening coefficient	0.20 ± 0.04
u1,B a	B-band quadratic limb-darkening coefficient	0.56 ± 0.05
u2,B a	B-band quadratic limb-darkening coefficient	0.25 ± 0.05
${u}_{1,g^{\prime} }$ a	$g^{\prime}$-band quadratic limb-darkening coefficient	0.53 ± 0.03
${u}_{2,g^{\prime} }$ a	$g^{\prime}$-band quadratic limb-darkening coefficient	0.18 ± 0.04
${u}_{1,r^{\prime} }$ a	$r^{\prime}$-band quadratic limb-darkening coefficient	0.34 ± 0.03
${u}_{2,r^{\prime} }$ a	$r^{\prime}$-band quadratic limb-darkening coefficient	0.23 ± 0.04
u1,R a	R-band quadratic limb-darkening coefficient	0.40 ± 0.05
u2,R a	R-band quadratic limb-darkening coefficient	0.22 ± 0.05
${u}_{1,i^{\prime} }$ a	$i^{\prime}$-band quadratic limb-darkening coefficient	0.32 ± 0.03
${u}_{2,i^{\prime} }$ a	$i^{\prime}$-band quadratic limb-darkening coefficient	0.23 ± 0.04
u1,I a	I-band quadratic limb-darkening coefficient	0.30 ± 0.05
u2,I a	I-band quadratic limb-darkening coefficient	0.20 ± 0.05
${u}_{1,z^{\prime} /{z}_{s}}$ a	$z^{\prime} /{z}_{s}$-band quadratic limb-darkening coefficient	0.29 ± 0.03
${u}_{2,z^{\prime} /{z}_{s}}$ a	$z^{\prime} /{z}_{s}$-band quadratic limb-darkening coefficient	0.19 ± 0.04
Derived Parameters
Rp	Planet radius (RJup)	1.347 ± 0.047
a	Semimajor axis (au)	0.01791 ± 0.00065
i	Orbital inclination (deg)	70.74 ± 0.37
q	Planet–star mass ratio	0.00331 ± 0.00052
Kp	Radial-velocity semiamplitude (km s−1)	0.86 ± 0.13
$\mathrm{log}{g}_{p}$	Planet surface gravity (cgs)	3.836 ± 0.071
ADopp	Doppler-boosting semiamplitude (ppm)	8.6 ± 1.3
${\left({R}_{p}/{R}_{* }\right)}^{2}$	Transit depth (ppm)	6651 ± 36
Dd,TESS	TESS-band secondary eclipse depth (ppm)	731 ± 46
Dn,TESS	TESS-band nightside flux (ppm)	9 ± 43
Tirr c	Irradiation temperature (K)	4340 ± 100
Teq c	Dayside-redistribution equilibrium temperature (K)	3646 ± 88
Tb,day c	Dayside brightness temperature (K)	3631 ± 69
Tb,night c	Nightside brightness temperature (K)	<2500

Notes.

^a These parameters were constrained by priors based on the tabulated limb-darkening coefficients from Claret (2017). See Table 6. ^b These parameters were constrained by priors based on modeling of the host star's SED and TRES spectra: M_* = 1.447 ± 0.077 M_☉, T_eff = 6530 ± 160 K, and $v\sin {i}_{* }=81.9\pm 1.7$ km s⁻¹. ^c The irradiation temperature is defined as ${T}_{\mathrm{irr}}\equiv {T}_{\mathrm{eff}}\sqrt{{R}_{* }/a}$. The dayside-redistribution equilibrium temperature assumes zero Bond albedo, a uniform dayside temperature, and T = 0 K on the nightside. The dayside brightness temperature is derived from a joint blackbody fit to the TESS- and K_s -band secondary eclipse depths, assuming zero geometric albedo. For the nightside temperature, the 2σ upper limit is given. See Section 5.1.

Download table as:  ASCIITypeset images: 1 2

To obtain the final results from our analysis of the TOI-2109 system, we carried out a joint fit of all available data sets—TESS photometry, ground-based transit and secondary eclipse light curves, RV measurements, and spectroscopic transits—to simultaneously measure all of the astrophysical quantities of interest. Given the mutual consistency between the individually measured planet–star radius ratios from the TESS and ground-based transit fits (Tables 5 and 7), we defined a single R_p /R_* parameter for all data sets.

By fitting the RV measurements jointly with the TESS light curve, we were able to separate the Doppler-boosting signal from the total phase-curve modulation at the sine of the orbital frequency and retrieve the phase offset ψ in the planet's atmospheric brightness modulation (see discussion in Section 4.1.1). Both the RV semiamplitude K_RV and the Doppler-boosting semiamplitude A_Dopp depend on the planet mass M_p (e.g., Loeb & Gaudi 2003; Perryman 2011; Shporer 2017):

$$\begin{eqnarray}&&{A}_{\mathrm{Dopp}}={\left[\displaystyle \frac{2\pi G}{{{Pc}}^{3}}\displaystyle \frac{{M}_{p}^{3}{\sin }^{3}i}{{({M}_{* }+{M}_{p})}^{2}}\right]}^{1/3}{\left\langle \displaystyle \frac{{{xe}}^{x}}{{e}^{x}-1}\right\rangle }_{\mathrm{TESS}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{K}_{\mathrm{RV}}={\left(\displaystyle \frac{2\pi G}{P}\right)}^{1/3}\displaystyle \frac{{M}_{p}\sin i}{{({M}_{* }+{M}_{p})}^{2/3}}.\end{eqnarray} \tag{ 10 }$$

In Equation (9), x ≡ hc/k λ T_eff, and the term in the angled brackets is integrated with respect to wavelength λ across the TESS bandpass, weighted by the transmission function of the instrument. The expression in Equation (10) assumes zero orbital eccentricity.

Instead of fitting for A_Dopp and K_RV independently, we included the planet mass M_p , stellar mass M_s , and stellar effective temperature T_eff as free parameters and self-consistently modeled both the RV trend and the Doppler-boosting signal. The stellar mass and effective temperature were constrained by Gaussian priors based on the results of the SED modeling (Table 4). All other orbital ephemeris, transit-shape, limb-darkening, phase-curve, stellar variability, and RV parameters were treated in an identical manner to the corresponding analyses of individual data sets (Sections 4.1–4.3).

The ellipsoidal distortion amplitude is also dependent on M_p . However, unlike in the case of Doppler boosting, the physical processes driving the stellar tidal response are strongly contingent upon the detailed characteristics of the stellar interior and atmosphere. Secondary effects from the stellar rotation and the interaction between the external tidal force and pulsation modes can often lead to significant discrepancies between the predicted behavior and the measured amplitude (see, for example, Burkart et al. 2012 and Wong et al. 2020c). Given this caveat, we did not model the ellipsoidal distortion using the planet mass in our joint fit, but instead kept the ellipsoidal distortion amplitude A_ellip as an independent fit parameter. A comparison of the predicted and measured amplitudes is provided in Section 5.2.

When modeling the spectroscopic transit observations, we followed the methodology of Zhou et al. (2016). Prior to fitting, the average line-broadening profile derived from the out-of-transit spectra (Sections 2.5.1) was subtracted from the full set of line-broadening profiles obtained during the two spectroscopic transits. The resultant differential rotational profiles from both nights of observation are displayed in the top panel of Figure 10. The Doppler shadow of the transiting planet was modeled as a Gaussian centered at v_p (t)—the projected rotational velocity of the region occulted by the planet; the width contains contributions from the spectral resolution of the instrument and a nonrotational broadening component v_nonrot, which covers both instrumental broadening and macroturbulence in the stellar atmosphere. The area of the signal at each timestamp is equivalent to 1 − λ_t (t), where λ_t (t) is the transit light curve evaluated at time t.

Zoom In Zoom Out Reset image size

The location and orientation of the Doppler transit signal depend on the sky-projected stellar rotational velocity $v\sin {i}_{* }$, sky-projected orbital obliquity λ, and impact parameter b. These parameters, in addition to the nonrotational broadening component v_nonrot, were included as free parameters in the joint fit, with $v\sin {i}_{* }$ constrained by a Gaussian prior based on the TRES-derived measurement: $v\sin {i}_{* }=81.9\pm 1.7$ km s⁻¹ (Section 2.5.1).

The total number of free parameters in the joint fit is 44. The photometric light curves were corrected for instrumental systematics prior to fitting, with the uniform per-point flux uncertainties fixed to their respective best-fit values from the individual analyses (see Tables 5 and 7). For the TESS light curve, the red-noise-inflated per-point uncertainty σ_r was used.

4.4.1. Results

The results of our joint MCMC fit are shown in Table 8; the values of relevant derived parameters are also provided. The best-fit full-orbit phase-curve, TESS- and K_s -band secondary eclipse, RV, and Doppler transit models are displayed in Figures 6, 9, 10, and 11, alongside the corresponding residuals.

Zoom In Zoom Out Reset image size

TOI-2109b has an orbital period of 0.67247414 ± 0.00000028 days and a radius of 1.347 ± 0.047 R_Jup. The addition of the ground-based transits in the joint fit greatly enhanced the precision of the orbital ephemeris and transit-shape parameters when compared to the TESS-only fit results in Table 5. Using the measured impact parameter and scaled semimajor axis—b = 0.7481 ± 0.0073 and a/R_* = 2.268 ± 0.021—we derived an orbital inclination of i = 7074 ± 037, indicating a moderately inclined viewing geometry.

A notable result is that the projected spin–orbit misalignment is consistent with zero. The tight constraints on the transit shape and the mid-transit time at the epoch of the spectroscopic transit observations yielded an extremely precise measurement of the projected orbital obliquity: λ = 17 ± 17). The faint Doppler shadow of the planet is discernible in the line-broadening profiles (Figure 10), which definitively confirms that the planet is orbiting the host star and excludes all blended scenarios. The low obliquity of TOI-2109b stands in stark contrast to the sizable subset of ultrahot Jupiters with large obliquities, including HAT-P-7b (1825 ± 94; Narita et al. 2009; Winn et al. 2009), KELT-9b (−88° ± 15°; Ahlers et al. 2020), WASP-33b (−1088 ± 10; Collier Cameron et al. 2010), and WASP-121b (−1022 ± 54; Delrez et al. 2016). Meanwhile, other ultrahot Jupiters with near-zero obliquities include HATS-70b (89 ± 51; Zhou et al. 2019), WASP-18b (4° ± 5°; Triaud et al. 2010), and WASP-19b (10 ± 12; Southworth & Mancini 2016).

From the joint fit, we obtained M_p = 5.02 ± 0.75 M_Jup, which is consistent with the mass estimate we obtained from the dedicated RV fit in Section 4.3 (4.77 ± 0.70 M_Jup); the derived planet surface gravity is $\mathrm{log}{g}_{p}=3.836\pm 0.071$ (in cgs units). When combined with the orbital parameters via Equation (9), the planet mass corresponds to a predicted Doppler-boosting amplitude of A_Dopp = 8.6 ± 1.3 ppm.

From the TESS light curve, we measured a high-S/N secondary eclipse depth of 731 ± 46 ppm and a nightside flux consistent with zero. The semiamplitude of the planet's atmospheric brightness modulation is A_atm = 362 ± 19 ppm, with a marginal eastward offset of ψ = 40 ± 23. The amplitude of the stellar ellipsoidal distortion signal is A_ellip = 245 ± 19 ppm. The two robustly detected phase-curve components are plotted separately in the middle panel of Figure 6. All of the phase-curve parameter values from our joint fit are statistically identical to the results we obtained from fitting the TESS light curve alone. Meanwhile, the K_s -band secondary eclipse depth is D_d,K = 2012 ± 80 ppm.

The two fitted stellar variability periods are Π₁ = 0.61395 ± 0.00055 days and Π₂ = 0.9674 ± 0.0013 days. Combining the sine and cosine coefficients at each period into a single value, we find that the variability amplitudes at the fundamentals are 154 ± 16 ppm and 224 ± 15 ppm, respectively. The total first harmonic amplitude at Π₁ is 93 ± 15 ppm. The best-fit stellar variability signals are shown in the right panels of Figure 4, phase-folded on the corresponding periods.

Given the rapid stellar rotation, we checked the MuSCAT3 $g^{\prime}$-band transit light curves for asymmetry that could arise due to the stellar oblateness and gravity darkening (e.g., Barnes et al. 2011). These light curves were chosen on account of their high S/N and the more pronounced gravity darkening of the star at bluer wavelengths. A measurement of such asymmetry can yield the full three-dimensional stellar rotation axis and the true obliquity of the system. Following the methodology of Ahlers et al. (2020), we used the formalism of Espinosa Lara & Rieutord (2011) to model the gravity-darkened transit, with the stellar inclination i_* and rotation period P_rot as input parameters. A prior on P_rot based on the measured $v\sin {i}_{* }$ was applied. The near-zero projected obliquity of the system makes significant asymmetries unlikely, with the effects of gravity darkening becoming largely degenerate with the impact parameter and limb-darkening coefficients. We indeed did not detect any statistically significant transit asymmetry in the light-curve modeling, and as such no constraints on the stellar inclination were obtained from the gravity-darkening analysis.

The detection of periodic stellar variability presents another method for estimating the stellar inclination. From the measured sky-projected rotational velocity $v\sin {i}_{* }$ and the stellar radius R_*, a given variability period Π can be converted to stellar inclination via ${i}_{* }={\sin }^{-1}({\rm{\Pi }}v\sin {i}_{* }/2\pi {R}_{* })$. Plugging in Π₁ and Π₂, we obtained two independent stellar inclination estimates: i_*,1 = 36° ± 2° and i_*,2 = 67° ± 6°. Given the roughly zero sky-projected obliquity and the measured orbital inclination of i = 7074 ± 037, we consider the second scenario, which corresponds to full spin–orbit alignment, to be more plausible. However, it should be acknowledged that the observed stellar variability need not be a manifestation of the star's rotation (e.g., from starspots) and may instead stem from intrinsic pulsation modes. In the latter case, the variability period is entirely unrelated to the measured sky-projected rotational velocity.

Assuming i_*,2 = 67° ± 6°, we computed a stellar oblateness of 1 − R_*,pole/R_*,eq ≃ 0.02. When accounting for the viewing geometry, the mean radius of the sky-projected disk is R_*,mean ≃ 1.716 R_☉, which differs from the R_* value derived from the SED fitting analysis by just ∼1% and is well within the 1σ confidence region. Therefore, the stellar oblateness is not expected to significantly bias the derived R_p and a values listed in Table 8.

We obtained secondary eclipse observations of TOI-2109b in the TESS and K_s bandpasses. The corresponding depths measured from the joint fit are 731 ± 46 and 2012 ± 80 ppm, respectively. Planets on ultrashort orbits are expected to be tidally locked, and so the depth of the secondary eclipse corresponds to the total brightness of the fixed dayside hemisphere relative to the stellar flux. When considering only thermal emission from the planet's atmosphere, we can express the relative planetary flux as

$$\begin{eqnarray}&&D={\left(\displaystyle \frac{{R}_{p}}{{R}_{* }}\right)}^{2}\displaystyle \frac{\int {B}_{p,\lambda }({T}_{b})\tau (\lambda )d\lambda }{\int {B}_{* ,\lambda }({T}_{\mathrm{eff}})\tau (\lambda )d\lambda },\end{eqnarray} \tag{ 11 }$$

where B_p,λ and B_*,λ are the emission spectra of the planet and star, respectively, and τ(λ) is the transmission function of the bandpass. The variable T_b denotes the brightness temperature of the planet in the bandpass.

To calculate the dayside brightness temperature of TOI-2109b, we assumed that the planet's emission spectrum is that of a blackbody. For the stellar flux, we utilized PHOENIX models (Husser et al. 2013) and calculated the band-integrated flux (i.e., the denominator in the previous equation) for a grid of models spanning the measured stellar parameters for TOI-2109 (Table 4). Then, following the technique described in Wong et al. (2020e), we fit for a polynomial function in (T_eff, $\mathrm{log}{g}_{* }$, [Fe/H]) that interpolates the grid points, enabling smooth sampling of the integrated stellar flux. The posterior distribution of the planet's brightness temperature was calculated using emcee, with Gaussian priors on D and R_p /R_*, and stellar parameters from the joint fit and the SED analysis.

From the TESS-band secondary eclipse depth, we derived a dayside brightness temperature of T_b,TESS = 3729 ± 82 K. A similar calculation with the K_s -band secondary eclipse depth yielded T_b,K = 3518 ± 81 K. Given the broad consistency between the TESS- and K_s -band brightness temperatures, we also carried out a joint fit to both secondary eclipse depths and obtained a dayside brightness temperature of T_b,day = 3631 ± 69 K. Only KELT-9b has a higher measured dayside temperature: 4566 ± 138 K (Bell et al. 2021). Through an analogous calculation, we converted the measured TESS-band nightside flux to an upper limit on the nightside brightness temperature. Due to the highly nonlinear relationship between the planet–star contrast ratio and the brightness temperature in the optical, we were only able to place a broad constraint: T_b,night < 2500 K (2σ).

The previous discussion ignored the possibility of reflected starlight across the dayside hemisphere. A nonzero geometric albedo A_g contributes an amount equal to ${A}_{g}\times {\left({R}_{p}/a\right)}^{2}$ to the relative planetary flux, which serves to lower the dayside brightness temperature needed to match the TESS-band flux. While the relatively broad constraints on the dayside brightness temperature from the K_s -band secondary eclipse alone and the marginally higher TESS-band dayside brightness temperature formally allow for a wide range of optical geometric albedos (up to ∼0.2 at 1σ), we consider the scenario of significant dayside reflectivity unlikely. The extremely high temperatures across the dayside of TOI-2109b preclude the formation of condensate clouds or hazes, as all major condensate species are expected to be in the vapor phase, or even dissociated (e.g., Wakeford et al. 2017; Lothringer et al. 2018). Indeed, higher-precision optical and near-infrared secondary eclipse measurements of other ultrahot Jupiters have revealed very low geometric albedos consistent with zero (see, for example, Bell et al. 2019; Shporer et al. 2019; Wong et al. 2020d, and references therein). Moreover, an apparent excess in planetary flux at optical wavelengths may indicate deviations in the dayside emission spectrum from a simple blackbody, which can be caused by high-temperature opacity sources such as H⁻ (e.g., Arcangeli et al. 2018). We address this possibility with detailed emission spectrum modeling of TOI-2109b in Section 5.4.

We can compare the measured dayside brightness temperature to the equilibrium temperature T_eq (Cowan & Agol 2011):

$$\begin{eqnarray}&&{T}_{{\rm{eq}}}={T}_{{\rm{eff}}}\sqrt{\displaystyle \frac{{R}_{\ast }}{a}}\xi {(1-{A}_{{\rm{B}}})}^{1/4}.\end{eqnarray} \tag{ 12 }$$

The first two terms are referred to as the irradiation temperature ${T}_{\mathrm{irr}}\equiv {T}_{\mathrm{eff}}\sqrt{{R}_{* }/a};$ using the parameters for the TOI-2109 system, we obtained T_irr = 4340 ± 100 K. The variable A_B represents the Bond albedo, while ξ is a factor that reflects the efficiency of day–night heat recirculation. For the dayside-redistribution equilibrium temperature, we assumed zero Bond albedo, a uniform dayside temperature, and no heat transport to the nightside atmosphere, which correspond to A_B = 0 and ξ = (1/2)^1/4. This equilibrium temperature (T_eq = 3646 ± 88 K) is consistent with the measured dayside brightness temperature (T_day = 3631 ± 69 K). Meanwhile, the maximum theoretical temperature, which assumes no heat circulation across the dayside atmosphere (ξ = (2/3)^1/4) is ${T}_{\max }\sim 3900$ K, significantly hotter than T_day.

We now place the dayside thermal energy budget of TOI-2109b in the context of other hot Jupiters. Bell et al. (2021) carried out a uniform analysis of all available Spitzer 4.5 μm phase curves and derived T_day for 16 planets based on their secondary eclipse depths. In Figure 13, we plot the ratio T_day/T_irr as a function of T_eq for these targets, alongside the corresponding values for TOI-2109b and the two recently discovered ultrahot Jupiters from TESS: TOI-1431b and TOI-1518b. The horizontal lines denote two benchmark scenarios: the dayside-redistribution equilibrium temperature, with a uniform dayside temperature and T = 0 K on the nightside (i.e., ξ = (1/2)^1/4), and the maximum theoretical dayside temperature (i.e., the no-recirculation limit, ξ = (2/3)^1/4).

Zoom In Zoom Out Reset image size

The first notable observation is that TOI-2109b lies intermediate between the dense cluster of well-characterized gas giants with T_eq ∼ 3000 K and the extreme end-member KELT-9b. It follows that intensive study of TOI-2109b's atmosphere may help connect the unique properties of KELT-9b to the broader underlying thermophysical processes driving the atmospheric dynamics and chemistry of hot gas giants. Next, we find a number of planets with T_day < T_eq, particularly at the ends of the temperature range shown. Such a scenario could be indicative of a nonzero Bond albedo and/or a significant amount of day–night heat transport. When examining the scatter plot for overall trends, we see that the T_day/T_irr value of TOI-2109b lies in between the average value of the cluster of gas giants with T_eq ∼ 3000 K, which are consistent with extremely low levels of day–night heat recirculation, and the lower value of KELT-9b. This behavior suggests a tenuous decreasing trend that is predicted by recent atmospheric modeling of ultrahot Jupiters, which has demonstrated that the dissociation of molecular hydrogen on the hot dayside and its recombination on the cooler nightside can efficiently transport energy across the planet's surface, resulting in lower dayside temperatures and reduced day–night temperature contrasts in the hottest exoplanets (e.g., Bell & Cowan 2018; Komacek & Tan 2018; Parmentier et al. 2018; Tan & Komacek 2019; Roth et al. 2021).

Past models simulating the atmospheric circulation of ultrahot Jupiters similar to TOI-2109b have shown that significant dayside hotspot offsets (∼20°) are expected under conditions of weak atmospheric drag and short planetary rotation period (Tan & Komacek 2019). While the lack of a large eastward phase offset in the measured visible phase curves of cooler gas giants may be attributable to the contribution of advected reflective clouds from the nightside crossing over the western terminator, the extremely high dayside temperature of TOI-2109b renders the possibility of partial cloud cover highly unlikely. Instead, magnetohydrodynamic forces on the dayside atmosphere may be enhancing the atmospheric drag and reducing the magnitude of the hotspot offset (e.g., Rogers & Komacek 2014). Observational tests of the detailed atmospheric circulation can be achieved through spectroscopic thermal phase curves, which we discuss briefly in Section 5.4.

We measured a strong ellipsoidal distortion signal in the TESS light curve at almost 13σ significance. Having obtained the amplitude of this modulation independent of the constraint on planet mass from the RVs, we can now compare the measured value to the theoretical prediction. The expected ellipsoidal distortion amplitude at the first harmonic of the cosine of the orbital phase is related to the planet–star mass ratio q via the following expression (e.g., Morris 1985; Morris & Naftilan 1993; Shporer 2017):

$$\begin{eqnarray}&&{A}_{\mathrm{ellip}}={\alpha }_{\mathrm{ellip}}q{\left(\displaystyle \frac{{R}_{* }}{a}\right)}^{3}{\sin }^{2}i.\end{eqnarray} \tag{ 13 }$$

The prefactor α_ellip is a term that depends on the linear limb- and gravity-darkening coefficients u and g (see, for example, Morris 1985 and Shporer 2017). We defined Gaussian priors centered on the values tabulated in Claret (2017): u = 0.4497 and g = 0.2273; we used 0.05 for the width of the prior on u, while we selected 0.10 for g, given the steeper correlation between g and T_eff for F-type stars. After carrying out Monte Carlo sampling of the parameter distributions, we arrived at a predicted ellipsoidal distortion semiamplitude of 281 ± 52 ppm.

This prediction is statistically consistent with the measured value from our joint fit (A_ellip = 245 ± 19 ppm) and shows that the classical theory of stellar tidal deformation accurately describes the gravitational response of TOI-2109 to its orbiting companion. This agreement holds despite the rapid rotation of the host star, which is not accounted for in the theoretical formalism underpinning Equation (13).

The ultrashort orbit of TOI-2109b raises the possibility that the planet may be close to being broken apart by tidal forces. The minimum distance from the star at which an orbiting companion can reside without being catastrophically disrupted is called the Roche limit. For gaseous planets, this quantity can be converted to a minimum orbital period ${P}_{min}\simeq 0.40\,{\rm{days}}/\sqrt{{\rho }_{p}}$, where ρ_p is the bulk density of the planet in units of g cm⁻³ (Rappaport et al. 2013). For the TOI-2109 system, we calculated ${P}_{\min }\simeq 0.12$ days, which is significantly shorter than the measured orbital period (∼0.67 days). Therefore, we conclude that the high planet mass ensures that TOI-2109b is not threatened by tidal disruption, despite its extremely close orbit. For comparison, the roughly Jupiter-mass planets WASP-12b, WASP-19b, and WASP-121b all have orbital periods that are within three times their respective ${P}_{\min }$ limits.

Even when a planet is not close to the Roche limit, the strong tidal forces can lead to significant deformation of the planet's atmosphere and even mass loss through tidal stripping. A commonly used metric for evaluating the likelihood of atmospheric mass loss is the Roche lobe filling factor r ≡ R_p /R_Roche, where R_Roche is the radius of the effective spherical surface around the planet within which material is gravitationally bound to the planet. This ratio can be approximated as (Eggleton 1983)

$$\begin{eqnarray}&&r=\left(\displaystyle \frac{{R}_{p}}{{R}_{* }}\right){\left(\displaystyle \frac{a}{{R}_{* }}\right)}^{-1}\displaystyle \frac{0.6{q}^{2/3}+\mathrm{ln}(1+{q}^{1/3})}{0.49{q}^{2/3}}.\end{eqnarray} \tag{ 14 }$$

For TOI-2109b, we used the values of q, a/R_*, and R_p /R_* from the joint fit and obtained r = 0.50 ± 0.03.

The planet radius used in the previous calculation corresponds to the radial distance to which the atmosphere is optically thick at visible wavelengths. A low-density exosphere might still extend far beyond the planet's surface. For comparison, the two notable Jupiter-mass planets experiencing detectable atmospheric mass loss through Roche lobe overflow—WASP-12b (e.g., Haswell et al. 2012; Fossati et al. 2013; Bell et al. 2019) and KELT-9b (e.g., Yan & Henning 2018; Wyttenbach et al. 2020)—have r ∼ 0.8 and r ∼ 0.5, respectively. Meanwhile, outflow of metastable helium has been reported for WASP-107b (Spake et al. 2018; Allart et al. 2019), a hot sub-Saturn with a Roche lobe filling factor of just r ∼ 0.3. The specific conditions necessary for large-scale atmospheric mass loss and the corresponding properties of the outflowing gas are still poorly understood. The combination of high surface gravity and high atmospheric temperature makes TOI-2109b a particularly compelling target for detecting signatures of gaseous outflow and evaluating models of tide- and irradiation-driven atmospheric mass loss.

The strong gravitational interaction and tidal dissipation in ultrashort-period systems can cause the orbital period to measurably decay on year-long or decade-long timescales. The predicted rate of orbital decay depends on the host star's modified tidal quality factor $Q{{\prime} }_{* }$ (Goldreich & Soter 1966; Rasio et al. 1996; Sasselov 2003):

$$\begin{eqnarray}&&\dot{P}=-\displaystyle \frac{27\pi q}{2Q{{\prime} }_{* }}{\left(\displaystyle \frac{{R}_{* }}{a}\right)}^{5}.\end{eqnarray} \tag{ 15 }$$

Studies of tidal dissipation in late-type main-sequence stellar binaries and planet-hosting stars indicate that typical values of $Q{{\prime} }_{* }$ for F-type stars lie in the range 10⁵–10⁷ (e.g., Ogilvie & Lin 2007; Lanza et al. 2011). Using values of q and a/R_* from our joint fit (Table 8), we found orbital decay rates that span ∼10–740 ms yr⁻¹.

For comparison, the most recently published measurement of orbital decay rate for WASP-12b is 32.5 ± 1.6 ms yr⁻¹ (Turner et al. 2021); that system is the only incontrovertible case of tidal orbital decay hitherto discovered. The inferred stellar tidal quality factor for the late-F-type host star WASP-12 is $Q{{\rm{{\prime} }}}_{\ast }=(1.39\pm 0.15)\times {10}^{5}$, which is consistent with the lower end of the previously cited range.

Figure 14 shows the predicted orbital decay rate for all known hot Jupiters, assuming $Q{{\prime} }_{* }={10}^{5}$. For a given $Q{{\prime} }_{* }$, TOI-2109b has by far the fastest expected orbital decay rate among planetary mass companions, making this system the most promising candidate for probing orbital evolution and constraining the stellar tidal quality factor in the coming years. In second place is the supermassive gas giant WASP-18b, which also orbits an F-type host star. Curiously, no statistically significant variation in orbital period has been detected for WASP-18b, despite its high mass and the availability of high-precision transit-timing measurements spanning more than a decade (Patra et al. 2020). If $Q{{\prime} }_{* }$ of TOI-2109 is comparable to that of WASP-12, and the cadence and precision of follow-up transit timings are similar to what was obtained for WASP-12b, then we may expect to detect a nonlinear ephemeris within 2–3 yr. TESS is slated to reobserve this system during the extended mission in Sector 52 (currently scheduled for 2022 May 18–June 13), at which point there should be sufficient time baseline to detect orbital decay rates at the upper end of the aforementioned range.

Zoom In Zoom Out Reset image size

In addition to long-term transit monitoring, collecting further secondary eclipse timings will be crucial for distinguishing between tidal orbital decay and apsidal precession (e.g., Yee et al. 2020): in the former case, the time interval between consecutive events decreases in tandem for both transits and secondary eclipses, while in the latter case, the directions of the timing deviations proceed in contrary motion. Likewise, continued RV monitoring is necessary to check for line-of-sight systemic acceleration due to the presence of a long-period bound companion (i.e., the Rømer effect).

The brightness of the host star (V = 10.3 mag, K = 9.1 mag) and the highly irradiated dayside atmosphere make TOI-2109b very amenable to intensive atmospheric characterization with current and near-future facilities. The rapid stellar rotation and associated rotational broadening of the stellar lines are conducive to ground-based high-dispersion spectroscopy, which has been used extensively in recent years to detect individual atomic, ionic, and molecular species in the atmospheres of ultrahot Jupiters (e.g., Hoeijmakers et al. 2018; Casasayas-Barris et al. 2019). Emission spectroscopy, both in eclipse and across the full orbital phase, is particularly promising. While the high surface gravity of TOI-2109b ($\mathrm{log}{g}_{p}=3.836\pm 0.071$) severely depresses the atmospheric scale height and therefore attenuates any spectral features observed in transmission relative to lower-mass hot Jupiters, the prospects for emission studies are not affected by the surface gravity.

We can quantify a planet's amenability to emission spectroscopy by calculating the Emission Spectroscopy Metricemission spectroscopy metric (ESM; Kempton et al. 2018). We tailored the formulation of the ESM to be more suitable for hot gas giants by evaluating the planet–star contrast ratio in the K band, instead of at 7.5 μm as originally defined. Figure 15 shows the ESM values for our sample of hot Jupiters. TOI-2109b has the seventh highest ESM, behind HIP 65Ab (an inflated hot Jupiter around a K-dwarf with R_p /R_s ∼ 0.25; Nielsen et al. 2020), three ultrahot Jupiters orbiting A-stars (WASP-33b, KELT-9b, and KELT-20b), HD 189733b, and WASP-121b.

Zoom In Zoom Out Reset image size

To explore the prospects for atmospheric characterization in detail, we generated model dayside emission, nightside emission, and transmission spectra of TOI-2109b using the open-source radiative transfer code HELIOS (Malik et al. 2017, 2019). This model framework takes as input the stellar spectrum (interpolated from PHOENIX models) and the planet–star radius ratio and self-consistently computes the planet's temperature–pressure (TP) profile and emission/transmission spectra under radiative–convective equilibrium. Atmospheric opacities are calculated using the HELIOS-K module (Grimm & Heng 2015; Grimm et al. 2021), while the equilibrium chemistry is computed using FastChem (Stock et al. 2018). Line lists for over 600 species are included, ⁵³ along with continuum opacities from H⁻ (John 1988) and collision-induced absorption from H₂–H₂, H₂–He, and H–He (Abel et al. 2011; Karman et al. 2019). Following the analogous application of HELIOS to modeling the atmosphere of KELT-9b (Wong et al. 2020e), we fixed the atmospheric metallicity to solar, set the Bond albedo to zero, and assumed uniform dayside, nightside, and terminator TP profiles. To produce different versions of the spectra for comparison, we varied the day–night heat recirculation parameter (ranging from 0 to 1, as defined in Cowan & Agol 2011), nightside interior temperature T_night, and terminator temperature T_term.

Figure 16 shows a compilation of HELIOS model spectra for TOI-2109b assuming various values of , T_night, and T_term. The dayside emission spectra are largely featureless throughout the wavelength range 0.5–5.3 μm. The notable trend is the increasing H⁻ continuum opacity with decreasing (i.e., increasing dayside temperature), which is manifested by the emergence of the lobe-shaped excess emission feature shortward of ∼1.8 μm. This behavior has been identified and extensively characterized in previous theoretical work (e.g., Kitzmann et al. 2018; Lothringer et al. 2018; Parmentier et al. 2018). The measured TESS-band and K_s -band secondary eclipse depths are both consistent with the = 0.3 model, indicating a moderate level of dayside–nightside heat recirculation. The H⁻ emission feature helps explain the somewhat higher TESS-band brightness temperature when compared to the K_s -band value (Section 5.1). The = 0.3 model has an effective temperature of ∼3700 K, in agreement with the measured dayside brightness temperature T_day = 3631 ± 69 K.

Zoom In Zoom Out Reset image size

Given the poorly constrained nightside brightness temperature of TOI-2109b from our TESS-band nightside flux measurement (T_night < 2500 K at 2σ), we considered a wide range of plausible T_night values when generating the HELIOS emission spectra. As illustrated in the middle panel of Figure 16, even moderate-precision spectroscopy in the near-infrared region (δ ∼ 100 ppm at ∼0.1 μm resolution) is expected to reveal prominent absorption features due to H₂O (∼1.4–1.6, 1.7–2.2, 2.4–3.4 μm) and CO (∼4–5 μm) for T_night values greater than ∼1500 K. Detecting and modeling these spectral features at high S/N and fine spectral resolution will provide robust constraints on the atmospheric metallicity, C/O, and O/H ratios, with broad implications for the formation and evolution of TOI-2109b (see, for example, the review by Madhusudhan et al. 2016).

The bottom panel of Figure 16 shows a set of model transmission spectra for TOI-2109b ranging in T_term from 1000 to 3000 K; all of the curves are normalized to match the measured TESS-band transit depth. Generally, the molecular absorption features observed in the nightside emission spectra are also seen in transmission, though with significantly smaller amplitudes of ∼30 ppm or less. The 3000 K model has less prominent H₂O absorption features in the near-infrared, reflecting the marked decrease in water abundance due to dissociation at high temperatures. Meanwhile, the negative spectral slope throughout the red-optical and near-infrared is primarily determined by the wavelength-dependent H⁻ opacity (John 1988). Although the high surface gravity of TOI-2109b leads to relatively weak transmission features, the absorption features in the model spectra are nonetheless detectable by instruments such as NIRSpec on the James Webb Space Telescope (JWST). Moreover, for limb temperatures above 2000–2500 K, cloud formation is strongly disfavored (e.g., Wakeford et al. 2017; Lothringer et al. 2018), allowing for an unfettered view into the atmosphere. Obtaining a snapshot of the high-altitude chemistry and TP profile along the day–night terminator can provide additional insight into the overall atmospheric composition and dynamics.

In recent years, near-infrared emission spectroscopy across the entire orbital phase has been successfully carried out using the Hubble Space Telescope (HST) for a small handful of exoplanets, including WASP-18b (Arcangeli et al. 2019), WASP-43b (Stevenson et al. 2014), and WASP-103b (Kreidberg et al. 2018). The ability to track the composition and TP profile of the atmosphere at all longitudes makes these studies immensely impactful for our understanding of global atmospheric properties. These data sets enable more sophisticated modeling to constrain the detailed characteristics of day–night atmospheric circulation and chemical gradients across the surface. Looking to the near future, the enhanced capabilities of JWST have placed spectroscopic phase curves at the forefront of efforts at atmospheric characterization for gas-giant exoplanets. TOI-2109b is an extremely promising candidate for phase-resolved emission spectroscopy. In particular, the planet's equilibrium temperature makes it an attractive target for probing fundamental trends in gas-giant atmospheric dynamics.

Figure 17 shows the measured day–night brightness temperature contrast 1 − T_night/T_day for hot and ultrahot Jupiters, alongside theoretical predictions from the three-dimensional global circulation model (GCM) MITgcm, as computed in Tan & Komacek (2019) and Mansfield et al. (2020). Two different levels of atmospheric drag were considered, parameterized by the Rayleigh drag timescale τ_drag: weak (τ_drag = 10⁷ s) and strong (τ_drag = 10⁴ s). Scenarios with and without hydrogen dissociation/recombination were modeled in order to assess its dynamical impact on the global heat transport. These effects were incorporated into the GCM by tracking the local atomic hydrogen mixing ratio and computing the recombination-driven heating/cooling that results when the atomic hydrogen mixing ratio decreases/increases, along with local changes in specific heat and mean molecular weight due to the spatially varying ratio of atomic to molecular hydrogen (Tan & Komacek 2019). For the GCM simulations at dayside-redistribution equilibrium temperatures below 4500 K, which are taken directly from Tan & Komacek (2019), we uniformly assumed a planetary radius of 1.47 R_Jup, a surface gravity of 11 m s⁻², and a rotation period of 2.43 days. Meanwhile, the GCM simulations of KELT-9b used a radius of 1.89 R_Jup, a surface gravity of 19.95 m s⁻², and a rotation period of 1.48 days (Mansfield et al. 2020).

Zoom In Zoom Out Reset image size

Crucially, Figure 17 shows that the majority of extant phase-curve observations of hot and ultrahot Jupiters cannot distinguish between GCM predictions with and without hydrogen dissociation/recombination. This is because the impacts of hydrogen dissociation/recombination on day–night heat transport are largest for gas giants with dayside-redistribution equilibrium temperatures in the range 3500 K ≲ T_eq ≲ 4500 K, where the hydrogen fraction varies most sharply with temperature and the spatial variations in the atomic hydrogen mixing ratio are large (Bell & Cowan 2018; Tan & Komacek 2019). At the extreme irradiation level of KELT-9b, the effects of hydrogen dissociation/recombination are largely "saturated," yielding a small day–night brightness temperature contrast that is consistent with observations (Mansfield et al. 2020; Wong et al. 2020e; Bell et al. 2021). As discussed in Section 5.1, TOI-2109b is the first ultrahot Jupiter discovered in the temperature range between KELT-9b and the cooler gas giants with T_eq < 3500 K, i.e., in the transitional zone where the differential effects of hydrogen dissociation/recombination on global heat transport are largest. Future spectroscopic phase-curve observations of TOI-2109b will empirically test our current understanding of the dynamical impacts of hydrogen dissociation/recombination on the global heat transport of ultrahot Jupiters.

TOI-2109b also presents an enticing case for studying atmospheric variability. A handful of previous works have reported modulations in the phase-curve amplitude and phase offset on ∼10 day timescales for a few exoplanets (e.g., Armstrong et al. 2016; Jackson et al. 2019; but see also Lally & Vanderburg 2020). The rapid rotation of TOI-2109b due to its ultrashort orbit can make its atmosphere particularly susceptible to hydrodynamic instabilities, which generate zonal propagation of waves at mid-to-high latitudes (Tan & Komacek 2019; Tan & Showman 2020). Other mechanisms that may yield orbit-to-orbit evolution in the global atmospheric properties include transient waves near the equator (Komacek & Showman 2020) and time-varying magnetohydrodynamic drag arising from the coupling of the planet's magnetic field with the partially ionized dayside atmosphere (Rogers 2017). The possibility of phase-curve variability can be explored when TESS returns to observe this system in the extended mission, and potentially in additional extended missions if approved by NASA.