EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 702 (October 2025) A&A, 702 (2025) A85 Full HTML
Press Release
Open Access
Issue
A&A
Volume 702, October 2025
Article Number A85
Number of page(s) 30
Section Planets, planetary systems, and small bodies
DOI https://doi.org/10.1051/0004-6361/202554419
Published online 24 October 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

The most common type of star in the galaxy are M dwarfs, representing 60–75% of stars within 10 pc from the Solar System (see, e.g., Henry et al. 2006; Reylé et al. 2021). This prevalence, along with their small sizes, low masses, and cool temperatures, facilitates planetary searches and enables detailed characterisation of the exoplanets orbiting them, especially terrestrial ones. This situation is referred to as the M-dwarf opportunity (Charbonneau & Deming 2007; Triaud 2021), and in recent years it has placed M-dwarfs in the spotlight of many surveys and initiatives, such as SPECULOOS (Gillon 2018; Sebastian et al. 2021; Zúñiga-Fernández et al. 2025), CARMENES (Ribas et al. 2023), MEarth (Dittmann et al. 2016), TIERRAS (Garcia-Mejia et al. 2020), Project EDEN (Gibbs et al. 2020; Dietrich et al. 2023), and Red Dots1 (Anglada-Escudé et al. 2016; Dreizler et al. 2020), among others. In synergy with these projects, the Kepler/K2 missions (Borucki et al. 2010) and the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) have also significantly contributed to this endeavour, reaching the current number of ~170 small planets with sizes ≤1.5 R⊕ found orbiting M dwarfs2.

These efforts have provided a robust framework for advancing the detailed characterisation of Earth-like planets. Some examples are the Rocky Worlds DDT program used with the Hubble Space Telescope (HST) and JWST (Redfield et al. 2024), whose main objective is to determine the existence of atmospheres in these planets and establish limits for the cosmic shoreline (Zahnle & Catling 2017), and the TRAPPIST-1 JWST initiative (TRAPPIST-1 JWST Community Initiative 2024), which is drawing a road map to efficiently characterise the M8.5 star TRAPPIST-1 and its seven Earth-like planets. The aim of these new scientific initiatives is to elucidate the processes governing the formation and evolution of such planets, thereby laying the foundation for assessing their potential habitability with greater precision, a topic that is actively debated and largely speculated within the community (see, e.g., Dressing & Charbonneau 2015; Shields et al. 2016; O’Malley-James & Kaltenegger 2019; Childs et al. 2022).

Stars typically form as part of multiple stellar systems, with binary systems being the most common configuration (Offner et al. 2023). Despite this, stellar multiplicity has often been overlooked in statistical studies of planet occurrence rates, even though it plays a pivotal role in shaping the architectures and evolution of planetary systems (Marzari & Thebault 2019; Bonavita & Desidera 2020). The gravitational influence of a nearby stellar companion can significantly alter the physical conditions under which planets form by perturbing the protoplanetary environment. In particular, close companions can truncate circumstellar discs, reduce their masses and lifetimes, and induce misalignments, thereby influencing both the efficiency and the outcome of planet formation processes (see, e.g., Jang-Condell 2015).

Recent high-resolution observations from ALMA (Atacama Large Millimeter/submillimeter Array) and the VLT (Very Large Telecope) have markedly expanded our understanding of circumstellar disc evolution in binary systems. For instance, surveys of nearby star-forming regions (Manara et al. 2019; Akeson et al. 2019; Zurlo et al. 2020) have revealed that circumstellar discs are not only prevalent in many Class II binaries but can also be hosted by both the primary and secondary stars. Notable systems include IRAS04158+2805 (Ragusa et al. 2021), AS205, EM*SR24, and FUOri (Weber et al. 2023) as well as HBC494 (Nogueira et al. 2023). These observations support the theoretical expectation that binaries with separations of a few tens to several hundred AU commonly sustain circumstellar discs around each component (Manara et al. 2023). The diversity in disc morphologies and lifetimes across such systems underscores the complexity of planet formation in a multi-stellar context.

While stellar companions can disrupt disc evolution and trigger dynamical perturbations that influence planetary migration or even ejection (Kaib et al. 2013), numerical simulations have shown that stable planetary configurations can still exist, especially when planets reside within a few tenths of the binary separation (Holman & Wiegert 1999; Quarles et al. 2020). Observationally, this is supported by the discovery of over 500 exoplanets in multi-stellar systems, many of which orbit main-sequence binaries3. However, the presence of a stellar companion introduces additional challenges for detection and characterisation, particularly in transit surveys. Flux contamination from unresolved companions can dilute transit depths, leading to systematic underestimation of planet radii, masses, and bulk densities (Ciardi et al. 2015; Furlan & Howell 2020), and in some cases, it can hinder the detection of Earth-sized planets altogether (Lester et al. 2021). Moreover, since the habitable zone depends on the luminosity of the host star, identifying whether a transiting planet lies within the habitable zone becomes more complex in multi-stellar systems (Savel et al. 2020).

This study examines the nearby binary system TOI-2267 (see Table 1 for further designations), which consists of an M5V and a M6V star with a projected separation of approximately 8 au. The system hosts at least two warm transiting Earth-sized exoplanets with orbital periods near a 3:2 mean motion resonance. These findings identify TOI-2267 as the coolest binary system with the smallest stellar projected separation known to host planets. The planetary candidate TOI-2267.02 is also assessed, although its planetary nature cannot be fully confirmed at this stage. If future observations validate TOI-2267.02 as a planet, N-body simulations indicate that the three planets cannot orbit the same stellar component. Therefore, TOI-2267.02 would likely orbit a different star than the b-c planets, which would make TOI-2267 the first binary system identified with transiting planets around both stellar components.

This paper is structured as follows: in Section 2 we describe our stellar characterisation efforts, and in Section 3 we summarise all participating facilities and our observations. Next, in Section 4, we check for possible observational bias and false positive scenarios and present the statistical validation of the system. Section 5.1 encompasses our global analysis of the system and the resulting planetary system parameters. In Section 6 we describe the procedure for our independent planetary searches and how we established detection limits. In Section 7, we evaluate the system architecture, and finally, we discuss our findings and draw our conclusions in Sections 8 and 9, respectively.

Table 1

TOI-2267 stellar designations and astrometric properties.

2 Stellar characterisation

2.1 Spectral energy distribution

To determine the basic stellar parameters, we performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia DR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016); Stassun et al. (2017, 2018). We pulled the JHKS magnitudes from 2MASS, the W1–W4 magnitudes from WISE, the iy magnitudes from Pan-STARRS, and the GBPGRP magnitudes from Gaia. Together, the available photometry spans the full stellar SED over the wavelength range 0.4–20 μm (see Fig. 1).

Due to the presence of the close companion observed in the speckle imaging (see Sect. 4.3), which contributes flux in all of the observed combined-light broadband photometric measurements, we performed a two-component fit following Stassun & Torres (2016). To this end we used NextGen stellar atmosphere models, with the free parameters being the effective temperatures (Teff) and radii (R⋆), constrained by the Gaia distance and Ic-band flux ratio (FI) measured from the speckle imaging (see Sect. 4.3.2). We set the extinction AV ≡ 0 due to the proximity of the system to Earth. Simple integration of the best-fit SEDs yields the bolometric fluxes at Earth (Fbol).

The resulting fit (Fig. 1) has a reduced χ2 of 1.3 and best fit parameters of Fbol,1 = 2.35 ± 0.19 × 10−10 erg s−1 cm−2, Fbol,2 = 4.74 ± 1.83 × 10−11 erg s−1 cm−2, R⋆,1 = 0.2075 ± 0.0225 R⊙, R⋆,2 = 0.13 ± 0.03 R⊙, Teff,1 = 3030 ± 100 K, Teff,2 = 2930 ± 160 K, d = 22.54 ± 0.18 pc, and FI = 0.300 ± 0.048. We note that the nature of the fit and the available constraints makes the fitted parameters highly correlated. Hence, for completeness, we provide the full parameter correlation matrix (see Appendix B).

Using a Monte Carlo approach, we employed the solar metallicity NextGen stellar atmosphere models to calculate synthetic photometry with associated uncertainties for each stellar component. The Species toolkit (Stolker et al. 2020) was utilised to interpolate the NextGen grid for the stellar parameters Teff, R⋆, parallax, and log g, and to compute synthetic photometry using the filter profiles of the selected magnitudes. The parameter distributions for Teff and R⋆ were derived from the SED fitting described earlier, while the parallax was taken from Gaia DR3 system values, and log g values were obtained from Table 2. The filter profiles used in this analysis were sourced from the SVO Filter Profile Service (Rodrigo & Solano 2020) and included Gaia DR3, 2MASS, and WISE filters. As with the SED fitting, interstellar extinction was neglected due to the system’s proximity. For each stellar component, we performed 1000 random samples of each parameter assuming Gaussian distributions, interpolated the NextGen grid, and calculated synthetic magnitudes for the relevant filter profiles. Since the upper boundary of the model grid for log g is 5.5 dex, we modified the random sampling process to cap any values exceeding this limit at 5.5 dex. The mean and standard deviation of the resulting magnitude distributions were adopted as the synthetic magnitudes and their uncertainties for each filter profile. The computed synthetic magnitudes and uncertainties are presented in Table 3.

Table 2

Derived properties of the primary and secondary stars of the TOI-2267 system.

thumbnail Fig. 1

Spectral energy distribution of TOI-2267. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit NextGen stellar atmosphere model for the two stellar components (hot component in blue, cool component in red, combined light in black).

2.2 Spectroscopic analysis

To extract the radial velocities (RVs) and projected rotational velocities (v sin i) of TOI-2267, we obtained the NIRSPEC high-resolution near-infrared spectra (R~35 000) on the W. M. Keck II Telescope (McLean et al. 1998, 2000; Martin et al. 2018). The NIRSPEC spectra were taken on 2022 March 13 (UT) under good seeing (~0.6″) and ~2 mm of precipitable water vapour. We chose the “Kband-new” filter with the 0⋅′′432×12′′$\[0^{\prime\prime}_\cdot432{\times}12^{\prime\prime}\]$ slit, and set the echelle and cross-disperser angles of 62.95° and 35.69°, respectively. We then took the A-B nodding sequence for TOI-2267, each with 300 s exposure, with the median signal-to-noise ratios per pixel S/N~23 and 25, respectively. We note that the double traces were shown due to a focusing issue at a high declination/elevation, not due to the nature of the binarity of TOI-2267. We used the A0V star HD 25175 to serve as our wavelength solution calibrator. We also acquired the associated dark, flat lamp data for our data reduction. Our NIRSPEC data were reduced using a modified version of NIRSPEC Data Reduction Pipeline (Tran et al. 2016), with the changes detailed in Hsu et al. (2021); Theissen et al. (2022).

Our observed Keck/NIRSPEC spectra and best-fit models are shown in Fig. 2 (further details see Appendix C). The stellar parameters for both nodes are consistent with the model uncertainties. The best-fit parameters after combining the results of both spectra are RV = −19.0 ± 0.3 km s−1, v sin i = 17.6 ± 0.3 km s−1 (Table 2), effective temperature Teff = 3167−13+14$\[3167_{-13}^{+14}\]$ K, and surface gravity log g = 5.18−0.07+0.08$\[5.18_{-0.07}^{+0.08}\]$ cgs dex. We note that our Teff and log g are merely for reference, as the high-resolution spectra only focused on narrow spectral regions (Hsu et al. 2021, 2024), while the robust and accurate Teff and log g measurements are determined using the SED and Sect. 2.4 (Table 2). Finally, we also examined if the secondary component was detected as combined light spectra in our NIRSPEC data. We forward-modelled a binary template, following Hsu et al. (2023), but we did not find strong evidence of detection of the secondary.

Table 3

Photometric properties of TOI-2267.

thumbnail Fig. 2

Keck/NIRSPEC K-band spectra of TOI-2267 around 2.3 μm. Top: observed NIRSPEC spectrum in the nodding position A (shown with the black line) with the data noise in the grey-shaded region. The stellar models with and without telluric absorption features are indicated with the blue and green lines, respectively. The residual (data – model) is denoted with the magenta line. Bottom: same as the top panel for NIRSPEC spectra and model fit but in the nodding position B.
thumbnail Fig. 3

Periodogram obtained with TESSExtractor tool for TOI-2267 from TESS Sector 53.

2.3 Rotational periods from TESS photometry

We performed a Lomb–Scargle periodogram analysis on TESS light curves using the Lightkurve package (Lightkurve Collaboration et al. 2018). We detected three main peaks in the periodogram, with the period at maximum power being 0.6958 d, consistent with the 0.695 d rotational period published by Newton et al. (2016). Using this period and the Lomb–Scargle model method from Lightkurve, we detrended the light curve and then repeated the periodogram analysis. The period at maximum power, in this case, was 0.4936 d, which could correspond to the rotational period of the secondary star in our binary system. We repeated the detrending and periodogram procedure as described before and found that the highest peak corresponded to a period of 0.3471 days, which aligns with the second harmonic of the main rotational period. As an independent approach, we used the TESSExtractor tool4 (Brasseur et al. 2019; Serna et al. 2021), which performs a correction of systematic effects in the light curve and a periodogram analysis on TESS data for a given target and sector. The results that we obtained from TESSExtractor tool are in agreement with the results of our analysis (see Fig. 3).

2.4 Estimated masses and evolutionary modelling

We first relied on evolutionary modelling to obtain the stellar masses from the models of very low-mass stars and young brown dwarfs presented in Fernandes et al. (2019). We used as constraints the luminosity computed from the bolometric flux Fbol and distance d derived from SED fitting (see Sect. 2.1), the metallicity [Fe/H]= 0.164 ± 0.11 (Dittmann et al. 2016), and assuming an age of ≳1 Gyr. Since very low-mass stars evolve so slowly, they keep a constant luminosity once the star has turned on core H-burning and has reached the main sequence. This implies we cannot infer the stellar age from evolutionary modelling. Another consequence is that we could model the two stars separately without considering co-evolution through binary evolution modelling.

We obtained a stellar mass of M⋆,1 = 0.1710 ± 0.0079 M⊙ for the M5V star and M⋆,2 = 0.0989 ± 0.0130M⊙ for the M6V star. The quoted uncertainties reflect the error propagation on the stellar luminosity and metallicity, but also the uncertainty associated with the input physics of the stellar models (Van Grootel et al. 2018). The inferred radii from models are R⋆,1 = 0.204 ± 0.021 R⊙ and R⋆,2 = 0.125 ± 0.005 R⊙, within 1-σ agreement with the radii derived from SED fitting (Sect. 2.1). Considering these SED radius estimates and the stellar masses inferred from evolutionary modelling, we derived stellar surface gravities of log g⋆,1 = 4.99 ± 0.05 and log g⋆,2 = 5.28 ± 0.18 (cgs), and stellar mean densities of ρ⋆,1 = 23.0−4.0+5.3$\[23.0_{-4.0}^{+5.3}\]$ g cm−3 and ρ⋆,2 = 83−39+113$\[83_{-39}^{+113}\]$ g cm−3.

3 Photometric observations

This section details all of the observations and facilities used to study the TOI-2267 system. This effort corresponds to time-series photometry measured in nine TESS sectors plus observations gathered by four ground-based telescopes used for the photometric follow-up.

3.1 TESS

TESS observed TOI-2267 with a 2 min cadence during the primary mission in Sector 19 (midpoint December 2019), Sector 20 (midpoint January 2020), Sector 25 (midpoint May 2020), and Sector 26 (midpoint June 2020). The Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016) reduced the image data to produce the photometric time series. These time series were used to search for transiting planet signatures with a noise-compensating matched filter (Jenkins 2002; Jenkins et al. 2010, 2020), leading to the detection of a transit signature with an orbital period of ~3.5 d. The signal passed all diagnostic tests (Twicken et al. 2018; Li et al. 2012) and was released as the planetary candidate TOI-2267.01 on September 30, 2020 (Guerrero et al. 2021). No other signal was released at that time. We conducted an independent transit search using our public pipeline SHERLOCK (see Sect. 6), which allowed us not only to recover the signal corresponding to TOI-2267.01 but also to find two other low S/N signals with orbital periods of 2.29 d and 2.03 d. We then triggered a ground-based follow-up campaign of this system, trying to confirm or refute the three signals detected. We discuss this follow-up effort in the following subsections. During the TESS mission extensions, TOI-2267 was re-observed in several sectors. First, it was re-observed in Sector 40 (midpoint July 2021). The reanalysis of all the observations gathered until that date by SPOC drove to the release of the candidate TOI-2267.03 with an orbital period of 2.29 d on February 28, 2022. Then, TESS kept observing the star on Sector 52 (midpoint May 2022), Sector 53 (midpoint June 2022), Sector 59 (midpoint December 2022), and Sector 60 (midpoint January 2023). The reanalysis of these nine sectors together by SPOC led to the detection of a new transiting signal with an orbital period of ~2.03 d, which was released on August 15, 2023, as TOI-2267.02. The independent findings by SPOC of TOI-2267.03 and .02, matching our early detections, reinforced their credibility. Moreover, TESS observed TOI-2267 again in Sector 73 (midpoint December 2023), Sector 79 (midpoint June 2024), and Sector 86 (midpoint December 2024).

The field-of-view of these sectors is displayed in Fig. G.1, showing the aperture photometry used in each case with red squares and the stars inside it. In all cases, we observed that very few faint stars (Δmag > 5), one or two, are in the aperture photometry, at more than one TESS pixel away (>21 arcsec). Then, to perform our analyses, we retrieved the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP) fluxes, which are corrected from crowding and systematic effects (Stumpe et al. 2012, 2014; Smith et al. 2012). These light curves are displayed in Fig. F.1, with the location of the transits for TOI-2267.01, .02, and .03 candidates highlighted.

In addition, visual inspection reveals the existence of flares in the data, which, combined with the rapid rotation of both stars, reported in Sect. 2, suggest that TOI-2267 A and B are magnetically active. A comprehensive analysis of activity levels, including flare rates and energy distribution, will be addressed in a subsequent publication.

3.2 Ground based photometry

In the following subsubsections, we describe all the ground-based observations. The photometric observations logs and transit detection are also summarised in Table D.1.

3.2.1 SPECULOOS-North/Artemis and SAINT-EX

We observed a total of 11 transits of TOI-2267.03 and 12 of TOI-2267.01 with SPECULOOS-North/Artemis (SNO/Artemis; Burdanov et al. 2022) in the Sloan-z′ and -r′ filters with an exposure time of 10 and 55 s, respectively. In addition, we gathered 5 transits of TOI-2267.03 and 1 transit of TOI-2267.01 with SAINT-EX (Demory et al. 2020) in the I + z and Sloan-z′ filters with an exposure time of 10 s. SNO/Artemis and SAINTEX are 1.0-m Ritchey–Chrétien telescopes located at the Teide Observatory in the Canary Islands (Spain) and at the Observatorio Astronómico Nacional in San Pedro Mártir (Mexico), respectively. Both are equipped with a thermoelectrically cooled 2K×2K Andor iKon-L BEX2-DD CCD camera with a pixel scale of 0.35″, and a field-of-view (FoV) of 12′×12′. These facilities are twins of the SPECULOOS-South telescopes (Delrez et al. 2018; Zúñiga-Fernández et al. 2024). Data reduction, calibration, and photometry measurements were performed using the PROSE5 pipeline (Garcia et al. 2022).

3.2.2 TRAPPIST-North

Two transits of the outer planet TOI-2267.01 were observed with the TRAPPIST-North (Jehin et al. 2011; Gillon et al. 2011; Barkaoui et al. 2019) telescope. TRAPPIST-North is a 60-cm robotic telescope located at Oukaimeden Observatory in Morocco since 2016. It is equipped with a thermoelectrically cooled 2K×2K Andor iKon-L BEX2-DD CCD camera with a pixel scale of 0.6″, resulting in a fov of 20′×20′. Both transits were observed in the I + z′ filter with an exposure time of 30 s. The data was processed using the PROSE pipeline.

thumbnail Fig. 4

Archival images for TOI-2267. From left to right: 1955 DSS/POSS-I image, 1996 DSS/POSS-II image, 2011 PanSTARRS image, and 2025 Artemis-1m0 image. The red circle corresponds to the target position at 2022 November 09.

3.2.3 LCOGT

A full transit of the inner planet TOI-2267.03 was observed in the Sloan-i′ filter using Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0-m at McDonald Observatory. Complementary, three full transits of the outer planet TOI-2267.01 were also observed in the Sloan-i′ filter and one full transit of TOI-2267.02 with the z′ filter using LCOGT 1.0-m at Teide Observatory. The science images were calibrated using the standard LCOGT BANZAI pipeline (McCully et al. 2018), and photometric measurements were extracted using AstroImageJ6 (Collins et al. 2017) software.

3.2.4 Sierra Nevada Observatory-T150

Using the 1.52-m telescope at Sierra Nevada Observatory (OSN/T150), we carried out a full-transit observation of TOI-2267.03 and TOI-2267.01, and two full-transits of TOI-2267.02. We employed the Johnson-Cousin I filter in both observations with an exposure time of 90 s. The OSN/T150 is a Ritchey-Chrétien telescope equipped with a thermoelectrically cooled 2K×2K Andor iKon-L BEX2DD CCD camera with a fov of 7.9′ × 7.9′ and pixel scale of 0.232″. The photometric data were extracted using the AstroImageJ package.

3.2.5 Observatoire du Mont-Mégantic

We observed a full transit of TOI-2267.01 with the PESTO EMCCD camera (1024×1024 pixel) mounted on the 1.6-m telescope at Observatoire du Mont-Mégantic (OMM), Canada. PESTO has an image scale of 0.466″ per pixel, providing an onsky fov of 7.95′ × 7.95′. The observation sequence was taken in the i′ filter with exposures of 30 s. The data was processed using the AstroImageJ package.

3.2.6 Gran Telescopio Canarias

We observed one full transit of TOI-2267.02 with the 10.4-m Gran Telescopio Canarias (GTC), using the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS; Cepa et al. 2000), a tunable imager and spectrograph installed at the Nasmyth-B focus of the telescope. OSIRIS is equipped with a mosaic of two Marconi CCD42-82 detectors, providing an unvignetted field-of-view of 7.8′×7.8′ and a pixel scale of 0.127″ per pixel in 1×1 binning mode. The detectors are cooled using a continuous-flow cryostat and offer high quantum efficiency over the 365–1000 nm wavelength range. For our observations, we used OSIRIS in imaging mode with the Sloan-i filter and an exposure time of 1 s. Data reduction, calibration, and photometry measurements were performed using the AstroImage package.

4 Validation of the planets

4.1 Archival imaging

We used archival science images of TOI-2267 to exclude background stellar objects that could be blended with the target at its current position. Such an object might introduce the same transit event we observed in our data and skew the system’s physical p, which we obtained from the global analysis. TOI-2267 has a high proper motion of 280±0.44 mas/yr. We used images from DSS/POSS-I (Minkowski & Abell 1963) in 1955 and DSS/POSS-II 1996 in the red filter, PanSTARRS in 2011 in the y filter, and SNO/Artemis in 2025 in the z′ filter, spanning 70 years with our current observations. The target has moved by 19.6″ from 1995 to 2025. There is no source in the current day positions of TOI-2267; see Fig 4.

4.2 Seeing-limited photometry – SG1

One of the critical steps during the validation process of any TOI is the confirmation of the events in the target star at transit predicted times (see, e.g., Kostov et al. 2019; Günther et al. 2019). This strategy is motivated by the TESS’ large pixel size of 21 arcsec and the associated point-spread function (PSF) that could be as large as 1 arcmin. These two characteristics of TESS observations imply a larger probability of contamination by a nearby eclipsing binary (NEB). Indeed, due to dilution effects, deep eclipses in a faint NEB might mimic a shallow transit observed on the target star. Hence, it becomes critical to confirm the transits in the target star and, when those are very shallow for ground-based facilities, explore the potential contamination by relatively distant neighbours to rule out the contamination from an NEB.

On the one hand, for TOI-2267.01 and TOI-2267.03, we confirmed the transits at predicted times using ground-based facilities (see Table D.1 and Figs. 7 and 8), allowing us to rule out the presence of any NEB that might cause these transits and strengthening the planetary interpretation. On the other hand, for TOI-2267.02, we gathered three observations at transit times using the OSN/T150 and the LCO-Teide-1 m. Unfortunately, due to the shallowness of this transit of ~1.3 ppt, we could not confirm any event in the target star. We used these observations to search for potential NEBs up to 5′ from the target star in the OSN-1.5 m observations and up to 2.5′ in the LCO-Teide-1 m. None of these observations suggests that an NEB is the origin of the signal. In our effort to confirm this event in the target star, we got a Director’s Discretionary Time to use GTC with the Osiris camera (GTC/Osiris). The optimum instrument would have been HiPERCAM7 (Dhillon et al. 2021), which allows for multi-band simultaneous observations using u′, g′, r′, i′, and z′ bands; unfortunately, the instrument was not mounted when the planet transited. Still, while simultaneous multi-band observations cannot be conducted with GTC/Osiris, its photometric precision should be good enough to detect the transit in the target star. However, an inappropriate election for the exposure time remarkably decreased the photometric quality of the light curve, yielding a non-conclusive observation. We used this observation to search for potential NEBs again in the fov, and we found nothing.

Therefore, while we have not confirmed the candidate TOI-2267.02 in the target star, we cleaned the fov and ruled out the presence of any NEB that might induce the transit-like signal detected in the TESS data. Still, we prudently do not claim this candidate as a validated planet but as pending of a more robust follow-up, which is in our roadmap for this system.

4.3 High angular resolution imaging

Spatially close stellar companions can confound exoplanet discoveries such that the detected transit signal might be a false positive. In addition, even for real planet discoveries, the close companion will yield incorrect stellar and exoplanet parameters if unaccounted for (Ciardi et al. 2015; Furlan & Howell 2017; Furlan & Howell 2020).

4.3.1 Gemini-North-8.0 m/‘Alopeke

TOI-2267 was observed on 2021 December 09 UT using the ‘Alopeke speckle instrument on the Gemini North 8-m telescope (Scott et al. 2021). ‘Alopeke provides simultaneous speckle imaging in two bands (562 nm and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections. Eighteen sets of 1000×0.06 sec images were obtained and processed using our standard reduction pipeline (see Howell et al. 2011). Fig. 5 shows our final contrast curves and the 832 nm reconstructed speckle image. We detect a close companion to TOI-2267 residing at a position angle of 279.7 degrees and a separation of 0.384 arcsec. The companion star is 1.6 magnitudes fainter than the primary target as measured in the 832 nm filter. No additional close companions brighter than 4–5 magnitudes below that of the target star from the 8-m telescope diffraction limit (20 mas) out to 1.2″ were detected. At the distance of TOI-2267 (~22.55 pc), these angular limits correspond to spatial limits of ~0.45 to 27 AU.

4.3.2 SAI-2.5m

TOI-2267 was observed with the speckle polarimeter (Strakhov et al. 2023) on the 2.5-m telescope at the Caucasian Observatory of Sternberg Astronomical Institute (SAI) of Lomonosov Moscow State University on three dates, details are given in Table 4. For observations in 2020 and 2021, EMCCD Andor iXon 897 was used as a detector. The last observation was obtained with low-noise CMOS detector Hamamatsu ORCA-quest. During all observations, the atmospheric dispersion compensator was active, which allowed use of the Ic band. The respective angular resolution is 83 mas. Exposures and total number of accumulated frames were 60 ms and 3342, 60 ms and 6053, 23 ms and 5222, for three observational epochs, respectively.

In all epochs, a stellar companion was reliably detected (Table 4). Binary parameters were determined using an approximation of the average power spectrum by the method described in Strakhov et al. (2023). The 180° ambiguity in position angle characteristic for speckle interferometry was resolved by use of the bispectrum. While the flux ratio demonstrates consistency between observations, separation clearly decreased during the past 3.8 yr. It is interesting that at the same time, the position angle of the binary changed only a little, indicating almost edgeon orientation of the binary orbit. Further high angular resolution and radial velocity observations of the object in the coming years will allow us to constrain the orbit. The epoch with the best seeing conditions (2024 August 09) was used to put limits on additional companions: ΔIc = 2.9m and 4.2m at distances 0.25 and 1.0″ from the main star, respectively (see Fig. 6).

thumbnail Fig. 5

‘Alopeke speckle imaging 5σ contrast curves along with the reconstructed 832-nm image.
Table 4

Speckle interferometric observations of TOI-2267 at SAI-2.5 m.

4.4 Statistical validation

The final step of our validation procedure was statistically validating the three planetary candidates. To this end, we used the TRICERATOPS package (Giacalone et al. 2021), which is widely used to validate TESS planets statistically. TRICERATOPS uses the folded light curve of the candidates at their given epochs and periods and runs a Bayesian fit of several different possible astrophysical scenarios (see Table 1 and Sect. 2.2 from Giacalone et al. 2021), obtaining for each one the corresponding marginal likelihoods. These likelihoods are normalised once all the scenarios are examined so their sum equals one. Under this assumption, TRICERATOPS computes two parameters to elucidate if a given candidate is a validated planet, likely a planet, or a nearby false positive. These parameters are the False Positive Probability (FPP) and the Nearby False Positive Probability (NFPP) with the following thresholds for each case: (1) FPP <0.015 and NFPP <0.001 for validated planets, (2) FPP <0.5 and NFPP <0.001 for likely planetary signals, and (3) NFPP >0.1 for likely nearby false positives.

In our case, using the high-resolution images described in Sect. 4.3 we learned that the system is composed of two stars, both carefully characterised in Sect. 2. Since this information was not in GAIA, the default information imported by TRICERATOPS was incomplete. Fortunately, the package allows for updating the stellar information, adding new stars when needed, and appropriately considering them when computing the marginal likelihoods for each astrophysical scenario. Hence, using our results presented in Tables 2 and 3, we included in our validation process the two stars that compound the TOI-2267 system.

Moreover, TRICERATOPS allows for adding extra constraints depending on the follow-up observations gathered. In particular, it allows for dropping those astrophysical scenarios that are known to not be at work; for example, the scenarios TP, EB, and EBx2P correspond to the case where the target star does not have an unresolved companion; what we know is false in our case. Then, we dropped them. Moreover, by examining the archival images presented in Sect. 4.1, we know that there are no stars in the background at the current location of TOI-2267 that might be inducing the transit signals; hence, these scenarios can also be dropped; they are: DTP, DEB, DEBx2P, BTP, BEB, BEBx2P.

It is important to note that in this study, we are unable to determine unambiguously around which star each planetary candidate orbits (see Sects. 5.1 and 9). Thus, for the purposes of validating the planetary nature of the signals in this binary system, we adapted the original TRICERATOPS methodology (see equations 4 and 5 from Giacalone et al. 2021) by adding to the FPP computation the probabilities corresponding to transiting planets around the bound companion scenarios (STP and NTP), and correspondingly removing the NTP scenario from the NFPP calculation. Hence, our adjusted FPP and NFPP values accurately reflect the probabilities relevant to planetary validation in this particular binary system configuration.

Taking all these considerations into account, we obtained for TOI-2267.01 values of FPP=(5.6 ± 3.4) × 10−5 and NFPP=(4.8 ± 3.1) × 10−5. For TOI-2267.02, we derived FPP=(1.6 ± 0.6) × 10−4 and NFPP=(1.5 ± 0.6) × 10−4. Finally, for TOI-2267.03, the calculated values were FPP=(1.2 ± 0.5) × 10−5 and NFPP=(9.9 ± 4.6) × 10−6. These computed values for FPP and NFPP are notably low, clearly satisfying the established validation criteria for planetary candidates with high statistical confidence. Specifically, their corresponding false-positive probabilities translate to confidence levels exceeding 99.99% regarding their planetary nature.

thumbnail Fig. 6

SAI-2.5m contrast curve for the observation obtained on 2024 August 09 UT. The bump at 0.3″ occurs due to the companion. The autocorrelation function is given in the inset.

5 Global transit analysis

5.1 Derivation of the system parameters

We carried out our light-curve analyses using the allesfitter package (Günther & Daylan 2019, 2021), which allowed us to model planetary transits using the ellc package (Maxted 2016) while accounting for other phenomena such as stellar flares, spots, and variability. allesfitter also allows several ways to model the correlated noise, including polynomials, splines, and Gaussian processes (GPs; Rasmussen 2004), which are implemented through the celerite package (Foreman-Mackey et al. 2017a,b). The parameters of interest are retrieved using a Bayesian approach implementing a Markov chain Monte Carlo method (see, e.g., Hastings 1970; Ford 2005) using the emcee package (Foreman-Mackey et al. 2013), or the Nested Sampling inference algorithm (see, e.g., Feroz et al. 2009, 2019) using the dynesty package (Speagle 2020). We used the Dynamic Nested Sampling algorithm to estimate the Bayesian evidence in this study directly. This strategy allows us to compare a diverse set of orbital configurations robustly.

For each planet, we fitted the ratio of planetary radius over stellar radius (Rp/R⋆), the sum of the stellar and planetary radius scaled to the orbital semi-major axis ((Rp + R⋆)/a), the cosine of the orbital inclination (cos ip), the mid-transit time (T0), the orbital period (P), and, when considering eccentric scenarios, jointly the eccentricity and the argument of pericenter (e cos⁡ ω$\[\sqrt{e} ~\cos~ \omega\]$ and e sin⁡ ω$\[\sqrt{e} ~\sin~ \omega\]$). Given that the host star is part of a stellar binary system, the photometric data sets are contaminated by the companion star; in that regard, we included a dilution factor (D0). The dilution factor prior distributions were calculated from the flux ratio at each bandpass. The dilution factor for the Ic filter was obtained directly from high angular resolution observation (see Table 4), and we leave it fixed as a reference. The flux ratios for the rest of the bandpasses were estimated from the SED model. We report the flux-ratio values and the corresponding dilution factor priors in Table 5.

We used a hybrid spline to model the baseline for our ground-based observations. For more complex systematics, present in the data obtained by TESS, we used GPs to model the correlated noise using the stochastically driven damped harmonic oscillator (SHO) kernel, described by its three hyper-parameters; the frequency of the undamped oscillator (ω0), the quality factor of the oscillator (Q) and S0, which is related to the amplitude of the variability (further detail see, Foreman-Mackey et al. 2017a). allesfitter also fits an error scaling term to account for white noise in the data.

To reduce the number of free parameters in the models, we fixed the quadratic limb-darkening (LD) coefficients during the fitting process (see, e.g., Kipping et al. 2017; Günther et al. 2019). We computed the values of the LD coefficients in the physical μ-space, u1 and u2, for each bandpass used in the data interpolating from the tables of Claret et al. (2012). For the non-standard I + z filter, we took the averages of the values for the standard filters Ic and Sloan-z’. Then, following the relations of Kipping (2013)8, we converted u1 and u2 to the transformed q-space, q1 and q2 used in allesfitter. We report these values in Tables E.1 and E.2.

Before performing a global model accounting for all the available data, we analysed each light curve independently to estimate the GPs hyper-parameters and the white noise parameter (see for detailed description, Günther et al. 2019; Pozuelos et al. 2023). For the TESS data, we used the orbital periods and transit times obtained from the SPOC pipeline to refine the transit locations by performing a preliminary fit across all TESS sectors, applying wide, uniform priors. Next, we masked eight-hour windows around each identified transit midpoint and, using out-of-transit data, fit for noise and GPs hyper-parameters with wide uniform priors. Finally, we refined the planetary and orbital parameters by propagating the out-of-transit GPs posteriors as priors into a fit of the full TESS dataset, applying Gaussian distributions. This allowed us to sample the planetary and orbital parameters from wide uniform priors. We tested the configuration with eccentricity set to zero and as a free parameter in the TESS dataset fit. In both cases, we obtained practically the same results and no significant preference from the Bayesian evidence. Then, to reduce the number of free parameters, we fixed e cos⁡ ω$\[\sqrt{e} ~\cos~ \omega\]$ and e sin⁡ ω$\[\sqrt{e} ~\sin~ \omega\]$ to zero during the fitting process. In the context of planets in multiple-star systems, recent population analyses indicate that orbital eccentricities primarily correlate with the ratio between the star-star projected separation and the planetary semi-major axis, s/a, rather than with the absolute separation (González-Payo et al. 2024). For TOI-2267 planets, the s/a ratios range from 300 to 600, placing the system in the high s/a regime where no eccentricity enhancement is expected, consistent with our near-circular solutions. The presence of Earth-sized planets in this close binary therefore adds to the small but growing set of compact low-eccentricity architectures known in multiple systems.

The follow-up observations could present high levels of red noise. To determine which observations could effectively refine the transit parameters and which were dominated by red noise at scales larger than the transit signals, we independently fit each dataset using a pure-noise model and a transit-and-noise model, recording the Bayesian evidence for each. For the pure-noise model, we fit the light curves by setting the planetary radius ratio (Rp/R⋆) to zero (indicating no transit) and use a hybrid spline to model the red noise. In the transit-and-noise model, we fit the planetary and orbital parameters by sampling the posteriors from the TESS-only model using uniform distributions. We then computed the Bayes factor for each follow-up observation as Δ ln Z = ln Ztransit – ln Znoise. Strong evidence for a transit signal was identified in the data when Δ ln Z > 5 (Trotta 2008). Ground-based observations that met this criterion were directly included in the global analysis. The complete list of follow-up observations, along with their corresponding Δ ln Z values, is presented in Table D.1.

During the modelling procedure, allesfitter uses the stellar mass and radius obtained from our stellar characterisation (see Sect. 2.4) to compute a normal prior on the stellar density (see Table 2). During the fitting process, the stellar density is calculated at each Nested Sampling step from the fitted parameters via ρ⋆≈3πGP2(aR⋆)3$\[\rho_{\star} \approx \frac{3 \pi}{G P^{2}}\left(\frac{a}{R_{\star}}\right)^{3}\]$ (Seager & Mallén-Ornelas 2003) and compared with the value provided in the prior. The fits are penalised when these two values disagree.

Given the particular architecture of this system, before performing the global model combining TESS and ground-based observation, we masked each transit in TESS data to fit each candidate independently. Each of these single planets TESS-only fits were tested with the primary and the secondary star as a host. For each candidate fit, we computed the Bayes factor as Δ ln Zhost = ln Zprimary – ln Zsecondary. The TOI-2267.01 and TOI-2267.03 candidates showed a Bayes factor slightly favourable with the primary star as a host, while TOI-2267.02 showed a Bayes factor favourable with the secondary star. Despite these results, the Bayesian evidence was not strong enough to confirm the host star for each candidate (−1 < Δ ln Zhost < 1). As we describe in Sect. 7, the three planets orbiting the same star would lead to a highly unstable architecture. Given the close 3:2 mean motion resonance for TOI-2267.01 and TOI-2267.03, and the lack of ground-based confirmation for TOI-2267.02, we decided to joint-fit TOI-2267.01 and TOI-2267.03, including the ground-based observation. Then, the orbital fit of TOI-2267.02 was performed isolated from the other candidates using just TESS data.

The result of the global fit shows a preference for the configuration where planets TOI-2267.01 and .03 orbit the primary star, with a Δ ln Zhost ~ 0.5 favourable to primary star host scenario and host density from orbit in agreement at 1σ level with the host density from stellar characterisation (see Figs. I.1 and I.2). These results reveal a tentative configuration of the planetary system where TOI-2267.01 and .03 are genuine planets orbiting the primary star of the binary, TOI-2267A. Hence, hereafter, we denoted these planets as TOI-2267 b and TOI-2267 c. However, since this result is based on the limited Bayes evidence and 1σ agreement host density from orbit but not directly observed, future follow-up studies will focus on confirming this hypothesis. In Figs. 7 and 8, we display our best-fitting model using the TESS and the ground-based data. In Fig. 9, we show our best-fitting model using TESS for TOI-2267.02, which remains a planetary candidate due to the lack of independent confirmation using ground-based observations. The fitted and derived physical parameters for primary and secondary host cases for planets b, c, and the candidate TOI-2267.02 are reported in Tables H.1 and H.2.

Table 5

Flux ratios and dilution factors for TOI-2267.

thumbnail Fig. 7

Phase-folded and detrended photometry of TOI-2267 b transits along with the best-fit transit model (solid black line). The unbinned data points are shown in grey, while the blue circles with error bars correspond to 9-min bins.
thumbnail Fig. 8

Phase-folded and detrended photometry of TOI-2267 c transits along with the best-fit transit model (solid black line). The unbinned data points are shown in grey, while the red circles with error bars correspond to 9-min bins.
thumbnail Fig. 9

Phase-folded and detrended photometry of TOI-2267.02 transit along with the best-fit transit model (solid black line, secondary star host). The unbinned data points are shown in grey, while the dark blue circles with error bars correspond to 9-min bins.

5.2 Check for transit chromaticity

The early analyses of ground-based observations showed no evidence for wavelength-dependent transit depths. To verify this behaviour, we checked the diluted and undiluted transit depth values from our global analysis of the TOI-2267 b, and c planets. We compared the posterior distributions of the derived transit depths corresponding to the ground-based observations and TESS (see Tables 6 and J.1), finding that all undiluted depths, after dilution factor correction, agree at the 1σ level (see Figs. 10 and 11 for primary host, and Figs. J.1 and J.2 for secondary host). Hence, we confirm that the transits for both planets do not show any chromatic dependence, thus reaffirming the results and conclusions presented in the previous section.

6 Planet searches and detection limits from the TESS photometry

As discussed in Section 3.1, we executed our pipeline SHERLOCK when only the candidate TOI-2267.01 had been identified by the SPOC pipeline. SHERLOCK is a dedicated Python package designed to explore space-based photometric data in search of planetary signals, originally introduced in Pozuelos et al. (2020) and Demory et al. (2020). We refer the reader to Dévora-Pajares et al. (2024) for the latest version of the pipeline and a comprehensive overview of its functionalities.

During our first data exploration using Sectors 19, 20, 25, and 26, SHERLOCK identified a strong sinusoidal signal attributable to a fast rotator using a similar procedure as described in Section 2.3, with a maximum peak at 0.695 d (see Fig. 3). Such variability might hinder any planetary detection; hence, to remove it previously to the planetary search, SHERLOCK automatically used the cosine function provided by the wötan package (Hippke & Heller 2019). Once this correction was applied, SHERLOCK properly identified a signal with an orbital period of 3.5 d. This detection corresponds to the released TOI-2267.01, allowing us to recover this candidate independently. In the subsequent runs, we identified two unknown extra signals: one with an orbital period of 2.28 d and the other with an orbital period of 2.03 d. We executed our vetting module and found no evidence of any potential false positive origin. Each time a new sector became available, we executed SHERLOCK to check if more signals appeared; however, until our last execution, including the twelve sectors currently available with a 2 min cadence, only the three original signals were detected.

The apparent absence of additional candidates using TESS data, whether analysed by SHERLOCK or SPOC, could be attributed to one of the following scenarios (see, e.g., Wells et al. 2021; Schanche et al. 2022; Pozuelos et al. 2023): (1) no other planets exist in the system; (2) other planets exist but do not transit; (3) other planets exist and transit, but their orbital periods are longer than those explored in this study; or (4) other planets exist, and transit, but the photometric precision of the data is insufficient to detect them. Scenarios (1) and (2) could be further investigated through high-precision radial velocity follow-up, which is out of the scope of this study. Scenario (3) can be tested by extending the observational baseline. To evaluate the final scenario, we performed injection-and-recovery experiments using the MATRIX code9 (Dévora-Pajares & Pozuelos 2022).

This code conducts a three-dimensional parameter space exploration, generating a grid of orbital periods, planetary radii, and transit epochs. Each of these sets of parameters generates a synthetic scenario that is injected into the original light curve; in our case, we set up a grid of 60 periods (from 0.5 to 15 d), 30 radii (from 0.5 to 3 R⊕), and 10 different epochs, making a total of 18 000 scenarios. We note that the upper limit of 15 days was chosen not only to ensure computational feasibility but also because, under the assumption of a co-planar system, planets with longer periods would likely not transit (see, e.g., Jenkins et al. 2019). For each of these scenarios, the rotational signal at 0.70 d was corrected as we did with SHERLOCK, and then we applied a detrend with a bi-weight filter using a window size of 0.25 d, which was found to be the optimal length during the SHERLOCK exploration. Then, the light curves are processed in the search for planets using the same algorithm as SHERLOCK. A synthetic planet is considered as retrieved when its found period and epoch differ by at most 1% and up to 1 hour from the injected values, respectively. The execution results are shown in Fig. 12. On the one hand, the smallest detectable planets have sizes of 0.6 R⊕ and would be placed in the innermost orbits below ~1 d, as expected. That is, any planet smaller than that value would be undetectable for any orbital period. On the other hand, for any transiting planets with sizes larger than 1.5 R⊕ with orbital periods up to 15 d, we got 100% of recovery rate; hence, their existence can be discarded. Above an orbital period of 8 d, transiting Earth-sized planets begin to be missed at some epochs, becoming undetectable with very low recovery rates at periods ~11 d and above.

Table 6

Measured transit depths at each bandpass for the global fit assuming the primary star as the host.

thumbnail Fig. 10

Measured transit depth vs wavelength for TOI-2267A b. The dark grey horizontal line indicates the depth of the Johnson-Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); other filters are highlighted by circles for comparison. Top: diluted transit depths. Bottom: undiluted transit depths corrected by the dilution factor.
thumbnail Fig. 11

Measured transit depth vs wavelength for TOI-2267A c. The dark grey horizontal line indicates the depth of the Johnson–Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); other filters are highlighted by circles for comparison. Top: diluted transit depths. Bottom: undiluted transit depths corrected by the dilution factor.
thumbnail Fig. 12

MATRIX injection-and-recovery experiment conducted to establish the detection limits using the ten TESS sectors described in Sect. 3.1. We explored a total of 18 000 different scenarios, where each pixel evaluated 40 of them. Colours indicate the recovery rate: bright yellow is for a high recovery, and dark purple and black are for a low recovery. The solid blue line indicates the 95% recovery contour, the dashed white line marks the 50% recovery contour, and the solid white line corresponds to the 5% recovery contour. Blue, red, and green points show planets b and c and candidate .02, respectively. Black-edged circles denote planets orbiting the primary star; edge-free circles indicate the secondary star.

7 System architecture

7.1 Instability of the three-planet solution

Throughout the previous sections, we fully validated the planetary nature of the candidates TOI-2267 b, and c, but the candidate TOI-2267.02, while statistically validated, still remains elusive from a ground-based confirmation. According to our results, if this candidate becomes a validated planet, its orbital period would be 2.03 d, dangerously close to TOI-2267 b, whose orbital period is 2.28 d. Then, we decided to test the system’s stability, considering that three planets orbit the same star using the nominal values presented in Tables H.1 and H.2, and circular orbits for the three bodies. The planets’ masses are derived from the mass-radius relationship obtained using the Springht package (Parviainen et al. 2024). Since we do not know the actual physical distance between the two stellar components of TOI-2267, only the projected distance between them, we conducted two analyses with single-star scenarios. In the first case, the three planets orbit the primary star, and in the second, they orbit the secondary. For stability tests, we used the most recent version of the Stability Orbital Configuration Klassifier10 (SPOCK 2.0.0; Tamayo et al. 2020), a machine-learning model capable of classifying the stability of compact 3+ planetary systems over 109 orbits of the innermost planet, which in this case translates into ~107 yr (see, e.g., Pozuelos et al. 2023; Dreizler et al. 2024). To further enhance the robustness of our study, we performed a Monte Carlo analysis using 1000 realisations of the planetary masses drawn from their Spright-derived distributions. We then computed the corresponding stability probabilities with SPOCK. We found that for the scenario in which the three planets orbit the primary star, none of the simulations resulted in a stability probability greater than 0.9, and only 3% exceeded a threshold of 0.5. In the scenario where the planets orbit the secondary star (M2), this percentage increased slightly to 6.5%, but none reached a stability greater than 0.9. Hence, for both scenarios, we found highly unstable architectures, hinting that the three-planet configurations orbiting the same star are highly disfavoured. We foresee two potential solutions: the TOI-2267.02 is not a genuine planet, or, if real, the three planets do not orbit the same star.

7.2 Stability of the two-planet solution

Due to the high instability of the three-planet configuration orbiting a single star, this subsection examines the dynamical stability of two-planet combinations: TOI-2267.02 and b, TOI-2267.02 and c, and TOI-2267 b and c. Each scenario considers the two planets orbiting either the primary or secondary stellar component. As seen in Sect. 7.1, due to the lack of well-constrained separation between the two stars, each case is modelled as a single-star system hosting two planets.

To this end, we employ the Mean Exponential Growth factor of Nearby Orbits, Y(t) (MEGNO; Cincotta & Simó 1999, 2000; Cincotta et al. 2003), a well-established diagnostic for assessing the dynamical stability of planetary systems (see, e.g., Hinse et al. 2015; Jenkins et al. 2019; Delrez et al. 2021). The MEGNO parameter was computed using the REBOUND N-body integrator (Rein & Liu 2012), with the Wisdom-Holman WHfast symplectic algorithm (Rein & Tamayo 2015). The time-averaged value, ⟨Y(t)⟩, enhances stochastic variations in the orbital evolution, allowing a clear distinction between quasi-periodic motion (when ⟨Y(t → ∞)⟩ ≈ 2) and chaotic trajectories (when ⟨Y(t → ∞)⟩ → ∞).

For each two-planet configuration, 1000 realisations were generated by sampling orbital parameters from the posterior distributions in Tables H.1 and H.2. Planetary masses were sampled from the distributions provided by the SPRIGHT package. Two eccentricity regimes were considered: strictly circular orbits (e = 0) and orbits with eccentricities up to e ≤ 0.025. Each configuration was integrated for 107 orbital periods of the outermost planet, using a timestep equal to 5% of the innermost planet’s orbital period. A realisation was classified as dynamically unstable if its MEGNO value exceeded 2.1.

Our MEGNO analysis shows that configurations including planets .02 and b orbiting the same star are highly unstable across the explored parameter space. For strictly circular orbits, ~78% of the realisations are dynamically unstable, while allowing eccentricities up to e ≤ 0.025 increases the fraction to ~95%. In contrast, configurations pairing planet .02 with c, or b with c, are found to be highly stable. In both cases, the fraction of unstable realisations is 0% in the circular and eccentric regimes, indicating that these architectures remain stable for the entire integration timespan. These results strongly suggest that the actual system architecture consists of one of these two pairs orbiting the same stellar component, while the remaining planet (b or .02, respectively) orbits the other star.

The orbital periods of planets b and c lie close to a 3:2 mean-motion commensurability, whereas the .02–c pair is far from any low-order resonance. Proximity to such resonances can be the outcome of the system’s evolutionary history, as first-order resonances are frequently observed in compact multi-planet systems and are often interpreted as the signature of convergent migration within a protoplanetary disk (see, e.g., Bryden et al. 2000; Ramos et al. 2017; Huang & Ormel 2023; Wong & Lee 2024). In this scenario, planets exchange angular momentum with the surrounding disk material, and the smooth inward migration facilitates resonant capture when the orbits converge (e.g. Lee & Peale 2002; Batygin 2015). Resonant locking can then act as a stabilising mechanism, protecting the planets from close encounters despite their tight orbital spacing.

Given the stability of the b–c configuration, we extended our simulations to explore higher eccentricity regimes before investigating the possible resonant behaviour. Following Demory et al. (2020), we computed a MEGNO stability map of 150 × 150 pixels in the eb–ec parameter space, ranging from 0 to 0.3. The resulting stability map shows that the mutual eccentricities must satisfy eb ≤ 0.65–1.10 ec for long-term stability, implying that low eccentricities are strongly favoured. In particular, maximum mutual values of about eb ≃ ec ≃ 0.025 are allowed.

Using these constraints, we then examined the dynamical behaviour of the b–c system to test whether it is indeed in 3:2 resonance. We computed the apsidal evolution over 1 Myr for eccentricities between 0.001 and 0.025. Apsidal motion in interacting planetary systems can manifest as libration or circulation, separated by a secular separatrix (see, e.g., Barnes & Greenberg 2006b,c,a; Kane & Raymond 2014; Kane 2019). In all simulations, the b–c system exhibited apsidal libration around a fixed point offset from the origin (Fig. 13), with trajectories lying entirely in the positive ebec cos Δω direction, indicative of an aligned configuration. This behaviour confirms that planets b and c are likely locked in a stable, aligned libration around the 3:2 mean motion resonance (MMR), consistent with an evolutionary pathway involving resonant capture during disk-driven migration.

Although the specific stellar component hosting each planet cannot be determined unambiguously, the dynamical analysis supports a configuration in which planets b and c orbit the same star in a first-order mean-motion resonance, with planet .02 orbiting the companion star. This arrangement prevents rapid dynamical instability and aligns with a formation pathway involving convergent migration within a protoplanetary disk. However, the presence of two low-mass stars may increase the complexity of the system’s formation and subsequent dynamical evolution.

thumbnail Fig. 13

Polar plot of the apsidal trajectory (ebec vs Δω) for the TOI-2267 b and c planets over a period of 1 Myr. In this particular scenario, the planetary eccentricities were chosen as 0.025. The figure shows that the apsidal modes are librating around the 3:2 MMR (see Sect. 7.2).

7.3 Transit timing variations

As we have described previously, TOI-2267 b and c are close to the first order 3:2 MMR, an orbital architecture that might produce significant TTVs (see, e.g., Lithwick et al. 2012) allowing us to measure planetary masses, which is especially critical to deepening our understanding for a given planetary system when other techniques such as radial velocities are challenging or impracticable (Agol et al. 2021). To assess whether a dedicated photometric campaign could yield precise masses, we generated 10 000 synthetic system realisations following the approach of Pozuelos et al. (2023) for TOI-2096. Orbital periods and mid-transit epochs were drawn from Gaussian distributions centred on the values in Table H.1, planetary masses were sampled from the mass–radius priors produced by the Spright package, and eccentricities uniformly from 0 to 0.025 to ensure dynamical stability. For simplicity, and because the third candidate, TOI-2267.02, remains not fully validated and, if real, may orbit the companion star instead, we restrict our analysis to two two-planet hypotheses: one in which planets b and c orbit the primary component, and another in which they orbit the secondary.

We observed that the TTVs for each planet do not follow a symmetric distribution. Hence, we employed a Gaussian Kernel Density Estimation (KDE) to characterise these distributions using the gaussian-kde function from the scipy package (Virtanen et al. 2020). This non-parametric approach allowed us to accurately estimate the probability density functions (PDFs) of the TTVs without assuming an underlying distribution. When assuming that planets b and c orbit the primary star, the KDE modes occur at ~0.7 min and ~0.6 min, respectively, with PDF means of ~1.4 min (b) and ~1.0 min (c), and maximum TTV excursions of ~4 min (see Fig. 14). By contrast, under the alternative hypothesis that both planets orbit the secondary component, we obtain slightly larger amplitudes: KDE modes of 0.83 min (b) and 0.75 min (c), and means of 2.53 min (b) and 2.26 min (c), and maximum TTV excursions of ~10 min.

Based on these results, we conclude that measuring the planetary masses of TOI-2267 b and c using TTVs might be very challenging due to it requiring extremely high precision photometry with mid-transit times of about 0.5 min, which is out of reach of small- and mid-sized ground-based telescopes. In a recent work by Gillon et al. (2024), the authors reported the discovery of SPECULOOS-3b, an Earth-sized planet orbiting an M6 star that is a similar star-planet system to TOI-2267. In such a study, the HiPERCAM mounted on the GTC was used, which yielded a mid-transit precision of ~5.2 s, well below the required precision that we would need for TOI-2267 b and c. Hence, we conclude that GTC/HiPERCAM would be an ideal instrument for accurately measuring the TTVs and deriving the planetary masses, setting the stage for further studies using space-based facilities such as the James Webb Space Telescope.

thumbnail Fig. 14

Expected TTV amplitudes for planets TOI-2267 b (upper panel) and TOI-2267 c (lower panel), assuming both planets orbit the primary star. Dashed and solid vertical lines correspond to the mode and mean of the PDFs, respectively.

8 Discussion and prospects

8.1 TOI-2267 in the context of the binary systems with planets

Stellar multiplicity is a common outcome of star formation, with multiplicity fractions ranging from 40–50% for FGK stars and 20–30% for M-dwarfs (e.g. Duchêne & Kraus 2013; Winters et al. 2019; Clark et al. 2022a). However, among stars hosting planets, the occurrence of close stellar companions (projected separations <100 au) appears strongly suppressed. For instance, the POKEMON survey (Clark et al. 2022b, 2024b,a) found that M-dwarf binaries hosting planets tend to have projected separations peaking around 200 au, while those without known planets peak at ~6 au. Similar trends have been observed for solar-type stars (Hirsch et al. 2021). These findings suggest that close stellar companions may hinder planet formation or survival. Currently, approximately 8% of confirmed exoplanets reside in binary systems, with 50 such planets orbiting M-dwarfs, spanning 37 distinct systems11. Figure 15 shows these systems sorted by projected separation, using complementary estimates from Pecaut & Mamajek (2013) when some parameters of the secondary component were unavailable. Two systems, Kepler-779 and LHS 1678, were excluded from our final list due to insufficient information regarding the secondary component.

In this context, TOI-2267 is remarkable for its exceptionally tight configuration: a projected separation of only 8 au between its components, the smallest in our compiled sample. The only possible rival, LHS 1678, has an uncertain companion whose nature remains debated (planet vs. brown dwarf; Silverstein et al. 2022, 2024). Thus, TOI-2267 stands as the clearest example to date of an extremely compact low-mass binary system hosting planets, challenging the current view that such architectures inhibit planet formation or long-term stability.

thumbnail Fig. 15

List of binary systems where the stellar component hosting planets is a low-mass star with Teff < 4000 K. The systems are arranged according to the projected separation between the two stellar components of the binary, with TOI-2267 standing out as the closest binary system found to harbour planets. The sizes of the stars are to scale relative to each other, and their colours represent their stellar temperatures. Blue circles indicate transiting planets, which are also scaled relative to each other. In contrast, blue triangles denote non-transiting planets (detected by radial velocity only), all set to an arbitrary size. For each system, the orbital distances from the host star where the insolation flux ranges from 10 to 0.1 S⊕ are highlighted.

8.2 Planet formation scenario

At the dynamical level, circumstellar discs in binary systems exhibit significant differences from those around single stars (Cuello et al. 2025). Most notably, they tend to dissipate more rapidly, thereby shortening the timescale for planet formation to occur (Cieza et al. 2009; Kraus et al. 2012). In addition, solid material in these discs drifts inward more quickly, potentially leading to substantial mass loss unless mitigated by pressure traps or instabilities (Rosotti & Clarke 2018; Zagaria et al. 2023). Binarity also leads to disc truncation, yielding smaller and less massive discs (Pichardo et al. 2005). Consequently, if planet formation proceeds similarly in all environments, planets forming in binaries are expected to be more compact and less massive than those around single stars (Akeson et al. 2019; Quarles et al. 2020). For instance, under plausible assumptions about disc properties, sufficient material may have existed to form several rocky planets in the Alpha Centauri system (Cuello & Sucerquia 2024). These findings suggest that planet formation in binaries may be less inhibited than previously thought, and that rocky planets in such systems could be common.

Overall, the study of S-type exoplanets (those orbiting one star in a binary) is still maturing. Comprehensive reviews by Marzari & Thebault (2019) and Bonavita & Desidera (2020) highlight the progress and limitations in this field. Current observational techniques remain biased, particularly in their ability to detect rocky planets in binaries, which hinders a comprehensive understanding of this population. Nonetheless, the potential for discovering multi-planet systems in binaries is growing thanks to an increase in sensitivity and longer monitoring campaigns. Our results on the TOI-2267 binary system, with at least one transiting planet candidate around each component, mark a significant step forward. This demonstrates the feasibility of detecting such systems and offers new opportunities to test planet formation models. Enhancing detection capabilities and correcting for survey incompleteness will be crucial for developing a more comprehensive view of planet occurrence in binary systems.

8.3 Disentangling the planetary hosts

As described in Section 5.1, a comparison of Bayesian evidence from transit fits accounting for the binary’s blended light dilution factor and the implied host star densities suggests that TOI-2267 b and c orbit the primary star. However, this indirect statistical evidence cannot unambiguously assign the planets to their host. One potential method to assign the TOI-2267 planets to a host star is by directly detecting the transits through high-precision time-series photometry with simultaneous extremely high spatial resolution (<0.4″) that resolves the binary components. This method has been demonstrated previously for the Kepler-13 AB system (Howell et al. 2019), which hosts much brighter stars and has deeper transits. We scaled the photometric precision obtained for Kepler-13 with the ‘Alopeke instrument on Gemini-N, and found that it would be impossible to resolve the binary and directly detect the shallow TOI-2267 transits using this technique for even the largest current ground-based telescopes. Instead, this direct transit detection could be made from space, as HST/WFC3 or JWST/NIRCam have sufficient technical capabilities. Both instruments have a sufficiently wide FOV and angular resolution to resolve the full PSFs of TOI-2267A and B, and a high enough photometric precision to detect the transits. Any space-based transit observations of the TOI-2267 planets would also provide exquisite transit timing measurements that could be combined with ground-based follow-up to better constrain the planetary masses through TTVs.

If TOI-2267.02 is confirmed to be a genuine planet, observations to assign the stellar host could reveal the first-ever double-transiting binary system architecture. There are a few known binary star systems with planets orbiting both components (González-Payo et al. 2024), but none in which both components have transiting planets assigned. Hence, the TOI-2267 system may present the only feasible opportunity to characterise a double transiting binary system, with significant implications for studying planet formation scenarios in close binary systems.

9 Conclusion

Through this work, we have identified, validated, and provided a preliminary characterisation of the nearby system TOI-2267, which is composed of two very low-mass stars with at least two warm Earth-sized exoplanets. We first conducted a detailed characterisation of this peculiar cool binary. To do so, we leveraged SED fitting, spectroscopic analyses, stellar atmosphere modelling, and evolutionary models. These analyses revealed that the binary is composed of M5V and M6V stars, with effective temperatures of about 3030 and 2930 K, respectively. The system has a projected separation of around 8 au. This makes it the coolest binary with the smallest stellar projected separation known to host planets.

We then validated the planetary nature of two of the three planetary candidates independently found by our team and TESS/SPOC. To achieve this, we combined photometry from ground-based facilities, data from 12 TESS sectors, archival information, high angular resolution imaging, and statistical considerations. As a result, we validated two planets, TOI-2267 b and TOI-2267 c, though we cannot yet unambiguously determine which star they orbit. For the third Earth-sized candidate, TOI-2267.02, our comprehensive validation analysis strongly supports its planetary nature, but it was not detected using ground-based observations. Despite favourable statistical evidence, we prudently keep it as a planetary candidate pending further confirmation.

When studying the system’s architecture, we found an interesting result. If TOI-2267.02 is a genuine planet, placing all three planets around the same star becomes rapidly unstable. In addition, planets b and .02 cannot orbit the same star. The most stable solutions are the pairs b-c or .02-c orbiting the same star, while the other planet orbits the other component. Since planets TOI-2267 b and c are close to the 3:2 first-order MMR, the scenario with b-c orbiting one star and .02 orbiting the other is the most plausible.

This hypothesis makes the TOI-2267 system a benchmark case for further studies, as it would represent the first detection of transiting planets orbiting both components of an extremely compact binary. This unique configuration provides an opportunity to investigate how planetary formation and orbital dynamics operate under the influence of two close stellar hosts.

References

  1. Agol, E., Dorn, C., Grimm, S. L., et al. 2021, Planet. Sci. J., 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
  2. Akeson, R. L., Jensen, E. L. N., Carpenter, J., et al. 2019, ApJ, 872, 158 [Google Scholar]
  3. Allard, F., Homeier, D., & Freytag, B. 2012, Philos. Trans. Roy. Soc. London Ser. A, 370, 2765 [Google Scholar]
  4. Aller, A., Lillo-Box, J., Jones, D., Miranda, L. F., & Barceló Forteza, S. 2020, A&A, 635, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437 [Google Scholar]
  6. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  8. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  9. Barkaoui, K., Burdanov, A., Hellier, C., et al. 2019, AJ, 157, 43 [NASA ADS] [CrossRef] [Google Scholar]
  10. Barnes, R., & Greenberg, R. 2006a, ApJ, 652, L53 [Google Scholar]
  11. Barnes, R., & Greenberg, R. 2006b, ApJ, 638, 478 [Google Scholar]
  12. Barnes, R., & Greenberg, R. 2006c, ApJ, 647, L163 [Google Scholar]
  13. Batygin, K. 2015, MNRAS, 451, 2589 [NASA ADS] [CrossRef] [Google Scholar]
  14. Blake, C. H., Charbonneau, D., & White, R. J. 2010, ApJ, 723, 684 [Google Scholar]
  15. Bonavita, M., & Desidera, S. 2020, Galaxies, 8, 16 [NASA ADS] [CrossRef] [Google Scholar]
  16. Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [Google Scholar]
  17. Brasseur, C. E., Phillip, C., Fleming, S. W., Mullally, S. E., & White, R. L. 2019, Astrocut: Tools for creating cutouts of TESS images, Astrophysics Source Code Library [record ascl:1905.007] [Google Scholar]
  18. Brown, T. M., Baliber, N., Bianco, F. B., et al. 2013, PASP, 125, 1031 [Google Scholar]
  19. Bryden, G., Różyczka, M., Lin, D. N. C., & Bodenheimer, P. 2000, ApJ, 540, 1091 [Google Scholar]
  20. Burdanov, A. Y., de Wit, J., Gillon, M., et al. 2022, PASP, 134, 105001 [CrossRef] [Google Scholar]
  21. Burgasser, A. J., Logsdon, S. E., Gagné, J., et al. 2015, ApJS, 220, 18 [Google Scholar]
  22. Cepa, J., Aguiar, M., Escalera, V. G., et al. 2000, SPIE Conf. Ser., 4008, 623 [Google Scholar]
  23. Charbonneau, D., & Deming, D. 2007, arXiv e-prints [arXiv:0706.1047] [Google Scholar]
  24. Childs, A. C., Martin, R. G., & Livio, M. 2022, ApJ, 937, L41 [Google Scholar]
  25. Ciardi, D. R., Beichman, C. A., Horch, E. P., & Howell, S. B. 2015, ApJ, 805, 16 [NASA ADS] [CrossRef] [Google Scholar]
  26. Cieza, L. A., Padgett, D. L., Allen, L. E., et al. 2009, ApJ, 696, L84 [NASA ADS] [CrossRef] [Google Scholar]
  27. Cincotta, P., & Simó, C. 1999, Celest. Mech. Dyn. Astron., 73, 195 [NASA ADS] [CrossRef] [Google Scholar]
  28. Cincotta, P. M., & Simó, C. 2000, A&ASS, 147, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D Nonlinear Phenomena, 182, 151 [NASA ADS] [CrossRef] [Google Scholar]
  30. Claret, A. 2017, A&A, 600, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Claret, A., Hauschildt, P. H., & Witte, S. 2012, A&A, 546, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Clark, C. A., van Belle, G. T., Ciardi, D. R., et al. 2022a, AJ, 163, 232 [Google Scholar]
  33. Clark, C. A., van Belle, G. T., Horch, E. P., et al. 2022b, AJ, 164, 33 [Google Scholar]
  34. Clark, C. A., van Belle, G. T., Horch, E. P., et al. 2024a, AJ, 167, 174 [Google Scholar]
  35. Clark, C. A., van Belle, G. T., Horch, E. P., et al. 2024b, AJ, 167, 56 [Google Scholar]
  36. Collins, K. A., Kielkopf, J. F., Stassun, K. G., & Hessman, F. V. 2017, AJ, 153, 77 [Google Scholar]
  37. Cuello, N., & Sucerquia, M. 2024, Universe, 10, 64 [NASA ADS] [CrossRef] [Google Scholar]
  38. Cuello, N., Alaguero, A., & Poblete, P. P. 2025, arXiv e-prints [arXiv:2501.19249] [Google Scholar]
  39. Cutri, R. M. et al. 2014, VizieR Online Data Catalog: II/328 [Google Scholar]
  40. Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
  41. Delrez, L., Gillon, M., Queloz, D., et al. 2018, Proc. SPIE, 10700, 107001I [Google Scholar]
  42. Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., 5, 775 [NASA ADS] [CrossRef] [Google Scholar]
  43. Demory, B. O., Pozuelos, F. J., Gómez Maqueo Chew, Y., et al. 2020, A&A, 642, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Dévora-Pajares, M. & Pozuelos, F. J. 2022, https://doi.org/10.5281/zenodo.6570831 [Google Scholar]
  45. Dévora-Pajares, M., Pozuelos, F. J., Thuillier, A., et al. 2024, MNRAS, 532, 4752 [CrossRef] [Google Scholar]
  46. Dhillon, V. S., Bezawada, N., Black, M., et al. 2021, MNRAS, 507, 350 [NASA ADS] [CrossRef] [Google Scholar]
  47. Dietrich, J., Apai, D., Schlecker, M., et al. 2023, AJ, 165, 149 [NASA ADS] [CrossRef] [Google Scholar]
  48. Dittmann, J. A., Irwin, J. M., Charbonneau, D., & Newton, E. R. 2016, ApJ, 818, 153 [NASA ADS] [CrossRef] [Google Scholar]
  49. Dreizler, S., Jeffers, S. V., Rodríguez, E., et al. 2020, MNRAS, 493, 536 [Google Scholar]
  50. Dreizler, S., Luque, R., Ribas, I., et al. 2024, A&A, 684, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45 [Google Scholar]
  52. Duchêne, G., & Kraus, A. 2013, ARAA, 51, 269 [Google Scholar]
  53. Fernandes, C. S., Van Grootel, V., Salmon, S. J. A. J., et al. 2019, ApJ, 879, 94 [NASA ADS] [CrossRef] [Google Scholar]
  54. Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601 [NASA ADS] [CrossRef] [Google Scholar]
  55. Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, Open J. Astrophys., 2, 10 [Google Scholar]
  56. Ford, E. B. 2005, AJ, 129, 1706 [Google Scholar]
  57. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  58. Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017a, celerite: Scalable 1D Gaussian Processes in C++, Python, and Julia, Astrophysics Source Code Library [record ascl:1709.008] [Google Scholar]
  59. Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017b, AJ, 154, 220 [Google Scholar]
  60. Furlan, E., & Howell, S. B. 2017, AJ, 154, 66 [NASA ADS] [CrossRef] [Google Scholar]
  61. Furlan, E., & Howell, S. B. 2020, AJ, 898, 47 [Google Scholar]
  62. Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  63. Garcia-Mejia, J., Charbonneau, D., Fabricant, D., et al. 2020, SPIE Conf. Ser., 11445, 114457R [NASA ADS] [Google Scholar]
  64. Garcia, L. J., Timmermans, M., Pozuelos, F. J., et al. 2022, MNRAS, 509, 4817 [Google Scholar]
  65. Giacalone, S., Dressing, C. D., Jensen, E. L. N., et al. 2021, AJ, 161, 24 [Google Scholar]
  66. Gibbs, A., Bixel, A., Rackham, B. V., et al. 2020, AJ, 159, 169 [NASA ADS] [CrossRef] [Google Scholar]
  67. Gillon, M. 2018, Nat. Astron., 2, 344 [Google Scholar]
  68. Gillon, M., Jehin, E., Magain, P., et al. 2011, in European Physical Journal Web of Conferences, 11, 6002 [Google Scholar]
  69. Gillon, M., Pedersen, P. P., Rackham, B. V., et al. 2024, Nat. Astron., 8, 865 [Google Scholar]
  70. González-Payo, J., Caballero, J. A., Gorgas, J., et al. 2024, A&A, 689, A302 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39 [NASA ADS] [CrossRef] [Google Scholar]
  72. Günther, M. N., & Daylan, T. 2019, Allesfitter: Flexible Star and Exoplanet Inference From Photometry and Radial Velocity, Astrophysics Source Code Library [record ascl:1903.003] [Google Scholar]
  73. Günther, M. N., & Daylan, T. 2021, ApJSS, 254, 13 [Google Scholar]
  74. Günther, M., Pozuelos, F., Dittmann, J., et al. 2019, Nat. Astron., 3, 1099 [Google Scholar]
  75. Hastings, W. K. 1970, Biometrika, 57, 97 [Google Scholar]
  76. Henry, T. J., Jao, W.-C., Subasavage, J. P., et al. 2006, AJ, 132, 2360 [Google Scholar]
  77. Hinse, T. C., Haghighipour, N., Kostov, V. B., & Goździewski, K. 2015, ApJ, 799, 88 [NASA ADS] [CrossRef] [Google Scholar]
  78. Hippke, M., & Heller, R. 2019, A&A, 623, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  79. Hirsch, L. A., Rosenthal, L., Fulton, B. J., et al. 2021, AJ, 161, 134 [NASA ADS] [CrossRef] [Google Scholar]
  80. Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [Google Scholar]
  81. Howell, S. B., Everett, M. E., Sherry, W., Horch, E., & Ciardi, D. R. 2011, AJ, 142, 19 [Google Scholar]
  82. Howell, S. B., Scott, N. J., Matson, R. A., Horch, E. P., & Stephens, A. 2019, AJ, 158, 113 [Google Scholar]
  83. Hsu, C.-C., Burgasser, A. J., Theissen, C. A., et al. 2021, ApJS, 257, 45 [NASA ADS] [CrossRef] [Google Scholar]
  84. Hsu, C.-C., Theissen, C., Burgasser, A., & Birky, J. 2021, SMART: The Spectral Modeling Analysis and RV Tool, Tech. rep. [Google Scholar]
  85. Hsu, C.-C., Burgasser, A. J., & Theissen, C. A. 2023, ApJ, 945, L6 [Google Scholar]
  86. Hsu, C.-C., Burgasser, A. J., Theissen, C. A., et al. 2024, ApJS, 274, 40 [NASA ADS] [CrossRef] [Google Scholar]
  87. Huang, S., & Ormel, C. W. 2023, MNRAS, 522, 828 [NASA ADS] [CrossRef] [Google Scholar]
  88. Jang-Condell, H. 2015, ApJ, 799, 147 [NASA ADS] [CrossRef] [Google Scholar]
  89. Jehin, E., Gillon, M., Queloz, D., et al. 2011, The Messenger, 145, 2 [NASA ADS] [Google Scholar]
  90. Jenkins, J. M. 2002, ApJ, 575, 493 [Google Scholar]
  91. Jenkins, J. M., Chandrasekaran, H., McCauliff, S. D., et al. 2010, Proc. SPIE, 7740, 77400D [NASA ADS] [CrossRef] [Google Scholar]
  92. Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proc. SPIE, 9913, 99133E [Google Scholar]
  93. Jenkins, J. S., Pozuelos, F. J., Tuomi, M., et al. 2019, MNRAS, 490, 5585 [Google Scholar]
  94. Jenkins, J. S., Díaz, M. R., Kurtovic, N. T., et al. 2020, Nat. Astron., 4, 1148 [Google Scholar]
  95. Kaib, N. A., Raymond, S. N., & Duncan, M. 2013, Nature, 493, 381 [Google Scholar]
  96. Kane, S. R. 2019, AJ, 158, 72 [NASA ADS] [CrossRef] [Google Scholar]
  97. Kane, S. R., & Raymond, S. N. 2014, ApJ, 784, 104 [NASA ADS] [CrossRef] [Google Scholar]
  98. Kipping, D. M. 2013, MNRAS, 435, 2152 [Google Scholar]
  99. Kipping, D. M., Cameron, C., Hartman, J. D., et al. 2017, AJ, 153, 93 [NASA ADS] [CrossRef] [Google Scholar]
  100. Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32 [Google Scholar]
  101. Kraus, A. L., Ireland, M. J., Hillenbrand, L. A., & Martinache, F. 2012, ApJ, 745, 19 [Google Scholar]
  102. Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [Google Scholar]
  103. Lester, K. V., Matson, R. A., Howell, S. B., et al. 2021, AJ, 162, 75 [NASA ADS] [CrossRef] [Google Scholar]
  104. Li, Y. P., Gan, W. Q., & Feng, L. 2012, ApJ, 747, 133 [NASA ADS] [CrossRef] [Google Scholar]
  105. Lightkurve Collaboration (Cardoso, J. V. d. M., et al.) 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library [record ascl:1812.013] [Google Scholar]
  106. Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122 [Google Scholar]
  107. Manara, C. F., Tazzari, M., Long, F., et al. 2019, A&A, 628, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in Astronomical Society of the Pacific Conference Series, 534, Protostars and Planets VII, eds. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 539 [Google Scholar]
  109. Martin, E. C., Fitzgerald, M. P., McLean, I. S., et al. 2018, SPIE Conf. Ser., 10702, 107020A [NASA ADS] [Google Scholar]
  110. Marzari, F., & Thebault, P. 2019, Galaxies, 7, 84 [NASA ADS] [CrossRef] [Google Scholar]
  111. Maxted, P. F. L. 2016, A&A, 591, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  112. McCully, C., Volgenau, N. H., Harbeck, D.-R., et al. 2018, SPIE Conf. Ser., 10707, 107070K [Google Scholar]
  113. McLean, I. S., Becklin, E. E., Bendiksen, O., et al. 1998, SPIE Conf. Ser., 3354, 566 [NASA ADS] [Google Scholar]
  114. McLean, I. S., Graham, J. R., Becklin, E. E., et al. 2000, SPIE Conf. Ser., 4008, 1048 [Google Scholar]
  115. Minkowski, R. L., & Abell, G. O. 1963, in Basic Astronomical Data: Stars and Stellar Systems, ed. K. A. Strand, 481 [Google Scholar]
  116. Moehler, S., Modigliani, A., Freudling, W., et al. 2014, A&A, 568, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  117. Monet, D. G., Levine, S. E., Canzian, B., et al. 2003, AJ, 125, 984 [Google Scholar]
  118. Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93 [Google Scholar]
  119. Nogueira, P. H., Zurlo, A., Pérez, S., et al. 2023, MNRAS, 523, 4970 [NASA ADS] [CrossRef] [Google Scholar]
  120. Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, in Astronomical Society of the Pacific Conference Series, 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 275 [Google Scholar]
  121. O’Malley-James, J. T., & Kaltenegger, L. 2019, MNRAS, 485, 5598 [CrossRef] [Google Scholar]
  122. Parviainen, H., Luque, R., & Palle, E. 2024, MNRAS, 527, 5693 [Google Scholar]
  123. Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9 [Google Scholar]
  124. Pichardo, B., Sparke, L. S., & Aguilar, L. A. 2005, MNRAS, 359, 521 [NASA ADS] [CrossRef] [Google Scholar]
  125. Pozuelos, F. J., Suárez, J. C., de Elía, G. C., et al. 2020, A&A, 641, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  126. Pozuelos, F. J., Timmermans, M., Rackham, B. V., et al. 2023, A&A, 672, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  127. Quarles, B., Li, G., Kostov, V., & Haghighipour, N. 2020, AJ, 159, 80 [Google Scholar]
  128. Ragusa, E., Fasano, D., Toci, C., et al. 2021, MNRAS, 507, 1157 [NASA ADS] [CrossRef] [Google Scholar]
  129. Ramos, X. S., Charalambous, C., Benítez-Llambay, P., & Beaugé, C. 2017, A&A, 602, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. Rasmussen, C. E. 2004, Gaussian Processes in Machine Learning, eds. O. Bousquet, U. von Luxburg, & G. Rätsch (Berlin, Heidelberg: Springer Berlin Heidelberg), 63 [Google Scholar]
  131. Redfield, S., Batalha, N., Benneke, B., et al. 2024, arXiv e-prints [arXiv:2404.02932] [Google Scholar]
  132. Rein, H., & Liu, S.-F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  133. Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [Google Scholar]
  134. Reylé, C., Jardine, K., Fouqué, P., et al. 2021, A&A, 650, A201 [Google Scholar]
  135. Ribas, I., Reiners, A., Zechmeister, M., et al. 2023, A&A, 670, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  136. Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, SPIE Conf. Series, 9143, 20 [Google Scholar]
  137. Rodrigo, C., & Solano, E. 2020, in XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society, 182 [Google Scholar]
  138. Rosotti, G. P., & Clarke, C. J. 2018, MNRAS, 473, 5630 [Google Scholar]
  139. Savel, A. B., Dressing, C. D., Hirsch, L. A., et al. 2020, AJ, 160, 287 [NASA ADS] [CrossRef] [Google Scholar]
  140. Schanche, N., Pozuelos, F. J., Günther, M. N., et al. 2022, A&A, 657, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Scott, N. J., Howell, S. B., Gnilka, C. L., et al. 2021, Front. Astron. Space Sci., 8, 138 [NASA ADS] [CrossRef] [Google Scholar]
  142. Seager, S., & Mallén-Ornelas, G. 2003, ApJ, 585, 1038 [Google Scholar]
  143. Sebastian, D., Gillon, M., Ducrot, E., et al. 2021, A&A, 645, A100 [EDP Sciences] [Google Scholar]
  144. Serna, J., Hernandez, J., Kounkel, M., et al. 2021, ApJ, 923, 177 [NASA ADS] [CrossRef] [Google Scholar]
  145. Shields, A. L., Ballard, S., & Johnson, J. A. 2016, PhysRep, 663, 1 [Google Scholar]
  146. Silverstein, M. L., Schlieder, J. E., Barclay, T., et al. 2022, AJ, 163, 151 [NASA ADS] [CrossRef] [Google Scholar]
  147. Silverstein, M. L., Barclay, T., Schlieder, J. E., et al. 2024, AJ, 167, 255 [NASA ADS] [CrossRef] [Google Scholar]
  148. Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000 [Google Scholar]
  149. Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
  150. Stassun, K. G., & Torres, G. 2016, AJ, 152, 180 [Google Scholar]
  151. Stassun, K. G., & Torres, G. 2021, ApJ, 907, L33 [NASA ADS] [CrossRef] [Google Scholar]
  152. Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136 [Google Scholar]
  153. Stassun, K. G., Corsaro, E., Pepper, J. A., & Gaudi, B. S. 2018, AJ, 155, 22 [Google Scholar]
  154. Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102 [Google Scholar]
  155. Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182 [EDP Sciences] [Google Scholar]
  156. Strakhov, I. A., Safonov, B. S., & Cheryasov, D. V. 2023, Astrophys. Bull., 78, 234 [NASA ADS] [CrossRef] [Google Scholar]
  157. Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985 [Google Scholar]
  158. Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100 [Google Scholar]
  159. Tamayo, D., Cranmer, M., Hadden, S., et al. 2020, PNAS, 117, 18194 [CrossRef] [PubMed] [Google Scholar]
  160. Theissen, C. A., Konopacky, Q. M., Lu, J. R., et al. 2022, ApJ, 926, 141 [Google Scholar]
  161. Tran, H. D., Cohen, R., Colson, A., et al. 2016, Proc. SPIE, 9910, 99102E [Google Scholar]
  162. TRAPPIST-1 JWST Community Initiative (de Wit, J., et al.) 2024, Nat. Astron., 8, 810 [Google Scholar]
  163. Triaud, A. H. M. J. 2021, in ExoFrontiers; Big Questions in Exoplanetary Science, ed. N. Madhusudhan, 6 [Google Scholar]
  164. Trotta, R. 2008, Contemp. Phys., 49, 71 [Google Scholar]
  165. Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502 [Google Scholar]
  166. Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30 [Google Scholar]
  167. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  168. Weber, P., Pérez, S., Guidi, G., et al. 2023, MNRAS, 518, 5620 [Google Scholar]
  169. Wells, R. D., Rackham, B. V., Schanche, N., et al. 2021, A&A, 653, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  170. Winters, J. G., Henry, T. J., Jao, W.-C., et al. 2019, AJ, 157, 216 [NASA ADS] [CrossRef] [Google Scholar]
  171. Wong, K. H., & Lee, M. H. 2024, AJ, 167, 112 [NASA ADS] [CrossRef] [Google Scholar]
  172. Zagaria, F., Rosotti, G. P., Alexander, R. D., & Clarke, C. J. 2023, Eur. Phys. J. Plus, 138, 25 [NASA ADS] [CrossRef] [Google Scholar]
  173. Zahnle, K. J., & Catling, D. C. 2017, ApJ, 843, 122 [NASA ADS] [CrossRef] [Google Scholar]
  174. Zúñiga-Fernández, S., Burdanov, A., Delrez, L., Gillon, M., & Jehin, E. 2025, SPECULOOS Exoplanet Search and Its Prototype on TRAPPIST, eds. H. J. Deeg, & J. A. Belmonte (Cham: Springer Nature Switzerland), 1 [Google Scholar]
  175. Zúñiga-Fernández, S., Gillon, M., Pedersen, P. P., et al. 2024, in Ground-based and Airborne Telescopes X, 13094, eds. H. K. Marshall, J. Spyromilio, & T. Usuda, International Society for Optics and Photonics (SPIE), 1309413 [Google Scholar]
  176. Zurlo, A., Cieza, L. A., Pérez, S., et al. 2020, MNRAS, 496, 5089 [Google Scholar]

Appendix A Acknowledgements

We thank the anonymous referee for providing constructive comments and suggestions, which helped improve the clarity and quality of this manuscript. We sincerely thank Steven Giacalone for his helpful advice and guidance on properly using TRICERATOPS. F.J.P and P.J.A acknowledge financial support from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033 and Ministerio de Ciencia e Innovación through the project PID2022-137241NB-C43. Partially based on observations made at the Observatorio de Sierra Nevada (OSN), operated by the Instituto de Astrofísica de Andalucía (IAA-CSIC). M.G. is F.R.S.-FNRS Research Director. J.K. acknowledge support from the Swiss NCCR PlanetS and the Swiss National Science Foundation. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40182901 and 51NF40205606. J.K. gratefully acknowledges the support of the Swedish National Space Agency (SNSA; DNR 2020-00104) and of the Swedish Research Council (Project Grant 2017-04945 and 2022-04043) and of the Swiss National Science Foundation under grant number TMSGI2_211697. M.N.G. acknowledges support from the European Space Agency (ESA) as an ESA Research Fellow. The postdoctoral fellowship of KB is funded by F.R.S.-FNRS grant T.0109.20 and by the Francqui Foundation. N.C. and M.S. acknowledge funding support from the European Research Council (ERC) under the European Union Horizon Europe program (grant agreement No. 101042275, project Stellar-MADE). B.R-A acknowledges funding support from the ANID Basal project FB210003. The material is based upon work supported by NASA under award number 80GSFC24M0006. This material is based upon work supported by the National Aeronautics and Space Administration under Agreement No. 80NSSC21K0593 for the program “Alien Earths”. The results reported herein benefited from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. YGMC acknowledges support from UNAM-PAPIIT-101224. We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through project PID2021-125627OB-C32, and from the Centre of Excellence “Severo Ochoa” award to the Instituto de Astrofisica de Canarias. F. M. acknowledges financial support from the Agencia Estatal de Investigación del Ministerio de Ciencia, Innovación y Universidades (MCIU/AEI) through grant PID2023-152906NA-I00. J.d.W. and MIT gratefully acknowledge financial support from the Heising-Simons Foundation, Dr. and Mrs. Colin Masson and Dr. Peter A. Gilman for Artemis, the first telescope of the SPECULOOS network situated in Tenerife, Spain. H.P. acknowledges support from the Spanish Ministry of Science and Innovation with the Ramon y Cajal fellowship number RYC2021-031798-I, and funding from the University of La Laguna and the Spanish Ministry of Universities. Funding for the TESS mission is provided by NASA’s Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This paper includes data collected by the TESS mission that are publicly available from the Mikulski Archive for Space Telescopes (MAST). Based on data collected by the SPECULOOS-North Observatory at the Observatorio del Teide in Canary Islands, Spain. The ULiege’s contribution to SPECULOOS has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) (grant Agreement n° 336480/SPECULOOS), from the Balzan Prize and Francqui Foundations, from the Belgian Scientific Research Foundation (F.R.S.-FNRS; grant n° T.0109.20), from the University of Liege, and from the ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. This work is supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). This research is in part funded by the European Union’s Horizon 2020 research and innovation program (grants agreements n° 803193/BEBOP), and from the Science and Technology Facilities Council (STFC; grant n° ST/S00193X/1, and ST/W000385/1). This work is based on observations taken with the SAINT-EX telescope at the Observatorio Astronómico Nacional on the Sierra de San Pedro Mártir (OAN-SPM), Baja California, México, that is supported by the Swiss National Science Foundation (SNSF; PP00P2-163967 and PP00P2-190080, SPIRIT-216537), the Centre for Space and Habitability (CSH) of the University of Bern, the National Centre for Competence in Research PlanetS, supported by the SNSF, and UNAM PAPIITIG101224. The research leading to these results has received funding from the ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. TRAPPIST is funded by the Belgian Fund for Scientific Research (Fond National de la Recherche Scientifique, FNRS) under the grant PDR T.0120.21. TRAPPIST-North is a project funded by the University of Liege (Belgium), in collaboration with Cadi Ayyad University of Marrakech (Morocco). Based on observations obtained at the Observatoire du Mont-Mégantic, financed by Université de Montréal, Université Laval, the Canada Economic Development program and the Ministère de l’Economie et de l’Innovation du Québec. This work benefited from support of the Fonds de recherche du Québec – Nature et technologies (FRQNT), through the Center for Research in Astrophysics of Quebec. Some of the observations in this paper made use of the High-Resolution Imaging instrument ‘Alopeke and were obtained under Gemini LLP Proposal Number: GN/S-2021A-LP-105. ‘Alopeke was funded by the NASA Exoplanet Exploration Program and built at the NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. ‘Alopeke was mounted on the Gemini North telescope of the international Gemini Observatory, a program of NSF’s OIR Lab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). Based on observations made with the GTC telescope, in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, under Director’s Discretionary Time (GTC10-23BDDT, IP F.J.P). IAS acknowledges the support of M.V. Lomonosov Moscow State University Program of Development. This research made use of Lightkurve, a Python package for Kepler, and TESS data analysis (Lightkurve Collaboration et al. 2018). This work made use of Astropy:12, a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022). This work made use of allesfitter, a public and user-friendly astronomy software package for modelling photometric and radial velocity data (Günther & Daylan 2021, 2019). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region. This research has made use of data obtained from or tools provided by the portal exoplanet.eu of The Extrasolar Planets Encyclopaedia.

Appendix B SED modelling

Table B.1

Fitted parameters correlation matrix.

Appendix C NIRSPEC spectra modelling

We briefly summarise our modelling method, following the analysis detailed in Hsu et al. (2021). We employed the forward-modelling method (Blake et al. 2010; Burgasser et al. 2015) that fits our spectra with a single stellar template and telluric absorption simultaneously with the SMART (Spectral Modelling Analysis and RV Tool) package (Hsu et al. 2021). We chose the BT-Settl models (Allard et al. 2012) and the ESO earth atmosphere models (Moehler et al. 2014) to fit our NIRSPEC spectra. We focused our analysis on order 33 (2.29–2.34 μm) because this order has the CO (2–0) rotational-vibrational bandhead, ideal for measuring precise RVs and v sin is values. We used the Markov chain Monte Carlo sampling method to derive our best-fit parameters using the emcee package (Foreman-Mackey et al. 2013), for a total of 9 parameters, 50 chains, and 600 steps, with the first 200 steps as burn-ins.

Appendix D Ground-based photometric observations

Table D.1

Ground-based time-series photometric observations logs of TOI-2267.

Appendix E Quadratic limb darkening coefficients

Table E.1

Quadratic limb-darkening coefficients for TOI-2267A in each bandpass reported in Table D.1.

Table E.2

Quadratic limb-darkening coefficients for TOI-2267B in each bandpass reported in Table D.1.

Appendix F TESS PCDSAP light curves

thumbnail Fig. F.1

TESS photometric time series of TOI-2267 obtained for sectors 19, 20, 25, 26, 40, 52, 53, 59, 60, 73, 79, and 86. In all cases, the grey points correspond to the PDCSAP fluxes obtained from the SPOC pipeline. The blue, red and green points correspond to the location of the transits for the candidates TOI-2267.01, TOI-2267.02 and TOI-2267.03, respectively.

Appendix G TESS field of view

thumbnail Fig. G.1

TESS target pixel files (TPFs) showing the field of view of the 12 sectors used in this study, generated employing tpfplotter (Aller et al. 2020). The photometric apertures used to extract the light curves in each case are shown as red-shaded regions. The Gaia DR2 catalogue is overplotted, with all sources up to six magnitudes in contrast with TOI-2267 shown as red circles. We note that the symbol size scales with the magnitude contrast. The target star, TOI-2267, is highlighted with a white cross and numbered as 1 in all cases.

Appendix H Global fit

Table H.1

Parameters for the TOI-2267 b and TOI-2267 c planets.

thumbnail Fig. H.1

Posterior probability distributions for all the physical parameters fitted for TOI-2267 bc systems (Primary star host) using allesfitter nested sampling as described in Sect. 5.1. The vertical dashed lines represent the median and the 68% credible interval. The figure highlights the correlation (or absence thereof) between all the parameters.
thumbnail Fig. H.2

Posterior probability distributions for all the physical parameters fitted for TOI-2267 bc systems (Secondary star host) using allesfitter nested sampling as described in Sect. 5.1. The vertical dashed lines represent the median and the 68% credible interval. The figure highlights the correlation (or absence thereof) between all the parameters.
Table H.2

Parameters for the TOI-2267.02 candidate.

thumbnail Fig. H.3

Posterior probability distributions for all the physical parameters fitted for TOI-2267.02 using allesfitter nested sampling as described in Sect. 5.1. Upper panel for the primary star host and lower panel for the secondary star host case. The vertical dashed lines represent the median and the 68% credible interval. The figure highlights the correlation (or absence thereof) between all the parameters.

Appendix I Host density from orbital fit

thumbnail Fig. I.1

Comparison between host density obtained from stellar model and orbital fitting for TOI-2267 bc for primary (Top panel) and secondary (Bottom panel) host.
thumbnail Fig. I.2

Comparison between host density obtained from stellar model and orbital fitting for TOI-2267.02 for primary (Top panel) and secondary (Bottom panel) host

Appendix J Transit depths

thumbnail Fig. J.1

Measured transit depth vs wavelength for TOI-2267 b (secondary host). The dark grey horizontal line indicates the depth of the Johnson-Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); circles highlight other filters for comparison. Top: Diluted transit depths. Bottom: Undiluted transit depths corrected by the dilution factor.
Table J.1

Measured transit depths at each bandpass for the global fit assuming the secondary star as the host.

thumbnail Fig. J.2

Measured transit depth vs wavelength for TOI-2267 c (secondary host). The dark grey horizontal line indicates the depth of the Johnson-Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); other filters are highlighted by circles for comparison. Top: Diluted transit depths. Bottom: Undiluted transit depths corrected by the dilution factor.
2

According to NASA Exoplanet Archive as of January 2025.

3

According to the Encyclopaedia of exoplanetary systems as of June 2025; https://exoplanet.eu/planets_binary/

8

q1 = (u1 + u2)2 and q2 = 0.5u1/(u1 + u2).

9

The MATRIX (Multi-phAse Transits Recovery from Injected eXoplanets) code is open access on GitHub: https://github.com/PlanetHunters/tkmatrix

11

According to the NASA Exoplanet Archive consulted in July 2025.

All Tables

Table 1

TOI-2267 stellar designations and astrometric properties.

In the text
Table 2

Derived properties of the primary and secondary stars of the TOI-2267 system.

In the text
Table 3

Photometric properties of TOI-2267.

In the text
Table 4

Speckle interferometric observations of TOI-2267 at SAI-2.5 m.

In the text
Table 5

Flux ratios and dilution factors for TOI-2267.

In the text
Table 6

Measured transit depths at each bandpass for the global fit assuming the primary star as the host.

In the text
Table B.1

Fitted parameters correlation matrix.

In the text
Table D.1

Ground-based time-series photometric observations logs of TOI-2267.

In the text
Table E.1

Quadratic limb-darkening coefficients for TOI-2267A in each bandpass reported in Table D.1.

In the text
Table E.2

Quadratic limb-darkening coefficients for TOI-2267B in each bandpass reported in Table D.1.

In the text
Table H.1

Parameters for the TOI-2267 b and TOI-2267 c planets.

In the text
Table H.2

Parameters for the TOI-2267.02 candidate.

In the text
Table J.1

Measured transit depths at each bandpass for the global fit assuming the secondary star as the host.

In the text

All Figures

thumbnail Fig. 1

Spectral energy distribution of TOI-2267. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the best-fit NextGen stellar atmosphere model for the two stellar components (hot component in blue, cool component in red, combined light in black).
In the text
thumbnail Fig. 2

Keck/NIRSPEC K-band spectra of TOI-2267 around 2.3 μm. Top: observed NIRSPEC spectrum in the nodding position A (shown with the black line) with the data noise in the grey-shaded region. The stellar models with and without telluric absorption features are indicated with the blue and green lines, respectively. The residual (data – model) is denoted with the magenta line. Bottom: same as the top panel for NIRSPEC spectra and model fit but in the nodding position B.
In the text
thumbnail Fig. 3

Periodogram obtained with TESSExtractor tool for TOI-2267 from TESS Sector 53.
In the text
thumbnail Fig. 4

Archival images for TOI-2267. From left to right: 1955 DSS/POSS-I image, 1996 DSS/POSS-II image, 2011 PanSTARRS image, and 2025 Artemis-1m0 image. The red circle corresponds to the target position at 2022 November 09.
In the text
thumbnail Fig. 5

‘Alopeke speckle imaging 5σ contrast curves along with the reconstructed 832-nm image.
In the text
thumbnail Fig. 6

SAI-2.5m contrast curve for the observation obtained on 2024 August 09 UT. The bump at 0.3″ occurs due to the companion. The autocorrelation function is given in the inset.
In the text
thumbnail Fig. 7

Phase-folded and detrended photometry of TOI-2267 b transits along with the best-fit transit model (solid black line). The unbinned data points are shown in grey, while the blue circles with error bars correspond to 9-min bins.
In the text
thumbnail Fig. 8

Phase-folded and detrended photometry of TOI-2267 c transits along with the best-fit transit model (solid black line). The unbinned data points are shown in grey, while the red circles with error bars correspond to 9-min bins.
In the text
thumbnail Fig. 9

Phase-folded and detrended photometry of TOI-2267.02 transit along with the best-fit transit model (solid black line, secondary star host). The unbinned data points are shown in grey, while the dark blue circles with error bars correspond to 9-min bins.
In the text
thumbnail Fig. 10

Measured transit depth vs wavelength for TOI-2267A b. The dark grey horizontal line indicates the depth of the Johnson-Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); other filters are highlighted by circles for comparison. Top: diluted transit depths. Bottom: undiluted transit depths corrected by the dilution factor.
In the text
thumbnail Fig. 11

Measured transit depth vs wavelength for TOI-2267A c. The dark grey horizontal line indicates the depth of the Johnson–Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); other filters are highlighted by circles for comparison. Top: diluted transit depths. Bottom: undiluted transit depths corrected by the dilution factor.
In the text
thumbnail Fig. 12

MATRIX injection-and-recovery experiment conducted to establish the detection limits using the ten TESS sectors described in Sect. 3.1. We explored a total of 18 000 different scenarios, where each pixel evaluated 40 of them. Colours indicate the recovery rate: bright yellow is for a high recovery, and dark purple and black are for a low recovery. The solid blue line indicates the 95% recovery contour, the dashed white line marks the 50% recovery contour, and the solid white line corresponds to the 5% recovery contour. Blue, red, and green points show planets b and c and candidate .02, respectively. Black-edged circles denote planets orbiting the primary star; edge-free circles indicate the secondary star.
In the text
thumbnail Fig. 13

Polar plot of the apsidal trajectory (ebec vs Δω) for the TOI-2267 b and c planets over a period of 1 Myr. In this particular scenario, the planetary eccentricities were chosen as 0.025. The figure shows that the apsidal modes are librating around the 3:2 MMR (see Sect. 7.2).
In the text
thumbnail Fig. 14

Expected TTV amplitudes for planets TOI-2267 b (upper panel) and TOI-2267 c (lower panel), assuming both planets orbit the primary star. Dashed and solid vertical lines correspond to the mode and mean of the PDFs, respectively.
In the text
thumbnail Fig. 15

List of binary systems where the stellar component hosting planets is a low-mass star with Teff < 4000 K. The systems are arranged according to the projected separation between the two stellar components of the binary, with TOI-2267 standing out as the closest binary system found to harbour planets. The sizes of the stars are to scale relative to each other, and their colours represent their stellar temperatures. Blue circles indicate transiting planets, which are also scaled relative to each other. In contrast, blue triangles denote non-transiting planets (detected by radial velocity only), all set to an arbitrary size. For each system, the orbital distances from the host star where the insolation flux ranges from 10 to 0.1 S⊕ are highlighted.
In the text
thumbnail Fig. F.1

TESS photometric time series of TOI-2267 obtained for sectors 19, 20, 25, 26, 40, 52, 53, 59, 60, 73, 79, and 86. In all cases, the grey points correspond to the PDCSAP fluxes obtained from the SPOC pipeline. The blue, red and green points correspond to the location of the transits for the candidates TOI-2267.01, TOI-2267.02 and TOI-2267.03, respectively.
In the text
thumbnail Fig. G.1

TESS target pixel files (TPFs) showing the field of view of the 12 sectors used in this study, generated employing tpfplotter (Aller et al. 2020). The photometric apertures used to extract the light curves in each case are shown as red-shaded regions. The Gaia DR2 catalogue is overplotted, with all sources up to six magnitudes in contrast with TOI-2267 shown as red circles. We note that the symbol size scales with the magnitude contrast. The target star, TOI-2267, is highlighted with a white cross and numbered as 1 in all cases.
In the text
thumbnail Fig. H.1

Posterior probability distributions for all the physical parameters fitted for TOI-2267 bc systems (Primary star host) using allesfitter nested sampling as described in Sect. 5.1. The vertical dashed lines represent the median and the 68% credible interval. The figure highlights the correlation (or absence thereof) between all the parameters.
In the text
thumbnail Fig. H.2

Posterior probability distributions for all the physical parameters fitted for TOI-2267 bc systems (Secondary star host) using allesfitter nested sampling as described in Sect. 5.1. The vertical dashed lines represent the median and the 68% credible interval. The figure highlights the correlation (or absence thereof) between all the parameters.
In the text
thumbnail Fig. H.3

Posterior probability distributions for all the physical parameters fitted for TOI-2267.02 using allesfitter nested sampling as described in Sect. 5.1. Upper panel for the primary star host and lower panel for the secondary star host case. The vertical dashed lines represent the median and the 68% credible interval. The figure highlights the correlation (or absence thereof) between all the parameters.
In the text
thumbnail Fig. I.1

Comparison between host density obtained from stellar model and orbital fitting for TOI-2267 bc for primary (Top panel) and secondary (Bottom panel) host.
In the text
thumbnail Fig. I.2

Comparison between host density obtained from stellar model and orbital fitting for TOI-2267.02 for primary (Top panel) and secondary (Bottom panel) host
In the text
thumbnail Fig. J.1

Measured transit depth vs wavelength for TOI-2267 b (secondary host). The dark grey horizontal line indicates the depth of the Johnson-Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); circles highlight other filters for comparison. Top: Diluted transit depths. Bottom: Undiluted transit depths corrected by the dilution factor.
In the text
thumbnail Fig. J.2

Measured transit depth vs wavelength for TOI-2267 c (secondary host). The dark grey horizontal line indicates the depth of the Johnson-Cousin Ic filter (dilution factor fixed in the global fit, see Sect. 5.1); other filters are highlighted by circles for comparison. Top: Diluted transit depths. Bottom: Undiluted transit depths corrected by the dilution factor.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.