TIC 172900988: A Transiting Circumbinary Planet Detected in One Sector of TESS Data

TESS observed TIC 172900988 (hereafter TIC 1729) during Sector 21 (2020 January 21 through 2020 February 18). Preliminary inspection of the light curve revealed one primary and two secondary eclipses (the first of which is only partially covered) with depths of ≈40% and ≈35%, respectively, and indicated an EB orbital period of ≈19.7 days with a significantly nonzero eccentricity. The light curve exhibits two additional transit-like events at days 1883.37 (BJD–2,457,000) and 1888.30 (just before and right after the downlink data gap), with depths of ≈3000 ppm (0.3%) and durations of ≈0.4 days and ≈0.27 days, respectively. Figure 1 shows the FFI eleanor light curve, both full scale to present the stellar eclipses and also zoomed in to highlight the CBP transits, as well as the diagnostic tests we used to rule out artifacts associated with the latter. The variations in the durations of the transits, combined with the lack of significant photocenter shift or systematic artifacts during the two events rule out false-positive scenarios due to a background EB or instrumental artifacts, and indicate a circumbinary body.

Zoom In Zoom Out Reset image size

The CBP transits, the stellar eclipses in TESS data, and the configuration of the system at the transit times are shown in more detail in Figure 2. The two transits are present in all four light-curve extractions from eleanor ("Raw," "Corrected," "PCA" ⁶⁹ , and "PSF" ⁷⁰ Flux) and are independent of the aperture used for the extraction, and the contamination from nearby sources is low ( ≈2 × 10⁻⁴). The closest and brightest known source inside the 13 × 13 target pixels, TIC 172900986, has an angular separation of ≈37'' and is ≈8 mag fainter (see Figure 3), resulting in a maximum flux contamination of ≈500 ppm—much smaller than the measured transit depths. The latter, combined with the catalog stellar radius (2.02R_⊙ from the TIC, 2.1R_⊙ from Gaia ⁷¹ ), indicates that the size of the transiting body is ≈10R_⊕, further strengthening the CBP hypothesis. We note that there is a contact binary just at the edge of the default 13 × 13 target pixel array for TIC 1729 (TIC 172900982, near Column 101, Row 1185 on Figure 3) but it is not a source of contamination as it is separated from TIC 1729 by almost 10 pixels and does not overlap the PSF, and the timescale of the variability (P = 0.23 days) is unrelated to the CBP transits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a relatively bright target (V = 10.141 ± 0.006 mag), there is an extensive archive of photometric data on TIC 1729 from multiple sources and with long baselines. The target was observed by ASAS-SN (Shappee et al. 2014; Kochanek et al. 2017), Evryscope (Law et al. 2015; Ratzloff et al. 2019), KELT (Pepper et al. 2007), and SuperWASP (Pollacco et al. 2006), and the stellar eclipses are detected in all data sets (see Figure 4). Analysis of the archival data confirmed the orbital period of the EB estimated from TESS data and also showed a clear change in the phase of the secondary eclipse relative to the primary eclipse (see Figure 5). This phase change is indicative of apsidal precession where the argument of periastron of the binary changes with time, and constrains the mass of the CBP.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

SuperWASP-North was an array of eight cameras located at the Roque de los Muchachos Observatory on La Palma. Each 200 mm/f1.8 lens was backed by 2k × 2k CCDs operating with a wide-band, 400–700 nm filter (Pollacco et al. 2006). TIC 172900988 was observed in 2006 July, accumulating 6000 photometric data points, then again in 2008, obtaining another 2000 points, and again from 2011 December to 2014, accumulating 12,000 more data points. TIC 172900988 is the only bright star in the 48'' extraction aperture.

The Kilodegree Extremely Little Telescope (KELT, Pepper et al. 2007, 2012) is a ground-based, small-aperture, wide-field photometric survey searching for transiting exoplanets, using a nonstandard passband equivalent to a broad V + R + I filter. It has a pixel scale of 23'' pixel⁻¹, and the photometric reduction procedures are described in Siverd et al. (2012). The KELT-North telescope at Winer Observatory in Arizona observed TIC 1729 from 2006 October 27 to 2012 April 22, obtaining 8059 observations.

We used g'-band observations from Evryscope-North, located at Mount Laguna Observatory, approximately 44 miles (70 km) southeast of Palomar Observatory in southern California. The Evryscope pair of robotic all-sky telescopes, sited in Chile and California, are arrays of 6 cm telescopes designed to cover the entire accessible sky from each observing site in every exposure. Together, the Evryscopes survey an instantaneous field of view (FOV) of 16,512 square degrees at 2 minutes cadence with 132 pixel⁻¹ resolution. The Evryscope system design and survey strategy are described in Law et al. (2015) and Ratzloff et al. (2019). The data were filtered to exclude observations with signal-to-noise ratios less than 40. A total of 93 points were rejected (1.56% of the observations). Eight additional points near HJD 2458909.84 were identified as outliers and omitted; unfortunately these were during a secondary eclipse. A total of 5656 points remained, and the HJD dates were converted to BJD_TBD using Jason Eastman's website tool ⁷² (Eastman et al. 2010). The pipeline magnitudes were converted to normalized fluxes by dividing by the out-of-eclipse median, and the uncertainties in fluxes were obtained from the geometric mean of the upper and lower uncertainties in magnitudes. The error bars were then scaled by a factor of 2.53 so that the uncertainties matched the rms scatter in the out-of-eclipse portions of the light curve.

2.3.1. Photometry

2.3.2. High Angular Resolution Imaging

As part of our standard process for validating transiting exoplanets to assess the effect of any contamination by bound or unbound companions on the derived planetary radii (Ciardi et al.2015), we observed TIC 1729 with infrared high-resolution adaptive optics (AO) imaging at Palomar Observatory. The Palomar Observatory observations were made with the PHARO instrument (Hayward et al. 2001) behind the natural guide star AO system P3K (Dekany et al. 2013) on 2021 February 24 UT in a standard five-point quincunx dither pattern with steps of 5''. Each dither position was observed three times, offset in position from each other by 05 for a total of 15 frames.

The camera was in the narrow-angle mode with a full field of view of ∼25'' and a pixel scale of approximately 0025 per pixel. Observations were made in the narrowband Br γ filter (λ_o = 2.1686; Δλ = 0.0326 μm) with an integration time of 5.6 s per frame (118 s total).

The AO data were processed and analyzed with a custom set of IDL tools. The science frames were flat-fielded and sky-subtracted. The flat fields were generated from a median average of dark-subtracted flats taken on sky. The flats were normalized such that the median value of the flats is unity. The sky frames were generated from the median average of the 15 dithered science frames; each science image was then sky-subtracted and flat-fielded. The reduced science frames were combined into a single combined image using an intrapixel interpolation that conserves flux, shifts the individual dithered frames by the appropriate fractional pixels, and median-coadds the frames. The final resolution of the combined dither was determined from the FWHM of the point-spread function to be 0097 (Figure 6).

Zoom In Zoom Out Reset image size

The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 20° at separations of integer multiples of the central source's FWHM (Furlan et al. 2017). The brightness of each injected source was scaled until standard aperture photometry detected it with 5σ significance. The resulting brightness of the injected sources relative to the target set the contrast limits at that injection location. The final 5σ limit at each separation was determined from the average of all of the determined limits at that separation and the uncertainty on the limit was set by the rms dispersion of the azimuthal slices at a given radial distance (Figure 6).

A single faint companion was detected 159 ± 001 to the southwest at a position angle of 237° ± 1°. The companion star is ΔKs = 7.22 ± 0.07 mag fainter than the primary star, yielding an apparent magnitude of K₂ = 16.0 ± 0.1 mag. Because the companion was only detected in the final reduced data, observations in a second filter were not obtained for a determination of the color temperature of the star. Further, the companion star is not detected in Gaia DR3. As such, we are unable to determine if the companion is bound to TIC 172900988 or if it is a background/foreground star. Based upon the analysis of the inner binary star, the two stars are late-F and early-G main-sequence stars (see Section 4.1 below). Thus if the detected companion is bound, it has a relative K-band brightness consistent with the star being a late-M dwarf (∼M8V). If it is indeed an M8V, the TESS magnitude to infrared color should be Tmag −Kmag ≈ 3.5 mag, making the real apparent TESS magnitude 10.5 mag fainter than the combined brightness of the primary stars—far too faint to be responsible for the observed planetary transits of ≈3000 ppm. If the companion star is bound to the binary and the distance of the system is ≈240 pc (Gaia DR3 and see below), the M8V star (M ≈ 0.08M_⊙) has a projected separation of 380 au, which would correspond to an orbital period of ≈5000 yr. If the companion star is not bound and has a typical color Tmag − Kmag ≈ 2 mag, it is still too faint to be responsible for the observed transits.

2.3.3. Spectroscopy

We obtained four spectra of TIC 1729 using the Tull Coudé spectrograph (Tull et al. 1995) on the Harlan J. Smith 2.7 m telescope (HJST) at McDonald Observatory in West Texas. Three spectra with net exposure times of 1200 s each were obtained 2020 May 13 and one additional spectrum with an exposure time of 1200 s was obtained 2020 May 14. The total wavelength coverage was 3760–10000 Å with a resolving power of R = 60,000. An iodine (I₂) cell was used, and I₂ absorption lines are superimposed on the stellar absorption lines between about 4992 Å and 6200 Å. The spectra were reduced using custom scripts based on standard IRAF ⁷⁴ routines. The reduced spectra have signal-to-noise ratios of between about 80 at 4400 Å to 120 per pixel at 6500 Å.

Four spectra of TIC 1729 were obtained using the Tillinghast Reflector Echelle Spectrograph (TRES, Szentgyorgyi & Furész 2007; Buchave et al. 2010) on the 1.5m telescope at the Fred L. Whipple Observatory in southern Arizona in October and December of 2020. These spectra cover a wavelength range from about 3850 to 9100 Å with a resolving power of R ≈ 44,000, and have signal-to-noise ratios of 30 to 35 per resolution element of 6.8 km s⁻¹ near 5200 Å. The spectra were reduced using the standard TRES data reduction pipeline described by Buchhave et al. (2010).

We obtained spectra of TIC 1729 on 2020 June 3 and 9 (one spectrum from each night) using the ARC Echelle Spectrograph (ARCES; Wang et al. 2003) on the Astrophysical Research Consortium 3.5 m telescope at Apache Point Observatory (APO) in New Mexico. ARCES has a resolving power of R ≈ 32,500 and a wide spectral coverage of 3200–10000 Å. Each night we obtained three 5 minutes exposures just after sunset at relatively high airmass (≈2.5–3). Quartz lamps were taken at the beginning of each half night, and thorium-argon lamps were taken at the beginning and end of each half night to perform wavelength calibration. The images were bias- and dark-subtracted using the Python package ccdproc (Craig et al. 2017). The spectral extraction, initial wavelength calibration, and normalization were performed using custom Python routines. The initial wavelength solution was performed by fitting a fourth-order polynomial individually for each order and then refined with a two-dimensional dispersion solution using the echelle package of the IRAF software.

Finally, eight spectra of the target were also obtained on the high-resolution spectrograph SOPHIE (Perruchot et al. 2008), which is mounted on the 193 cm telescope at Observatoire de Haute Provence (OHP), in France, and was designed and constructed to detect exoplanets (e.g., Bouchy et al. 2009). The instrument has a long-term stability of 2 m s⁻¹. The spectra were all taken with an exposure time of 900 s between 2020 November 15 and December 6 to cover several phases of the binary period.

There are no indications of additional stars in any of the obtained spectra. To extract the radial velocities of each star from the spectra we used the "broadening function" (BF) technique (Rucinski 1992, 2002). A high signal-to-noise ratio spectrum of a slowly rotating star is needed as a template for the BF analysis. For the McDonald spectra, high signal-to-noise ratio spectra (≈1000 per pixel) of diffused daytime sky near the zenith are routinely obtained, and we used one such spectrum obtained on 2020 May 12 as the template. This spectrum was obtained with the "Solar Port," which is a ground quartz diffuser on the roof of the coudé slit room. This sends the light falling on it into the spectrograph. That light is a combination of direct sunlight (the dominant component) and general daytime skylight. A spectrum of HD 182488 taken 2018 May 31 was used as the template for the TRES spectrum. A spectrum of HD 185144 taken 2021 June 13 was used as the template for the Sophie spectra. Finally, a spectrum of the twilight sky was used as the template for the ARCES spectra. The BFs were extracted separately for each order (in preparation for this analysis, the spectra were normalized to their local continua using spline fitting as implemented in IRAF). Representative BFs are shown for each instrument in Figures 7–9. With the exception of two Sophie observations (those from 2020 October 31 and November 11), the BF peaks from the primary and secondary stars are cleanly resolved, and no other peaks from additional stars (either bound to the system or along the line of sight) are evident. The BFs from the ARCES spectra (Figure 9) show a third component near zero velocity from night sky background.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Stellar Parameters for TIC 172900988

Parameter	Value	Error	Source
Identifying Information
TIC ID	172900988		TIC
Gaia ID	709709093301608448		Gaia EDR3
Tycho Reference ID	2483-160-1		TRC
α (hh:mm:ss)	08:34:38.799192	0.000013	Gaia EDR3
δ (dd:mm:ss)	+31:33:14.314250	0.00010	Gaia EDR3
μα (mas yr−1)	−6.5329	0.0145	Gaia EDR3
μδ (mas yr−1)	−22.8868	0.0115	Gaia EDR3
π (mas)	4.0607	0.0149	Gaia EDR3
Distance (pc)	246.263	0.904	Gaia EDR3
Photometric Properties
V (mag)	10.141	0.006	TIC
Gaia (mag)	10.0477	0.0006	Gaia EDR3
T (mag)	9.6319	0.0078	TIC
Contamination	0.000453		TIC
Stellar Properties
Teff (K)	6030	100	TIC
Teff, 1 (K)	6050	100	This work
Teff, 2 (K)	5983	100	This work
[Fe/H]	+0.34	0.10	This work
Age (Gyr)	3.1	0.1	This work

Download table as:  ASCIITypeset image

To obtain the final velocities for each spectrum, the following procedures were used. For the TRES spectra, we used 26 orders covering roughly 4500–6700 Å. The BFs were fitted with a double Gaussian model to find the velocity shifts and widths. There were no apparent trends with central wavelength (Figure 10). A horizontal line was fitted to the velocities and the error bars on the individual measurements were scaled to give ${\chi }_{\min }^{2}={N}_{\mathrm{dof}}$ for this fit. After the errors were scaled, a weighted average was computed to get the velocity for the primary and the velocity for the secondary (Table 4). A similar procedure was used for the Sophie spectra, where we used 39 orders covering roughly 3950–6830 Å. The adopted radial velocities for the six double-lined spectra are reported in Table 4.

Zoom In Zoom Out Reset image size

Table 4. Radial Velocities for TIC 172900988

Time	Date	Primary	Secondary	Instrument
(BJD–2,455,000)		Velocity (km s−1)	Velocity (km s−1)
3982.63491	2020-May-13T03:14:15.8	−9.248 ± 0.060	61.677 ± 0.067	Tull
3982.65341	2020-May-13T03:40:54.9	−9.114 ± 0.064	61.534 ± 0.043	Tull
3982.67077	2020-May-13T04:05:54.7	−8.995 ± 0.057	61.460 ± 0.064	Tull
3983.63222	2020-May-14T03:10:24.2	−2.758 ± 0.056	54.759 ± 0.054	Tull
4003.61687	2020-Jun-03T02:48:17.6	0.000 ± 0.027	53.057 ± 0.027 a	ARCES
4009.62264	2020-Jun-09T02:56:35.8	0.000 ± 0.027	−44.992 ± 0.028 a	ARCES
4131.99887	2020-Oct-09T11:58:22.1	93.471 ± 0.111	−43.072 ± 0.066	TRES
4135.02023	2020-Oct-12T12:29:08.2	−9.340 ± 0.098	62.841 ± 0.051	TRES
4138.65484	2020-Oct-16T03:42:58.3	−16.383 ± 0.071	69.606 ± 0.077	Sophie
4147.67188	2020-Oct-25T04:07:30.9	52.738 ± 0.022	−1.369 ± 0.026	Sophie
4147.67188	2020-Oct-25T04:07:30.9	52.738 ± 0.022	−1.369 ± 0.026	Sophie
4153.67808	2020-Oct-31T04:16:26.4	⋯	⋯	Sophie
4164.70059	2020-Nov-11T04:48:50.7	⋯	⋯	Sophie
4171.63576	2020-Nov-18T03:15:29.8	91.962 ± 0.019	−41.769 ± 0.021	Sophie
4175.69111	2020-Nov-22T04:35:11.8	−23.475 ± 0.025	77.104 ± 0.027	Sophie
4181.59568	2020-Nov-28T02:17:47.0	7.353 ± 0.022	45.441 ± 0.024	Sophie
4190.69055	2020-Dec-07T04:34:23.7	92.689 ± 0.017	−42.547 ± 0.019	Sophie
4191.97004	2020-Dec-08T11:16:51.5	77.852 ± 0.096	−27.075 ± 0.061	TRES
4194.99147	2020-Dec-11T11:47:42.7	−22.177 ± 0.099	76.021 ± 0.063	TRES

Note.

^a Velocity difference between primary and secondary.

Download table as:  ASCIITypeset image

For the McDonald spectra, we used 22 orders covering roughly 3900–6700 Å. Orders 42 and 43 were skipped owing to a "picket fence" artifact caused by internal reflections, and orders 23–37 were skipped owing to the iodine lines (while the pipeline for extracting precise radial velocities for single-lined spectra taken with the I₂ cell is well developed, a similar pipeline for double-lined spectra is not quite perfected). The BFs were fitted with double Gaussian models to find the velocity shifts and widths. There was an odd trend with wavelength in that the six bluest orders had velocities that were shifted to slightly larger values (Figure 10). This systematic shift seemed to be stronger in velocities for the secondary. Given this systematic shift, we skipped the first six bluest orders when computing the final velocities. As before, a horizontal line was fitted to the remaining velocities, and the individual error bars were scaled to give ${\chi }_{\min }^{2}={N}_{\mathrm{dof}}$. After the errors were scaled, a weighted average was computed to get the velocity for the primary and the velocity for the secondary (Table 4).

The ARCES spectra were the most challenging to extract radial velocities from. There were 51 orders used, covering from about 4000 Å to about 6700 Å (a few of the orders with strong telluric absorption lines were skipped). After the BFs were found, it was apparent that there was residual sky contamination, and that the sky contamination was stronger in the bluer orders. Given the extra peak, a triple Gaussian model was used to fit the BFs (Figure 9). There was a clear and strong trend seen in the velocities from individual orders (Figure 11). Although the wavelengths of the sky peaks seemed to be roughly constant with the echelle order, the velocities for both the primary and the secondary showed an almost linear trend with the central wavelength of the echelle order. This trend was most likely a result of atmospheric refraction at high airmass as the slit width (16) was larger than the seeing-limited image of the star (≈1'') on the given nights, resulting in a chromatic shift in the stellar line-spread function not seen in the residual sky spectrum that evenly filled the slit. The velocities at the reddest orders were nearly 2 km s⁻¹ larger than the velocities at the bluest orders. However, the velocity difference between the two peaks stayed reasonably constant over wavelength. Thus, for each night, we measured the separation between the two peaks in a manner similar to what was done for the TRES spectra and used the velocity differences in the photodynamical fitting described below.

Zoom In Zoom Out Reset image size

The flux fraction, which we define as the secondary flux divided by the primary flux, is a useful constraint on the photodynamical model. The BFs in each echelle order were modeled either with a double Gaussian model (McDonald, Sophie, and TRES) or a triple Gaussian model (ARCES), and the ratio of the areas under each peak was taken to be the flux fraction for that echelle order. Figure 12 shows the individual flux fractions as a function of the central wavelengths. The two stars are pretty similar, and there are no apparent trends in the flux fractions over these wavelengths (≈3900 to ≈6800 Å). The average flux fractions are 0.8687 ± 0.0052 for TRES, 0.8555 ± 0.0018 for the McDonald spectra, 0.8632 ± 0.0043 for Sophie, and 0.8544 ± 0.0079 for ARCES. To adopt a flux fraction for the photodynamical model, we combined the flux fractions from the TRES and Sophie instruments and computed a weighted average of 0.8669 ± 0.0037. The McDonald flux fractions were excluded since the iodine-free echelle orders have very little overlap with the filters used for the photometry (the B band with a peak near ≈4650 Å and filters redward of this). The ARCES flux fractions were excluded owing the contamination from night sky lines. Since the stars are very similar, we adopted a flux fraction of 0.8669 ± 0.0037 for all bandpasses in the photodynamical analysis described below.

Zoom In Zoom Out Reset image size

We have light curves in eight distinct bandpasses: the TESS bandpass, B, V, R, Sloan i (see Table 1), Sloan g (ASAS-SN and Evryscope), KELT, and WASP. A secondary eclipse was partially covered at the start of the TESS observations, and these data were not used in the fitting owing to difficulties with normalizing the light curve. We only want to fit photometric data near an eclipse or a transit event, and the selection of observations near these events from the archival data requires a good ephemeris. Fortunately it was not difficult to find a suitable ephemeris (P = 19.658025, T₀ = 3878.33860, where the time convention is BJD−2,455,000) using the TESS data, the RV data, and some of the early WASP data where partially covered primary eclipses are seen. Using this ephemeris, we selected data with phases within 0.02 of phase zero (e.g., the primary eclipse) and also with phases between 0.575 and 0.645 (e.g., the secondary eclipse) from the KELT, WASP, and Evryscope data sets. All of the photometric data were converted to flux units as necessary, and then normalized to have an out-of-eclipse flux of 1.0. In an effort to facilitate the fitting, we also fit for the measured times of the primary eclipse, the secondary eclipse, and the two transits in the TESS data. Experience has shown us that it is important to include the transit times in the fitting because the χ² penalty of the light-curve fit does not change once the predicted model transit does not overlap with the observed transit. When including the times, the likelihood of a model increases as the predicted transits get closer to the observed ones. We also have five "sets" of radial velocities: McDonald, TRES, Sophie, ARCES #1, and ARCES #2 (in the case of ARCES, we are fitting for the velocity difference between the two components). The systemic velocities for each set were fit separately as nuisance parameters, with primary and secondary velocities initially forced to have the same systemic velocity offset within each set.

For the model setup for TIC 1729, we have a total of 54 free parameters. There are parameters related to the orbit of the planet: the orbital period P₂, the time of barycentric (inferior) conjunction T_conj,2, the eccentricity parameters $\sqrt{{e}_{2}}\cos {\omega }_{2}$ and $\sqrt{{e}_{2}}\sin {\omega }_{2}$, the inclination angle i₂, and the nodal angle Ω₂. There are similar parameters related to the binary orbit: the orbital period P, the time of barycentric conjunction T_conj, the eccentricity parameters $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$, and the inclination angle i (the nodal angle of the binary is set to Ω = 0). The mass and radius of the planet are specified by the mass ratio Q₃ ≡ (M₁ + M₂)/M₃ (i.e., the ratio of the binary mass to the planet mass) and the radius ratio R₁/R₃ (primary radius to planet radius). The component masses and radii for the binary are specified by the primary mass M₁, the binary mass ratio Q ≡ M₂/M₁, the primary radius R₁, and the radius ratio R₁/R₂. The effective temperatures of the two stars are specified by the primary temperature T₁ and the temperature ratio T₂/T₁. We have light curves in eight bandpasses, so consequently we have 32 sets of limb-darkening coefficients with two per star per bandpass. The standard quadratic limb-darkening law given by

$$\begin{eqnarray*}&&I/{I}_{0}=1-{u}_{1}{(1-\mu )+{u}_{2}(1-\mu )}^{2}\end{eqnarray*}$$

was used (where $\mu =\cos \theta$ is the projected distance from the center of the stellar disk), but with the "triangular" sampling technique of Kipping (2013) with coefficients given by ${q}_{1}={({u}_{1}+{u}_{2})}^{2}$ and ${q}_{2}=0.5{u}_{1}{({u}_{1}+{u}_{2})}^{-1}$. The remaining free parameters are the apsidal constants for the primary (k_2,1) and secondary (k_2,2), and a contamination parameter to account for the other sources of light in the TESS aperture.

ELC has several fitting algorithms available. For TIC 1729, the algorithms we used were a simple Markov Chain Monte Carlo (MCMC) algorithm outlined in Tegmark et al. (2004), a "differential evolution" MCMC algorithm (DE-MCMC, Ter Braak 2006), and the nested sampling algorithm outlined by Skilling (2006). In the case of TIC 1729, the size of the parameter space to search is relatively large given that most of the orbital parameters of the planet are not immediately known. The planet's orbit is presumably close to edge-on given the two observed transits, but other than that the range of possibilities for the planet's orbital period, eccentricity, and argument of periastron are fairly large. Regarding the orbital period of the planet, we can compute the minimum orbital period based on dynamical stability, and we find ${P}_{3,\min }\approx 129.7$ days (Quarles et al. 2018). As discussed in Kostov et al. (2020b), the orbital period of the planet can be estimated from the timings and durations of the two transits and the detailed knowledge of the binary's orbit. Using this method, we estimated that the planet's orbital period is likely in the range between 150 and 340 days for an eccentricity smaller than 0.2 (details presented below).

Our basic strategy was to use the nested sampling algorithm to find possible solutions, and then to use the MCMC and the DE-MCMC algorithms to refine them. In the nested sampling algorithm, one defines the parameter space by giving minimum and maximum values of the free parameters. One then randomly selects "live points," where the number of them (N_live) should be at least twice the number of free parameters (N_parm). The likelihood values for each live point are computed. The "sampling" part of the nested sampling algorithm is relatively straightforward. Select a point in parameter space at random, and evaluate the likelihood. If that likelihood is better than the likelihood of the worst live point, then replace that live point with the new point. The discarded live point goes into the list of "dead points." If the new point has a worse likelihood than that of the worst live point, then that new point is discarded and not considered again. To make the sampling practical in a reasonably short period of time, the N_parm-dimensional hyper-ellipsoid that bounds the live points is computed, and new possible points are randomly selected from the hyper-ellipsoid (for large dimensions, the volume of the hyper-ellipsoid can be orders of magnitude smaller than the volume of the hyper-cube defined by the ranges of each fitting parameter). As better live points are found, the volume of the hyper-ellipsoid shrinks and potential replacement points are sampled from increasingly relevant regions of parameter space. Skilling (2006) showed that if the volume of the hyper-ellipsoid is successively scaled by a factor proportional to 1/N_live, then the live points and dead points can be used to compute the Bayesian evidence and also to construct a posterior sample. Unfortunately, when the number of parameters N_parm is larger than about 20, the statistically scaled volume of the hyper-ellipsoid is usually orders of magnitude larger than the actual volume of the hyper-ellipsoid. As a result, it can take a very long time (a few weeks of wall-clock time or longer for cases like TIC 1729) to run the nested sampling algorithm to completion, even when using a parallelized version that runs on multiple computer cores since the random samples usually come from regions of parameter space that are far from the region that has the largest likelihood. In our implementation of the nested sampling algorithm, we have the option of omitting the statistical scaling of the volume of the hyper-ellipsoid. The advantage of this is that it is much faster to find replacements for the live points since the volume from which the random samples are drawn is generally much smaller. Although we forfeit the ability to compute a posterior sample or to compute the evidence, the increased speed at which the algorithm can find good solutions as starting points for other algorithms makes up for this shortcoming.

Our quest to find good model solutions proceeded by iteration. When running the nested sampling algorithm, we kept the limb-darkening coefficients fixed, which left us with 22 free parameters and 44 live points. It is usually the case that when modeling circumbinary systems, reasonably good parameters for the binary can be found quickly, whereas good parameters for the planet (its orbital parameters and its mass) take longer to find. Once a good solution for the binary's parameters was found in the first run of the nested sampling, we used the MCMC and the DE-MCMC algorithms to refine the solution, which included the limb-darkening coefficients. With better values for the limb-darkening parameters, the nested sampling was run again with the limb-darkening coefficients fixed at better values. For the second iteration of the nested sampling, we used 66 live points. As the progress of the second run was monitored, it was found that not all values of the orbital period of the planet near the best value were possible. Instead, there are distinct families of solutions as shown in Figure 13. The figure shows the evolution of the live points as they are replaced with better models. We note that this is a representative figure based on the 2020 data. The "fingers" persist after the inclusion of the 2021 data (discussed below) but recreating this figure with all available data proved to be computationally prohibitive. However, we confirmed the pattern by using brute-force steppers in the planet period and found that regions of high χ² quickly emerge between the fingers, and the best-fit planet periods shift by about a day (compare location of fingers here to the periods listed in Table 6). At first (small y-values in the plot), the planet's period among the live points fills the range between about 180 and 210 days nearly uniformly. As better live points are found (moving up the y-axis), distinct "fingers" develop, and at the end, six families of solutions remain (top of the y-axis). This represents a problem for the nested sampling algorithm. The range of the planet periods is still fairly large, so the hyper-ellipsoid cannot become very small. Moreover, many, if not most, of the randomly selected points will fall in between the fingers where the likelihood is not very high. As a result, the convergence of the algorithm is slowed. In cases like this, the "multinest" routine discussed by Feroz et al. (2009), might offer relief, but we do not have an implementation of multinest within ELC.

Zoom In Zoom Out Reset image size

As an aside, we do not have a good explanation of why we see these families of solutions with distinct values of the planet's period, as opposed to a smooth distribution of points with a range of ≈10%–20% around the best value that one might expect based on the discussion in Kostov et al. (2020b). It is tempting to think that these discrete solutions are somehow analogous to some kind of aliasing caused by cycle count ambiguities over the relatively long time baseline we have (about 5600 days or about 15 yr). However, the fingers still appear when we only include the TESS data and the follow-up photometric data that were taken after the transits were identified. We note that the CBP Kepler-34 produced a "1–2 punch" where a pair of transits about 5 days apart were observed with Kepler (Welsh et al. 2012), but when we fit only those two transits (along with the stellar eclipses and radial-velocity data) using the nested sampling in a similar manner that was used on TIC 1729, we found only a broad range of possible planet periods with no distinct and separate families of solutions. Thus, it appears that the appearance of these distinct families of solutions is related to the peculiarities of TIC 1729 itself rather than being a generic feature of only having a pair of observed transits with a few days between them.

These distinct families of solutions with regions of low likelihood between them also present problems for the MCMC and the DE-MCMC algorithms. We therefore broke up the parameter space into six regions where the range of the planet's period was restricted in such a way to isolate each family, and leaving the ranges of the other 53 parameters to be the same (the identifying numbers for each family are shown in Figure 13). Using the best nested sampling solution from each family as an initial seed (see Orosz et al. 2019, for a discussion on how the chains are initialized), the DE-MCMC algorithm was run on each family. After these preliminary runs, we noticed two minor issues. First, the residuals for the secondary star radial velocities were systematically larger than the residuals for the primary star radial velocities by ≈75 m s⁻¹. Second, a comparison of the stellar masses and radii with stellar evolution models (see Section 4.1) showed a poor fit. We therefore made two modifications. First, with a physical interpretation in mind such as convective blueshift and gravitational redshift, we allowed the systemic radial velocities of the primary and secondary star to vary independently for each radial-velocity set. With this modification we can no longer use the ARCES velocities since we were fitting for the difference between the primary and secondary velocities. Second, we used the flux fraction as an observed constraint. As discussed previously, there is no apparent trend in the flux fractions with wavelength, so we adopted a flux fraction of 0.8669 ± 0.0037 for each bandpass.

After these modifications, the DE-MCMC code was run seven times for each family, with each run lasting 15,000 generations. After a generous burn-in period of 7500 generations, posterior samples were taken from every 1500th generation for each run and combined. For each posterior distribution (each of which had 8064 samples), the numbers are sorted and we computed the median value X_med, the value X_15.85 that marks the point where 15.85% of the distribution is smaller, and the value of X_84.15 that marks the point where 15.85% of the distribution is larger (i.e., 68.3% of the distribution is between X_15.85 and X_84.15). The adopted parameter value is taken to be the median X_med, and the 1σ uncertainty is taken to be the larger of (X_med − X_15.85) or (X_84.15 − X_med). The results for the fitted parameters are shown in Table 5, and for some derived parameters of interest in Table 6. From looking at Table 5, the solution for Family 5 (where P₃ = 200.452 days) has the smallest overall χ², with ${\chi }_{\min }^{2}$ = 11,348.96. The solutions for the other families are nearly as good with ${\chi }_{\min }^{2}$ values within ≈280. At this time it is difficult to choose between these these families with high confidence.

Table 5. Fitted Parameters for TIC 172900988 ^a

Parameter	Family 1	Family 2	Family 3	Family 4	Family 5	Family 6
${\chi }_{\min }^{2}$	11,410.26	11,192.97	11,174.33	11,164.23	11,141.53	11,187.00
Planet orbital parameters
P2 (days)	188.835 ± 0.092	190.375 ± 0.053	193.965 ± 0.078	198.856 ± 0.099	200.452 ± 0.011	204.046 ± 0.057
Tconj,2 b	3976.3 ± 2.2	3932.5 ± 1.1	3877.5 ± 2.8	3954.2 ± 3.1	3937.5 ± 0.5	3932.0 ± 0.9
$\sqrt{{e}_{2}}\cos {\omega }_{2}$	−0.2779 ± 0.0048	−0.2510 ± 0.0070	−0.2677 ± 0.0068	−0.2810 ± 0.0045	−0.1536 ± 0.0063	−0.1150 ± 0.0126
$\sqrt{{e}_{2}}\sin {\omega }_{2}$	−0.0734 ± 0.0091	−0.0058 ± 0.0176	0.0329 ± 0.0144	−0.0951 ± 0.0083	0.0578 ± 0.0149	0.1168 ± 0.0213
i2 (deg)	88.20 ± 0.21	88.18 ± 0.19	89.67 ± 0.25	91.84 ± 0.20	91.79 ± 0.19	90.50 ± 0.32
Ω2 (deg)	−0.88 ± 0.22	−0.74 ± 0.17	−2.18 ± 0.24	0.72 ± 0.17	0.32 ± 0.11	1.65 ± 0.15
Planet parameters
M3 (M⊕)	823.9 ± 4.8	845.7 ± 4.7	837.0 ± 4.7	869.7 ± 5.0	942.0 ± 5.6	981.3 ± 5.7
R1/R3	13.44 ± 0.54	13.54 ± 0.50	13.36 ± 0.47	13.49 ± 0.52	13.42 ± 0.53	13.40 ± 0.47
Binary orbital parameters
P1 (days)	19.657444	19.658083	19.658171	19.656652	19.658173	19.658016
	±0.000027	±0.000035	±0.000030	±0.000033	±0.000025	±0.000040
Tconj,1 b	3878.1746	3878.3572	3878.3823	3877.9473	3878.3835	3878.3386
	±0.0079	±0.0100	±0.0085	±0.0095	±0.0072	±0.0114
$\sqrt{{e}_{1}}\cos {\omega }_{1}$	0.23338 ± 0.00008	0.23326 ± 0.00008	0.23324 ± 0.00008	0.23354 ± 0.00008	0.23323 ± 0.00008	0.23326 ± 0.00008
$\sqrt{{e}_{1}}\sin {\omega }_{1}$	0.62716 ± 0.00017	0.62718 ± 0.00017	0.62729 ± 0.00015	0.62721 ± 0.00015	0.62730 ± 0.00016	0.62728 ± 0.00015
i1 (deg)	89.4182 ± 0.0049	89.4198 ± 0.0051	90.5940 ± 0.0046	90.5795 ± 0.0046	90.5762 ± 0.0047	89.4077 ± 0.0055
Ω1 (deg)	0 c	0 c	0 c	0 c	0 c	0 c
Stellar parameters
M1 (M⊙)	1.2388 ± 0.0008	1.2388 ± 0.0007	1.2388 ± 0.0007	1.2388 ± 0.0007	1.2384 ± 0.0007	1.2384 ± 0.0007
Q ≡ M2/M1	0.9706 ± 0.0003	0.9706 ± 0.0003	0.9706 ± 0.0003	0.9705 ± 0.0003	0.9705 ± 0.0003	0.9706 ± 0.0003
R1 (R⊙)	1.3846 ± 0.0019	1.3841 ± 0.0020	1.3842 ± 0.0017	1.3843 ± 0.0017	1.3842 ± 0.0016	1.3844 ± 0.0016
R1/R2	1.0527 ± 0.0017	1.0535 ± 0.0017	1.0539 ± 0.0016	1.0532 ± 0.0015	1.0547 ± 0.0015	1.0558 ± 0.0015
T2/T1	0.9903 ± 0.0029	0.9887 ± 0.0035	0.9882 ± 0.0022	0.9864 ± 0.0024	0.9866 ± 0.0030	0.9833 ± 0.0027

Notes.

^a Instantaneous parameters valid at BJD 2,453,250.0. ^b BJD–2,455,000. ^c Fixed.

Download table as:  ASCIITypeset image

Figures 2, 14, and 15 show the various photometric data and the best-fitting models (for Family 5). There do not appear to be any large systematic problems with the model fits. Comparing Figure 5, where the data are phased with a linear ephemeris, and Figure 14, where the times are relative to the nearest corresponding eclipse, we see that the apsidal motion is accounted for in the model. Finally, Figure 16 shows the radial velocities and the best-fitting model (for Family 5). The velocity residuals are all less than 200 meters per second. The residuals for the McDonald spectra are even tighter still, with a scatter of around 20 meters per second.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each of the three spectroscopic data sets the fitted center-of-mass velocities for the primary and secondary components are systematically different, and are always larger for the secondary. The differences (secondary minus primary) are +0.030 ±0.023 km s⁻¹ for McDonald, +0.082 ± 0.005 km s⁻¹ for TRES, and +0.076 ± 0.008 km s⁻¹ for SOPHIE. While in principle this could result from instrumental effects, the fact that the sign is the same for all three instruments, and the magnitudes quite similar, suggests an alternate explanation of an astrophysical nature. The component properties (mass, radius, and temperature) differ by less than 5%, yet this is enough to cause perceptible differences in the gravitational redshift, which we expect to be larger for the secondary by 0.010 km s⁻¹. Although convective blueshifts are not as easy to quantify, based on the prescription by Dravins et al. (1999) we expect the effect to be roughly 0.060 km s⁻¹ smaller for the secondary, resulting in a more positive velocity for that component. The two effects combine to give a larger expected velocity for the secondary by about 0.070 km s⁻¹, which appears to explain the bulk of the shifts mentioned above.

As noted earlier, the solution for Family 5 is formally the best. With one exception, there does not seem to be any pattern in how the various parameters change from family to family. That one exception is the mass of the planet: it is the smallest (M₃ = 823.9 ± 4.8 M_⊕) for Family 1 (where the planet's period P₃ is the shortest) and largest (M₃ = 981.3 ± 5.7 M_⊕) for Family 6 (where the planet's period P₃ is the largest). This trend is easy to understand. The planet is responsible for most of the apsidal motion seen in the binary. As the planet's orbital period gets longer, the planet's mass needs to go up to have the same effect on the binary.

Future transit and eclipse observations will allow us to determine which family represents the correct solution. We used the models from the posterior samples for each family to compute the times (and their uncertainties) of future transits (see Table 7) and eclipses (see Tables 8–19 in the Appendix). As one might expect, the predicted times of the next pair of observable transits depends on the orbital period of the planet. A clear detection of one or more transits in the future will unequivocally tell us which family represents the true solution. If future transit observations are not available, then observations of future eclipses will eventually allow us to discriminate between the various families. Figure 17 shows differences in the predicted eclipse times between Family 5 and Family 3, and between Family 5 and Family 6. The differences between the predicted times can be on the order of one to two minutes, which should be easy to detect.

Zoom In Zoom Out Reset image size

Finally, we note that there is a distinction between the initial conditions of the best-fitting models for each family and the parameters given in Table 5, which come from the medians of the posterior distributions. The parameters in Table 5 represent our best estimates of the parameters and their statistical uncertainties. However, those parameters collectively will not necessarily produce an optimal fit to the data owing to the complexity of the model and correlations among the various parameters. For the purposes of computing an optimal model, the initial conditions for the best-fitting models should be used. The initial conditions for the best-fitting models for all of the families are given in Tables 20–25 in the Appendix.

The photodynamical model predicted transits of the planet occurring during a conjunction in early 2021 (February–March, depending on the family). To detect the transits, we obtained four observations from CHEOPS via a generous allocation of Director's Discretionary Time. These data are shown in Figure 18. We note that these observations were scheduled to be centered on the predicted transit ingress corresponding to the best-fit model for the respective family based on the radial-velocity data available at the time, i.e., less than half of the RV measurements. As we received more measurements, the model fit to the binary star in particular improved and, as a result, some of the times of the predicted transits changed. This is the reason the current best-fit transit models for Family 1, 2, and 3 do not overlap with the CHEOPS observations.

Zoom In Zoom Out Reset image size

The egress part of an eclipse was observed near day 4264, and this observation includes a small part of the out-of-eclipse light curve. A fit to this partial eclipse was used to normalize all of the CHEOPS data. CHEOPS data were obtained near the expected transits for the Family 4 and Family 5 solutions. Unfortunately, there are two problems that make the interpretation of these data difficult. First, there are small but significant flux offsets between each visit as seen in Figure 18 (e.g., the observations covering out-of-eclipse phases in panels 1, 2, and 3 have fluxes that are different by a few tenths of a percent). These offsets might be produced by either a genuine stellar variability (which would be undetectable in TESS data), or could represent a systematic instrumental effect. Second, the duration of each visit was relatively short. As a result, the data near day 4266 could be consistent with being at the bottom of the Family 4 transit, or could be consistent with the out-of-transit light curve for all other families. Likewise, the predicted transit from the Family 5 solution partially overlaps with the egress of an eclipse. The light curve at the end of the corresponding CHEOPS visit is flat, but that could be consistent with being the out-of-eclipse light curve for Families 1 through 4 and Family 6, or it could be consistent with the bottom of the transit for the Family 5 solution.

We attempted to separately fit the CHEOPS data near the Family 4 and Family 5 transits. ELC is currently limited to eight distinct filter bandpasses, so we removed the data for the B filter and replaced the specific intensities in the model atmosphere table with intensities for CHEOPS using the appropriate filter response.

Thanks to a fantastic response from the community, we also obtained several ground-based observations of the predicted transits. In particular, for the potential Family 4 transit we obtained observations from KeplerCam in the Sloan i filter. These data are consistent with a flat line, as seen in Figure 19, where we also show the results of the fitting of the potential Family 4 solution using both the CHEOPS (black points) and KeplerCam (red points). We note that it is possible to orient the planet's orbit in such a way that no transit is observable near day 4265.5, so the model near the CHEOPS data and KeplerCam data is nearly flat. However, this occurs at the expense of the fit to the transit of the secondary star observed by TESS (middle panel in the figure). As a result, we can likely rule out the Family 4 solution, although the confidence level is not as high as one would normally like.

Zoom In Zoom Out Reset image size

For the potential Family 5 transit, we also have KeplerCam observations that are nearly simultaneous with the CHEOPS observations, so we fit these data in a similar manner as above. The results are shown in Figure 20. Here, the families where there is no transit blended with the egress of the eclipse fit the CHEOPS data better than does the Family 5 model (χ² ≈ 249 versus χ² = 290 for Family 5). However, the Family 5 model fits the KeplerCam data better than the other families (χ² = 4089 for Family 5 versus χ² ≈ 4381 for the others). Given the relatively short time span of the follow-up observations and the noise present in the KeplerCam data (especially near the end when morning twilight was approaching), Family 5 remains a viable solution.

Zoom In Zoom Out Reset image size

We also attempted to catch some of the predicted transits from TRAPPIST-North but the observations did not reach the necessary precision.

We discuss here the determination of the effective temperatures of the binary components as well as the system metallicity, based on the TRES spectra of TIC 1729 described earlier. We followed a cross-correlation procedure similar to the one described by Torres et al. (2002) to optimize the match of the observed spectrum against synthetic spectra, except that in this case the spectrum is double lined, so instead of 1D cross-correlations we used TODCOR (Zucker & Mazeh 1994), which is a two-dimensional correlation algorithm. This allows us to account for the presence of the secondary star in determining the stellar properties. Calculated templates were taken from a library of pre-computed spectra based on model atmospheres by R. L. Kurucz (see Nordström et al. 1994; Latham et al. 2002) covering the region centered on the Mg i b triplet at 5187 Å. We ran a grid of correlations solving for a mean effective temperature (strictly, a luminosity-weighted mean) and the rotational broadening of components, ${V}_{\mathrm{rot}}\sin i$, assumed to be the same for the two stars. Initial tests indicated the rotational broadening is below our spectral resolution, so we adopted ${V}_{\mathrm{rot}}\sin i=0$ km s⁻¹. As discussed by Torres et al. (2012), the temperatures determined in this way tend to be strongly correlated with the surface gravity and with metallicity, with higher values of these parameters leading to hotter temperatures. Given that in this case the surface gravities are well determined from the photodynamical analysis, we fixed them to a common value of $\mathrm{log}g=4.25$, near the average for the two components. On the assumption that the metallicity is solar, the resulting mean temperature we obtained is 〈T_eff〉 = 5770 K. Increasing the metallicity to [Fe/H] = +0.5 gave a hotter mean temperature of 〈T_eff〉 = 6140 K, as expected. Estimated uncertainties for these values are 100 K. Individual temperatures for the components were then computed by using the temperature ratio as measured from the photodynamical analysis.

A preliminary comparison against stellar evolution models from the PARSEC 1.2S series (Chen et al. 2014) indicated a satisfactory fit in the mass–radius plane for either metallicity (at different ages), but revealed that the predicted temperatures were hotter than we obtained when assuming solar metallicity, and cooler than we obtained when adopting the higher composition of [Fe/H] = +0.5. We then explored intermediate metallicity values, and found good agreement in both the mass–radius and mass–temperature diagrams for a composition of [Fe/H] = +0.34 ± 0.10. The mean system temperature at this composition is 6030 K, and the individual values are 6050 and 5983 K for the primary and secondary, respectively, corresponding to spectral types of approximately F9 and G0. This fit is shown in Figure 21, in which model isochrones are plotted for ages between 2 and 4 Gyr, in steps of 0.5 Gyr. The best simultaneous agreement with all observations is achieved for an age of about 3.1 Gyr, according to these models.

Zoom In Zoom Out Reset image size

As an independent test on the absolute temperatures, we gathered photometry for TIC 1729 from the literature and constructed 14 different but nonindependent color indices. We de-reddened them using an estimate of the interstellar reddening of E(B − V) = 0.03 from Schlafly & Finkbeiner (2011), and then applied 19 different empirical color−temperature calibrations from Casagrande et al. (2010) and Huang et al. (2015). Adopting [Fe/H] = +0.34 from the spectroscopic analysis above, we obtained a mean photometric temperature of 〈T_eff〉 = 5940 ± 100 K, marginally lower than the spectroscopic value, but within the uncertainties.

We ran an additional test on the spectroscopic parameters by performing a fit of the spectral energy distribution (SED) of TIC 1729 using the methods of Stassun & Torres (2016). Brightness measurements in the GALEX near-ultraviolet, Tycho-2 (B_T, V_T), Gaia (G, G_BP, G_RP), 2MASS (J, H, K_S), and WISE (W1–W4) passbands were gathered from the VizieR database, and cover the range from 0.2 to 22 μm (Figure 22). In addition, we pulled the GALEX far-ultraviolet flux to provide an estimate of the chromospheric activity. With the $\mathrm{log}g$ and effective temperature values as above, we obtained an excellent fit with ${\chi }_{\nu }^{2}=0.85$ (excluding the far-ultraviolet flux, which appears in excess) and a best-fit extinction A_V = 0.08 ± 0.03. The best-fit metallicity is [M/H] = +0.25 ± 0.05, strongly constrained by the near-ultraviolet flux in particular; this is somewhat lower but comparable to the estimate of +0.34 above. The distance inferred from comparing the integrated bolometric flux to the bolometric luminosity given by the radii and temperatures is 241.0 ± 3.6 pc, consistent with the distance of 245.0 ± 2.7 pc from Gaia EDR3, adjusted by the parallax offset recommended by Lindegren et al. (2021).

Zoom In Zoom Out Reset image size

The above results clearly point to a supersolar composition for TIC 1729, which is consistent with expectations for a massive planet and late F- to early G-type stars of the masses we determine (see, e.g., Santos et al. 2017).

To illustrate the evolutionary state of the stars, Figure 23 compares the surface gravity and temperature determinations against a different set of models from the MIST series (Choi et al. 2016). Evolutionary tracks are shown for the measured masses, and a metallicity of [Fe/H] = +0.30 ± 0.10 that best fits the observations, at which the individual spectroscopic temperatures are 6030 and 5963 K for the primary and secondary (6000 K for the luminosity-weighted mean). The best-fit age is 3.3 Gyr. These results are consistent with those obtained from the PARSEC models. Both stars are seen to be evolved, and are near the midpoint of their main-sequence lifetimes. For the final stellar parameters we adopt the average of the determinations based on the PARSEC and MIST models (T₁ = 6040 K, T₂ = 5970 K, [Fe/H] = +0.32), and report them in Table 3.

Zoom In Zoom Out Reset image size

We end this section with a note about the expected rotation of the stars. Given the 19.7 days orbital period, one might expect the stars to have been tidally spun-down and thus be rotating slowly. The timescale for reaching spin–orbit synchronicity is approximately 1.5 Gyr, significantly shorter than the age of the system (though the timescale for orbital circularization is much longer). However, the binary is quite eccentric (e = 0.45) and this means the equilibrium pseudosynchronous spin period (Hut 1981) is only about 8.3 days. The corresponding ${V}_{\mathrm{rot}}\sin i$ is 8.4 km s⁻¹ for the primary star. This somewhat relatively rapid spin may induce stellar activity that could be detected (e.g., starspot-induced photometric modulations, or the Mount Wilson Ca ii H and K S-index). Such enhanced activity, persistent on a long timescale, is reminiscent of the "forever young" effect (Mason et al. 2013), where the stars remain more rapidly rotating, and hence active, than their ages would imply.

The observed far-ultraviolet excess in the SED (Figure 22) can be used as a check on these ideas. For simplicity we assume that both of the (nearly identical) stars in the system contribute equally to the observed excess. We then infer a chromospheric activity of $\mathrm{log}R{{\prime} }_{\mathrm{HK}}=-4.8\pm 0.2$ via the empirical relations of Findeisen et al. (2011), which in turn predicts a rotation period of 10.8 ± 2.0 days via the empirical rotation−activity relations of Mamajek & Hillenbrand (2008). The latter relations also imply an age of 3.5 ± 1.4 Gyr, consistent with that determined above.

Based on the analysis of all available data, the observed phase change amounts to a change in the binary's argument of periastron of 2.88 × 10⁻³ degrees per cycle. Precession due to the effects of General Relativity accounts for a change of 1.70 × 10⁻⁴ degrees per cycle, and precession due to tidal bulges (assuming zero obliquity) accounts for 2.32 × 10⁻⁵ degrees per cycle, assuming apsidal constants of k₂ = 0.01 for both stars (Torres et al. 2010). Thus we see that most of the apsidal precession must be accounted for by another source, which is the circumbinary body of planetary mass.

As mentioned above, the measured times of the two CBP transits, combined with the RV-derived stellar masses and EB orbital parameters, allow an analytical estimate of the planet's orbital period (e.g., Kostov et al. 2020b) independent of the photodynamical analysis. Assuming e₃ = [0.0, 0.2] for the eccentricity of the planet—the range for the known transiting CBPs (Welsh & Orosz 2018)—and taking into account orbital stability, the calculated period of the CBP, P_3,calc, ranges from ≈152 days to ≈340 days, with an average of 209.8 days and a median of 194.1 days (see Figure 24). This is in line with the CBP orbital periods corresponding to the families of solutions listed in Table 6, and further strengthens their validity.

Zoom In Zoom Out Reset image size

Gaia measured an astrometric excess noise above a single-star model for TIC 1729 of 0.069 mas, with a significance of σ = 5.707 (Gaia EDR3, Gaia Collaboration et al. 2020). According to the catalog, there are 291 good observations (N = astrometric_n_good_obs_al) so the target falls in the large N limit and the excess noise is significant. At the measured parallax of 4.0608 mas, the astrometric excess noise for TIC 1729 corresponds to a measured physical displacement of Δa_measured = 0.017 au (Belokurov et al. 2020). We note that the binary has a high eccentricity (e = 0.45) and the two stars spend most of the time near apastron (${r}_{\max }=a(1+e)=0.27$ au). Accounting for the argument of periastron of the binary (696), the projected apastron is 0.25 au and thus the expected photocenter-induced physical displacement is Δa_expected = 0.01 au ⁷⁵ —comparable to the measured displacement. Additionally, Stassun & Torres (2021) have shown using a set of benchmark eclipsing binaries that the renormalized unit weight error (RUWE) statistic can be quantitatively used as an estimator of photocenter motion. Their empirical relation for the RUWE value for TIC 1729 (0.882) suggests a photocenter semimajor axis of ∼0.064 mas, or 0.016 au at the distance of the system. Thus there is evidence for binarity based on Gaia's astrometric excess noise as well as on RUWE.

Compared with single-star transiting exoplanets, the CBPs are required to have longer orbital periods and larger planetary radii. The longer period is due to the stability requirement around the binary, coupled with the lack of CBPs around the shortest-period binaries (Armstrong et al. 2014; Martin & Triaud 2015; Welsh et al. 2014). The larger radius is possibly a detection bias against finding smaller planets due to the CBPs' being intrinsically more difficult to detect than planets orbiting a single star (see, e.g., Windemuth et al. 2019; Martin & Fabrycky 2021). Yet the CBPs have smaller radii, on average, than the hot Jupiter subset. The observed radius distribution of the known CBPs may be an interesting combination of orbital migration history and observational bias against detecting planets that have not migrated or have been scattered to longer periods (see, e.g., Pierens & Nelson 2008; Periens & Nelson 2013; Kley & Haghighipour 2015; Kley et al. 2019; Penzlin et al. 2021). In general, the currently known CBPs experience low insolation, and a sizeable fraction of them are in the habitable zone of their host binaries, likely a fortuitous consequence of the bias of the Kepler sample toward G- and K-type stars. A comparison of the radii and insolations of the CBPs versus those of other transiting exoplanets is shown in Figure 25 (for CBP habitability see also Kane & Hinkel 2013). As seen in the figure, CBPs tend to reside in the larger-radius and lower-insolation portion of the diagram. However, it is difficult to know if the CBPs have a different parent population because of the different detection biases and small sample size. For comparison, the figure also marks planets that have two or more host stars but are not circumbinary (i.e., S-type planets). S-type planets in binary star systems are outlined with circles and those in triple or higher-order stellar systems are outlined with squares. Of the 2761 planets with radius and insolation estimates shown in the figure, 137 (5.0%) are in multiple star systems. Currently there are 125 (4.5%) noncircumbinary planets in multiple star systems, roughly a factor of 10 more than the circumbinary planets. Unlike the circumbinary planets, these S-type planets do not appear to favor any particular region in the diagram; they seem to follow the single-star distribution. Certainly this population does not show a tendency for experiencing low insolation or residing in the habitable zone. TIC 1729 b itself is too hot to be in the habitable zone, given its mean insolation of ∼5.2 S_⊕. Figure 26 shows a schematic birds-eye view of the TIC 1729 system, with the habitable zone region shown in green.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With solar-like stars in a moderately high eccentricity (e = 0.45) 19.7 days orbit, the binary hosting TIC 1729 b falls neatly into the distribution of known CBP systems, albeit with a higher than typical metallicity. The planet's orbit has a small mutual inclination and the ∼200 days period is typical of the known CBP population. The planet is close to the critical (in)stability radius: a/a_crit = 1.34 (P/P_crit = 1.55) using either the more sophisticated interpolation scheme in Quarles et al. (2018) ⁷⁶ or the approximation given in Holman & Wiegert (1999). In Figure 27 we show TIC 1729 b in comparison with the other known CBPs. With the notable exception of Kepler 1647 and KIC 7821010, the figure shows the accumulation of the orbits of the currently known CBPs near their boundaries of stability. This figure includes five candidate CBPs presented at various conferences going back as far as the Extreme Solar Systems II meeting. ⁷⁷ One of the candidates, KIC 10753734, is transiting (Orosz et al. 2016) and four others do not transit: KIC 7821010, KIC 8610483, KIC 3938073 (presented as a group for the first time at the "Toward Other Earths II" conference ⁷⁸ ), and KIC 5095269 (discovered by Getley et al. 2017). For consistency with the other 16 systems shown in Figure 27, we use our own preliminary photodynamical values for KIC 5095269 (J. Orosz et al. 2021, in preparation). The nontransiting planets are detected by their gravitational perturbation of the binary star's orbit, as manifested by eclipse timing variations (ETVs). While such a planet's radius is unknown, its mass is measurable by the pattern it produces in the ETVs (e.g., Borkovits et al. 2015). For the nontransiting cases, we used an approximate radius obtained by employing the mass–radius relation given in Bashi et al. (2017). While TIC 1729 b has the second-largest radius of the published circumbinary planets, exceeded only by Kepler 1647 b, its radius is in no way atypical. Where TIC 1729 b does differ from the published CBP is its mass: with a mass ≈2.9 M_Jup, it is the most massive transiting CBP known, a factor of 2 times larger than the next most massive planet, Kepler 1647b. However, both candidate CBPs KIC 7821010 and KIC 5095269 likely have comparable or larger masses.

Zoom In Zoom Out Reset image size

Models from the posterior samples for the six families were integrated using the GRK12 routine in ELC starting at day −1700 and going out to day 200,000 (i.e., from BJD 2,453,300 to BJD 2,655,000, which is about 530 yr into the future) in order to investigate near-term variations in the orbital elements and transit probabilities. For each model, the conjunction times and impact parameters for various combinations of the bodies over that time span are computed (e.g., primary eclipse, secondary eclipse, transits of the primary, etc.) and saved. In addition, the instantaneous values of the orbital elements and the Cartesian coordinates of all of the bodies are recorded every 10 days. A Lomb–Scargle periodogram was computed for the planet's orbital inclination time series to find the precession period of the planet's orbital plane. The results are given in Table 6. The precessional periods are between about 44 yr for Family 1 to about 54 yr for Family 6. Both the binary's inclination and the planet's orbital inclination vary with the same ≈50 yr period (depending on the family), but are out of phase. The binary's inclination varies by only about 05, and as such the binary will always be eclipsing.

However, the inclination of the planet's orbit varies by ≈5°, and this is a large enough change that transits will not always be possible (see also Schneider 1994, Kostov et al. 2014; Martin 2017). Figure 28 shows the impact parameter b of the (inferior) conjunctions of the planet with the primary and with the secondary over the span of about 100 yr (using the Family 5 solution). Only ≈40% of the conjunctions have ∣b∣ < 1, which is a necessary condition to have an observable transit. As seen from Table 6, the transit fractions ${\eta }_{\mathrm{prim}}$ and ${\eta }_{\sec }$ are about 15% and 40%. This transit fraction is consistent with the rough average of ≈36% for the first 10 known CBP systems (Martin 2017). Given that TESS observed TIC 1729 for only one sector out of thirteen in the Northern Hemisphere sky during its second year of operation, we were fortunate that TESS was able to observe the transits. Overall, the expected yield of TESS CBPs producing two (or more) transits during a single conjunction from Cycles 1 and 2 is on the order of a hundred (Kostov et al. 2020b).

Zoom In Zoom Out Reset image size

We performed 3-body integrations to understand both the short- and long-term evolution of Family 5 better. We used the ias15 integrator in REBOUND (Rein & Liu 2012) for the short-term (1000 yr) integrations, which shows the planetary semimajor axis and eccentricity vary in time (Figures 29(a) and (b)). These variations are bounded with a small amplitude, which is indicative of long-term stability (i.e., no secular growth). As a result the x-component of the eccentricity (Figure 29(c)) and inclination (29(d)) vectors also show a well-behaved evolution, where the small-scale variations in the eccentricity (Figure 29(c)) are due to the 3-body interaction with the inner binary. As mentioned in the discussion of "transitability," the orbital plane of the planet precesses over several decades. Thus, we extended the initial short-term integrations to 10,000 yr and identified the dominant secular frequencies (transformed to periods) using the fft module from scipy. Figures 29(e) and (f) illustrate the periodograms associated with the evolution of the planetary eccentricity and inclination vectors, respectively. The secular cycle for the planetary eccentricity using Family 5 is ≈49.3 yr, while the inclination is a little longer (≈51.8 yr). There is a longer variation in the inclination over ≈4000 yr (not shown in Figure 29(d)), where this envelope expands slightly in the range of inclination variation to ≈888–923 without changing the mean value.

Zoom In Zoom Out Reset image size

The long-term evolution of Family 5 used a different approach, where we investigated the stability of many initial conditions using 10⁵ yr integrations with the well-established mercury6 code that is modified to efficiently evolve planets in binary systems (Chambers et al. 2002). These simulations keep the binary solution for Family 5 fixed while varying the planetary semimajor axis and eccentricity in the planetary solution in the ranges 0.45–1.0 au (0.001 au steps) and 0.0–0.36 (0.002 steps), respectively. A simulation is terminated if the planet collides with the inner binary (i.e., pericenter crosses the inner binary orbits), escapes the system (i.e., distance to the center of mass is greater than 10 au), or the full simulation time (10⁵ yr) elapses. We use these simulations to measure the amplitude of eccentricity variation (${e}_{\max }-{e}_{\min }$) of a planet for a given initial condition in semimajor axis a_p and eccentricity e_p . This technique is similar to using the chaos indicator MEGNO (Cincotta & Simó 2000), but is less computationally expensive. Figure 30 illustrates the nominal location of the binary N: 1 mean motion resonances (MMRs) with gray curves, which typically shape the parameter space that stable CBPs can inhabit (Mardling 2013). The white cells denote initial conditions that undergo an instability, where the colored cells show the amplitude of eccentricity oscillations. The dashed curve marks the estimate for the stability boundary accounting for planetary eccentricity using Quarles et al. (2018). The green dot marks the nominal parameters for Family 5 (just outside the 10:1 MMR), where the parameters for other Families (gray dots) also lie in the stable regions around the 10:1 MMR. Note that solutions for interior to Family 5 have elevated eccentricity. All six Families represent stable 3-body solutions as long as the planetary eccentricity remains low (e_p ≲ 0.15) despite the forced eccentricity (∼0.02, yellow strip in Figure 30) from the inner binary.

Zoom In Zoom Out Reset image size