EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 667 (November 2022) A&A, 667 (2022) A59 Full HTML
Issue
A&A
Volume 667, November 2022
Article Number A59
Number of page(s) 31
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/202244041
Published online 08 November 2022

© L. Delrez et al. 2022

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

One main goal of modern astronomy is the identification and atmospheric characterisation of temperate terrestrial exoplanets, to understand how frequently and under which conditions life may exist around other stars. Terrestrial planets transiting nearby late-type M dwarfs are key in this endeavor. Indeed, for a given planet (fixed radius, mass, and equilibrium temperature), the signal-to-noise ratio (S/N) of atmospheric spectral features probed by eclipse (transit or occultation) spectroscopy increases for smaller and cooler host stars (e.g. Kaltenegger & Traub 2009; de Wit & Seager 2013). As a consequence, fewer eclipse observations need to be co-added to achieve a significant atmospheric detection. In this regard, the low luminosity of the latest-type M stars is also a considerable advantage, as it results in more frequent planetary eclipses for a given stellar irradiation, meaning that it takes much less time to obtain a certain amount of ineclipse observations. The S/N of atmospheric spectral features also scales inversely with the distance from the Earth to the host star, making the nearest late-type M dwarfs the best targets to search for transiting temperate terrestrial planets whose atmospheres could be characterised with current or next-generation facilities.

These considerations motivated the development of the SPECULOOS project (Search for habitable Planets EClipsing ULtra-cOOl Stars; Gillon 2018; Burdanov et al. 2018; Delrez et al. 2018), an exoplanet transit survey targeting a volume-limited (40 pc) sample of about 1700 late-type dwarfs with spectral type M6 and later – mostly stars, but also ~5% brown dwarfs (Sebastian et al. 2021) – using a network of 1m-class robotic telescopes. While SPECULOOS started its scientific operations officially in 2019 (Jehin et al. 2018), it was initiated in 2011 as a prototype survey (Gillon et al. 2013) targeting fifty of the brightest southern late-type M dwarfs with the TRAPPIST-South telescope. This prototype survey led to the discovery of the TRAPPIST-1 system1, consisting of seven temperate Earth-sized planets transiting a nearby (12 pc) M8V star (Gillon et al. 2016, 2017; Luger et al. 2017). The discovery of this benchmark system has caused a flurry of theoretical and observational follow-up studies, so that the TRAPPIST-1 planets are today the best-studied terrestrial planets outside our Solar System (Agol et al. 2021). They are also the most favourable targets found so far in the temperate terrestrial regime for atmospheric characterisation with the recently-launched JWST (see, e.g., Lustig-Yaeger et al. 2019; Gillon et al. 2020) and are already approved for substantial (~200 h) Cycle 1 Guaranteed Time Observations and General Observers programs.

Besides TRAPPIST-1, only two other transiting systems are currently known around M dwarfs with spectral type M5 or later: LHS 3844 (Vanderspek et al. 2019) and LP 791-18 (Crossfield et al. 2019). LHS 3844 b is an ultra-short-period super-Earth (1.3 R⊕) that orbits its M5V star every 11 hours. LP 791-18 is an M6V star transited by a super-Earth (1.1 R⊕) and a sub-Neptune (2.3 R⊕) with respective orbital periods of 0.95 and 4.99 days. Both systems were detected by the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015). While TESS is efficient at detecting small transiting planets around early-to-mid-type M dwarfs (e.g. Günther et al. 2019; Kostov et al. 2019; Demory et al. 2020; Wells et al. 2021; Schanche et al. 2022; Gan et al. 2022), its detection potential drops sharply for objects later than ~M5 (Sullivan et al. 2015; Barclay et al. 2018; Sebastian et al. 2021), due to their faintness in TESS’s bandpass. The two systems mentioned above demonstrate that TESS is able to detect some short-period and/or ≳2 R⊕ planets around some bright late-type M stars. However, possible additional transiting planets in the same systems with longer orbital periods and/or smaller radii may be missed due to the limited photometric precision. With its near-infrared-optimised cameras and custom ‘I + z’ filters (transmittance >90% from 750 to ~1100 nm), SPECULOOS is specifically designed to achieve high photometric precision on late-type/ultra-cool dwarf stars and can detect those additional transiting planets, thus allowing an effective synergy between the two surveys.

In this paper, we present a discovery that leverages that synergy: the detection of two temperate super-Earths transiting the nearby M6-type dwarf star LP 890-9. The inner planet was first detected by TESS based on four sectors of data. The announcement of this planet candidate (identified by TESS as TOI-4306.01) triggered intensive photometric monitoring of the target with the SPECULOOS Southern Observatory, which led to the discovery of a second longer-period transiting planet (identified by SPECULOOS as SPECULOOS-2 c), previously undetected by TESS. The orbital period of this second planet was later confirmed by MuSCAT3 follow-up observations. We validated the planetary nature of the system through follow-up observations from several ground-based facilities, including high-precision multi-colour photometry, spectroscopy, high-angular-resolution imaging, and archival images.

The paper is structured as follows. In Sect. 2, we describe the contributing facilities and datasets used in the validation and characterisation of the system. We infer the properties of the host star in Sect. 3 and validate the planetary nature of the two transit signals in Sect. 4. In Sect. 5, we present our detailed analysis of the photometry, including a global transit analysis to derive the system’s properties (Sect. 5.1), a search for additional transiting planets in the TESS and SPECULOOS-South data and an assessment of their detection limits (Sect. 5.2), as well as a study of the stellar variability (Sect. 5.3). Section 6 describes our analysis of Subaru/IRD radial velocities and derivation of preliminary mass constraints. In Sect. 7, we present a dynamical analysis of the system, including studies of its tidal evolution and long-term stability. Finally, we discuss the system’s properties, potential habitability, and prospects for follow-up in Sect. 8, before concluding in Sect. 9.

2 Observations

In this section, we present all the observations of LP 890-9 obtained with TESS and ground-based follow-up facilities.

2.1 TESS photometry

LP 890-9 (TIC 44898913) is part of the TESS Cool Dwarf catalogue (Muirhead et al. 2018). It was identified as an attractive target for a transit search with TESS, and thus included in the TESS Candidate Target List (Stassun et al. 2018b). The star was observed at a two-minute cadence during sectors 4-5 (18 October-11 December 2018) of the primary mission and again two years later during sectors 31–32 (21 October-17 December 2020) of the extended mission. The observations were acquired with CCDs 1 (sectors 5 and 32) and 2 (sectors 4 and 31) on Camera 2. The data were processed with the TESS Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016) at NASA Ames Research Center and searched for periodic transit signals (Jenkins 2002; Jenkins et al. 2010, Jenkins et al. 2020). While the analysis of the first two sectors did not identify any convincing signal, the addition of two more sectors during the extended mission revealed a ~0.7%-deep transit-like signature at a period of 2.73 days, with a Multiple Event Statistic (MES) of 8.3. The TESS Science Office reviewed the SPOC Data Validation Report (Twicken et al. 2018; Li et al. 2019) for this signal and announced the planet candidate TOI-4306.01 on 21 July 2021 (Guerrero et al. 2021). Figure 1 shows the target pixel files and photometric apertures used by SPOC for each of the four sectors, with the locations of nearby Gaia DR2 sources, up to 6 magnitudes in contrast with LP 890-9. None of them lie within the SPOC photometric apertures.

We retrieved the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP, Stumpe et al. 2012, 2014; Smith et al. 2012) from the Mikulski Archive for Space Telescopes (MAST) and removed all data points for which the quality flag was not zero. The resulting light curves are shown in Fig. 2, together with the transit times. As noted in the data release notes2 for sector 4, communications between the instrument and spacecraft were interrupted between BJD times 2 458 418.54 and 2 458 421.21, during which time no data or telemetry were collected, hence a gap in the light curve. The photometry obtained after this instrument anomaly was impacted by some thermal effects, as the camera temperature increased by ~20 degrees before returning to nominal within three days. We thus decided to exclude this region of the light curve (marked in orange in Fig. 2) from our analysis.

thumbnail Fig. 1

TESS target pixel files of LP 890-9 for sectors 4 (upper left), 5 (upper right), 31 (lower left), and 32 (lower right). The pixel scale is 21″ per pixel. The electron counts are colour-coded. The target is indicated by a white cross. Nearby sources identified in Gaia DR2 (Gaia Collaboration 2018), up to 6 magnitudes in contrast with LP 890-9, are shown as red circles. The symbol size is proportional to the magnitude contrast with the target. The apertures used by the SPOC pipeline to extract the photometry are shown as red shaded regions. These plots were created with tpfplotter (Aller et al. 2020).

2.2 Ground-based follow-up photometry

We obtained follow-up photometry of multiple transit events of the candidate TOI-4306.01 with several ground-based facilities and different bandpasses, as part of TESS Follow-up Observing Program (TFOP) Sub-Group 1 (SG1) for Seeing-Limited Photometry. The goals of these observations were to confirm that the transit signal detected by TESS was on the expected star, assess its chromaticity, and obtain higher-precision transit light curves compared to the TESS data. In addition to observations of transits associated with TOI-4306.01, we also monitored the target intensively with the SPECULOOS Southern Observatory (see Sect. 2.2.1), to search for possible additional transiting planets that may not have been detected by TESS and study the variability of the star. As mentioned previously, this monitoring revealed a second longer-period transiting planet candidate (identified by SPECULOOS as SPECULOOS-2 c), whose orbital period was later confirmed by MuSCAT3 follow-up observations (see Sect. 2.2.4). We describe all the ground-based photometric observations in the following sections and summarise the transit light curves of both planets in Table 1.

2.2.1 SPECULOOS-South

The SPECULOOS Southern Observatory (SSO, Burdanov et al. 2018; Delrez et al. 2018; Gillon 2018; Jehin et al. 2018) consists of four 1-m telescopes (Io, Europa, Ganymede, and Callisto) located at ESO Paranal Observatory in Chile. Each telescope is equipped with a deep-depletion 2K × 2K CCD detector optimised for the near infrared, providing a 12′ – 12′ field of view with a pixel scale of 0.35″ per pixel.

Fourteen transit light curves of TOI-4306.01 were obtained with SSO in several filters (see Table 1): three in r′, two in i′, two in z′, and seven in the special ‘I + z’ filter. As mentioned above, we also performed an intense photometric monitoring of the system outside of these transits, gathering in total 614.45 hours of observations spread over 119 nights between 9 August 2021 and 20 January 2022. For this monitoring, we used the I + z filter with an exposure time of 39 s to optimise the S/N. These data revealed a second transit signal, with a transit depth of about 0.6%. Three transits of this second object (identified by SPECULOOS as SPECULOOS-2c) were observed with SSO (see Table 1), all in the I + z filter. Each pair of consecutive transits was separated by 16.914 days. The only fractional period alias that remained possible based on the SSO data was 8.457 days (16.914/2 days). Shorter period aliases were excluded, as they would have produced several other detectable transits in the SSO data. The period of 8.457 days, which was ~1.7 times more likely than the 16.914-days period based on geometric transit probability, was later confirmed thanks to photometric follow-up observations from Hawaii using MuSCAT3 (see Sect. 2.2.4).

The data were processed with the automatic SSO pipeline, which is described in detail in Murray et al. (2020). This pipeline first performs standard image reduction steps, applying bias, dark, and flat-field corrections. Star detection, astrometric solving, and aperture photometry are then conducted using the casutools package (Irwin et al. 2004). For the observations presented here, the selected aperture radii are between 8 and 16 pixels (2.8–5.6″), depending on the filter and night. The target’s light curve is then generated by an automated differential photometry algorithm, which uses a weighted ensemble of comparison stars to correct for most atmospheric and instrumental effects. Finally, this light curve is corrected for the effects of time-varying telluric water vapour (see also Pedersen et al., in prep.), which can significantly affect near-infrared differential photometry of very red objects such as LP 890-9.

thumbnail Fig. 2

TESS photometry of LP 890-9. For each of the four sectors, the 2-min data points (in grey) have been binned into 30-min intervals to produce the black points, with error bars corresponding to the root-mean-square of the uncertainties of the points in the bins. The transits of LP 890-9 b and c are indicated by red and blue dotted lines, respectively. The region marked in orange in sector 4 was impacted by thermal effects and thus excluded from our analysis.
Table 1

Ground-based transit light curves of LP 890-9 used in our global transit analysis.

2.2.2 TRAPPIST-South

We observed a transit of TOI-4306.01 on 10 August 2021 with the 0.6-m TRAPPIST-South telescope (Gillon et al. 2011; Jehin et al. 2011), located at ESO La Silla Observatory in Chile. The wide ‘blue-blocking’ filter (transmittance >90% from 520 to ~1100 nm) was used for these observations, with an exposure time of 120 s. The data were processed using the prose open-source3 Python framework, which is described in detail in Garcia et al. (2022). We used a photometric aperture of 8.7 pixels (5.6’’).

2.2.3 ExTrA

ExTrA (Bonfils et al. 2015) is a near-infrared (0.85–1.55 µm) multi-object spectrograph fed by three 60-cm telescopes located at La Silla Observatory in Chile. One full transit of TOI-4306.01 was observed on 20 September 2021 using two ExTrA telescopes. We used 8″ aperture fibres and the low-resolution mode (R ~ 20) of the spectrograph with an exposure time of 60 s. Five fibre positioners are used at the focal plane of each telescope to select light from the target and four comparison stars. As comparison stars, we observed 2MASS J04172127-2815315, 2MASS J04191176-2812449, 2MASS J04181010-2753179, and 2MASS J04164312-2759006 with 2MASS J-magnitudes (Skrutskie et al. 2006) and Gaia effective temperatures (Gaia Collaboration 2018) similar to the target. The resulting ExTrA data were analysed using a custom data reduction software, described in more detail in Cointepas et al. (2021).

2.2.4 MuSCAT3

We observed three transits of TOI-4306.01 in the g′, r′, i′, and zs bands simultaneously by using the multi-band imager MuSCAT3 on the Faulkes Telescope North of Las Cumbres Observatory on Haleakala, Hawaii (Narita et al. 2020). MuSCAT3 is equipped with four 2k × 2k CCD cameras for the four channels, each providing a pixel scale of 0.266″ pixel−1 and a field of view of 9.1′ × 9.1′. The first transit was observed on 28 September 2021 (UT) with exposure times of 120, 120, 110, and 100 s in the g′, r′, i′ and zs bands, respectively. The second transit was observed on 28 October 2021 (UT) with exposure times of 120, 120, 55, and 29 s in the g′, r′, i′ and zs bands, respectively. The third transit was taken on 8 November 2021 (UT) with the same exposure times as the second transit. All observations were conducted without defocussing. On the night of the first, second, and third transits, the lunar phase was 21, 22, and 3 days, and the angular distance between the target and the moon was 53, 70, and 125 degrees, respectively.

We also observed one transit of the second planet candidate detected by SSO (SPECULOOS-2 c) with MuSCAT3 on 17 January 2022 (UT). The goal of these observations was to firmly confirm the orbital period of 8.457 days, since the less likely 16.914-days period would not have produced a transit during that night. We used exposure times of 120, 120, 55, and 29 s for the g′, r′, i′ and zs bands, respectively. The lunar phase and the moon distance were 14 days and 64 degrees, respectively. The sky background was bright because the almost-full moon was located nearby. Thus, the photometric precision in the g′ and r′ bands was not as good as for the previous observations of TOI-4306.01. Nonetheless, the transit was clearly detected in the i′ and z′ bands (see Fig. 10), thus confirming the orbital period of 8.457 days. We note that we did not use the g′-band light curve in our global transit analysis (see Sect. 5.1), as its scatter (~2.9%) was too large for it to be useful.

The obtained images were calibrated by the BANZAI pipeline (McCully et al. 2018). We performed aperture photometry on the calibrated images using a custom pipeline (Fukui et al. 2011). We optimised the aperture radius and combination of comparison stars for each night and each band so that the dispersion of the light curve was minimised. The adopted aperture radii were 1014 pixels, depending on the band and night.

2.2.5 SAINT-EX

A partial transit of TOI-4306.01 was observed on 30 November 2021 with the 1-m SAINT-EX telescope, located at the Observatorio Astronómico Nacional in the Sierra de San Pedro Mártir in Mexico (Sabin et al. 2018; Demory et al. 2020). SAINT-EX is a twin of the SSO telescopes but the camera is operated using a different readout mode and higher gain so that it can observe brighter stars. For these observations, we used the I + z filter with an exposure time of 128 s. The data were reduced using the custom PRINCE pipeline, which is detailed in Demory et al. (2020). A weighted principal component analysis (PCA) approach (Bailey 2012) was used to correct the light curve for systematics that are common between the target and comparison stars.

2.3 Spectroscopy

In this section, we describe the spectroscopic observations of LP 890-9 that we obtained to characterise the star (low-resolution near-infrared and optical spectra, Sects. 2.3.1 and 2.3.2) and to measure radial velocities (high-resolution near-infrared spectra, Sect. 2.3.3).

2.3.1 IRTF/SpeX

We obtained a low-resolution near-infrared spectrum of LP 890-9 with the SpeX spectrograph (Rayner et al. 2003) on the 3.2-m NASA Infrared Telescope Facility (IRTF) on two nights, 30 August 2021 and 23 November 2021 (UT). The conditions on both nights were clear with seeing of 0.7−0.8″. We used the short-wavelength cross-dispersed (SXD) mode with the 0.3″ × 15″ slit aligned to the parallactic angle, giving a 0.8−2.5 µm spectrum with a resolving power of R ~ 2000. In August, we collected three ABBA nod sequences (12 exposures) with integration times of 179.8 s per exposure, giving a total integration time of 36 min. In November, we collected a single ABBA nod sequence with integration times of 599.6 s per exposure, giving a total integration time of 40 min. We collected a set of standard SXD flat-field and arc-lamp exposures immediately after the science exposures, followed by a set of six exposures of the A0 V star HD 27016 (V = 7.77) in August and the A0 V star HD 22687 (V = 8.55) in November for flux and telluric calibration. We reduced the data using Spextool v4.1 (Cushing et al. 2004), following the standard instructions in the Spextool User’s Manual4. The final spectra have a median S/N of ~100, with peaks in the J, H, and K bands of 126, 136, and 126, respectively.

2.3.2 Lick/Kast

We obtained a low-resolution optical spectrum of LP 890-9 with the Kast Double Spectrograph on the 3m Shane Telescope at Lick Observatory on 14 November 2021 (UT) in clear conditions with 2 seeing. This instrument provides simultaneous blue and red optical spectra with parallel spectrograph cameras. Our analysis focussed on red optical data obtained with the 600/7500 grating and 1.5″-wide slit, providing 6000-9000 Å spectra at an average resolution of λ/Δλ ≈ 2500. The source was observed with 6 exposures of 600 s each within an hour of transit and at an average airmass of 2.5. The flux calibrator Feige 110 (Hamuy et al. 1992, 1994) was observed on the same night, and dome illuminated flat field and arc lamp exposures were obtained at the start of the night for pixel response and wavelength calibration. Data were reduced using the KastRedux package5 (Burgasser et al. in prep.), and included image reduction (flat-fielding, bad pixel masking, and linearity correction), optimal extraction, wavelength calibration, and flux density calibration. No attempt was made to correct for telluric absorption. The data have a median S/N of ~50 at 7400 Å.

2.3.3 Subaru/IRD

To measure radial velocities (RVs) for LP 890-9, we obtained high-resolution near-infrared spectra with the InfraRed Doppler (IRD) spectrograph mounted on the Subaru 8.2-m telescope (Tamura et al. 2012; Kotani et al. 2018). IRD is a fibre-fed, temperature-stabilised spectrograph, having a wavelength coverage of 970–1730 nm (Y, J, and H-bands) with a spectral resolution of ≈70000. For the simultaneous wavelength reference, laser-frequency comb (LFC) light was routinely injected into IRD using a secondary fibre. A total of 14 frames were secured for LP 890-9 between 10 September 2021 and 10 January 2022 (UT). Due to instrumental trouble in one of the two IRD detectors, only H-band spectra (> 1450 nm) were obtained for the three frames taken in January 2022. Given the faintness of the target (J = 12.258), the exposure times were set to 1800 s for all frames.

We reduced the raw IRD data as per the reduction procedure in Hirano et al. (2020), and extracted both stellar and LFC spectra, whose wavelengths were calibrated using the Thorium-Argon hollow cathode lamp. The reduced stellar spectra had a typical S/N of ≈25 per pixel at 1000 nm. Using IRD’s analysis pipeline for RV measurements (Hirano et al. 2020), we computed relative RVs for the individual stellar spectra. In doing so, both YJ and H-band spectra were analysed for all frames except the ones (three frames) in January 2022, for which we were only able to measure RVs for H-band spectra. The resulting RV measurements are given in Table A.1. The typical measurement errors are 10−15 and 7−10 m s−1 for the YJ-band and the H-band spectra, respectively.

2.4 High-angular-resolution imaging

Exoplanet host stars can have spatially close companions which are bound or line-of-sight objects. These companions can create a false-positive transit signal if, for example, they are an eclipsing binary or other variable star. But mainly, close companions cause ‘third-light’ flux contamination that leads to an underestimated planetary radius if not accounted for in the transit model (Ciardi et al. 2015). Companion stars can also cause non-detections of small planets residing with the same exoplanetary system (Lester et al. 2021). Thus, to search for close-in (bound) companions unresolved in TESS or other ground-based follow-up observations, we obtained high-resolution imaging speckle observations of LP 890-9.

LP 890-9 was observed on 20 March 2022 (UT) using the Zorro speckle instrument on the Gemini South 8-m telescope (Scott et al. 2021). Zorro provides simultaneous speckle imaging in two bands (562 and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections (e.g. Howell et al. 2016). Twenty-five sets of 1000 x 0.06-s exposures were collected for this faint star and subjected to Fourier analysis in our standard reduction pipeline (Howell et al. 2011). Figure 3 shows our final 5σ contrast curves and the 832 nm reconstructed speckle image. We find that LP 890-9 is a single star with no companion brighter than 4−5 magnitudes below that of the target star from ~0.2″ out to 1.2″. At the distance of LP 890-9 (32.3 pc), these angular limits correspond to spatial limits of 6.5−39 au.

thumbnail Fig. 3

Zorro speckle imaging 5σ contrast curves at 562 (blue) and 832 nm (red), along with the 832 nm reconstructed image.

3 Stellar properties

LP 890-9 is a relatively low-activity late-type M dwarf located at 32 pc. In this section, we describe the methodology used to determine its properties, which are given in Table 2 together with its catalogued photometric and astrometric parameters.

3.1 Spectroscopic analysis

Figure 4 displays the reduced Kast spectrum compared to late-M dwarf SDSS optical spectral templates from Bochanski et al. (2007). The best overall match is M6. We note that this is somewhat later than the index-based classification reported in Cruz & Reid (2002), and index-based classifications derived from relations defined in Lépine et al. (2003) and Riddick et al. (2007), all of which yield a classification of M5. Our template comparison rules this classification out. The spectrum shows a weak Ha emission line with an equivalent width of −1.5 ± 0.3 Å corresponding to log10(LHα/Lbol) = −4.58 ± 0.10 based on the χ factor relation of Douglas et al. (2014), making LP 890-9 a relatively low-activity M dwarf (cf. Newton et al. 2017). Metallicity indices from Lépine et al. (2007, 2013); Dhital et al. (2012); and Zhang et al. (2019) all indicate a dwarf metallicity classification, with a modest subsolar metallicity ([Fe/H] ≈ −0.17 to −0.12 based on Mann et al. 2013).

We performed a parallel set of analyses on the infrared SpeX spectra of LP 890-9 using the SpeX Prism Library Analysis Toolkit (SPLAT, Burgasser & Splat Development Team 2017). Results from both observations were identical. We obtained a spectral classification of M6.0 ± 0.5 by comparison to spectral standards (Kirkpatrick et al. 2010), as shown in Fig. 5. With the SpeX spectrum, we can also obtain a metallicity estimate. We used SPLAT to measure the equivalent widths of the K-band Na I and Ca I doublets and the H2O–K2 index (Rojas-Ayala et al. 2012). We then estimated the stellar metallicity using the Mann et al. (2014) relation between these observables and [Fe/H]. We used a Monte Carlo approach to calculate the uncertainty in our estimate, drawing 106 samples from normal distributions given by the means and standard deviations of the measurements. We calculated the mean and standard deviation of the resulting values and, adding in quadrature the systematic uncertainty of the relation (0.07), we arrived at our final metallicity estimate of [Fe/H] = −0.028 ± 0.089 based on the August data, consistent at 1σ with our estimate from the optical spectrum. The November data gave identical results within the uncertainties.

Table 2

Properties of the host star.

3.2 SED fitting and evolutionary modelling

We performed an analysis of the broadband spectral energy distribution (SED) of the star together with the Gaia EDR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016) and Stassun et al. (2017, 2018a). We pulled the JHKS magnitudes from 2MASS (Skrutskie et al. 2006), the W1W2W3 magnitudes from WISE (Cutri et al. 2021), and the griz magnitudes from Pan-STARRS (Chambers et al. 2016). Together, the available photometry spans the full stellar SED over the wavelength range 0.4−10 µm (see Fig. 6).

We performed a fit using NextGen stellar atmosphere models, with the free parameters being the effective temperature (Teff) and metallicity ([Fe/H]). The remaining free parameter is the extinction AV, which we fixed at zero due to the star’s proximity. The resulting fit (Fig. 6) has a reduced χ2 of 1.7, with best fit Teff = 2850 ± 75 K and [Fe/H] = 0.0 ± 0.5. Integrating the model SED gives the bolometric flux at Earth, Fbol = 4.41 ± 0.15 × 10−11 erg s−1 cm−2. Taking the Fbol and Teff, together with the Gaia EDR3 parallax, gives the stellar radius, R★ = 0.1556 ± 0.0086 R⊙.

For the mass, we applied stellar evolutionary modelling, using the models for very low-mass stars presented in Fernandes et al. (2019). We used as constraints the luminosity derived from Fbol and the Gaia EDR3 parallax (L★ = 1.4406 ± 0.0375 × 10−3L⊙), the metallicity inferred in Sect. 3.1, and assuming an age ≳2 Gyr (see Sect. 3.3)6, we obtained a stellar mass of 0.118 ± 0.002 M⊙. This uncertainty reflects 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). Considering the stellar radius estimate from SED fitting, this gives a mean stellar density of g cm−3 .

thumbnail Fig. 4

Lick/Kast red optical spectrum of LP 890-9 (solid black lines) compared to M5 (dashed magenta line), M6 (solid magenta line), and M7 (dotted magenta line) optical spectral templates from Bochanski et al. (2007). The three spectral templates are offset by a constant to facilitate comparison. Key atomic and molecular spectral features are labelled, and the inset box shows a close-up of the region around the 6563 Å Ha emission line and 6708 Å Li I absorption line (not detected).
thumbnail Fig. 5

SpeX spectrum of LP 890-9 (blue) obtained in August 2021. The SpeX Prism spectrum of the M6.0 standard LHS 1375 (Kirkpatrick et al. 2010) is shown in grey for comparison. Prominent spectral features are highlighted, and regions of strong telluric absorption are shaded. The bottom line gives the measurement uncertainty.
thumbnail Fig. 6

Spectral energy distribution (SED) of LP 890-9. The red symbols represent the observed photometric measurements (the horizontal bars represent the effective width of the bandpass). The blue symbols are the model fluxes from the best-fit NextGen atmosphere model (black curve).

3.3 Estimated stellar age

The age of LP 890-9 was estimated by comparing its UVW kinematics and metallicity to local stars, an adaptation of the method used to age-date TRAPPIST-1 (Burgasser & Mamajek 2017). The comparison sample was drawn from the GALAH Data Release 3 catalogue (Buder et al. 2021), for which ages have been estimated using the Bayesian Stellar Parameter Estimation (BSTEP) code by Sharma et al. (2018). The match set assumes individual UVW velocities within 10 km s−1 of LP 890-9 and within 1σ of the metallicity. Figure 7 compares the distribution of ages in the full GALAH sample and the match set, which shows a broad peak at Gyr (i.e., 4–9 Gyr). This age is broadly consistent with the kinematics of the source, which are consistent with the thin disk Galactic population (8% probability thick disk) based on Bensby et al. (2003).

4 Planet validation

4.1 TESS data validation report

As an initial step for false-positive vetting and before obtaining the ground-based follow-up observations described in Sect. 2, we first closely examined the TESS Data Validation Report (Twicken et al. 2018; Li et al. 2019) combining all four sectors (4, 5, 31, and 32) provided by the SPOC pipeline. TOI-4306.01 successfully passed the eclipsing binary discrimination tests: the odd and even transit depths agreed to within 0.06σ, and no significant secondary eclipse was detected. TOI-4306.01 also passed the difference image centroid offset test, which did not reveal any significant offset between the transit source and the target (3.02 ± 3.41 arcseconds), that would otherwise suggest a blend of multiple stars (e.g., background star or stellar system) and possible source confusion. In addition, the statistical bootstrap test estimated at only 2.6 × 10−17 the probability that the signal is a false alarm due to noise fluctuations in the light curve (e.g., stellar variability or residual instrumental systematics).

We note that TOI-4306.01 failed the ghost diagnostic test, which involves correlating flux time series derived from photometric core and halo aperture pixels against the transit model. A higher correlation for the halo compared to the core may indicate that the transit signal does not originate from the target, but is caused instead by scattered light or a background object, such as a background eclipsing binary. For TOI-4306.01, the test returned a slightly higher correlation statistic7 for the halo (5.17) than the core (4.72). However, the recovery of the transit signal with ground-based telescopes (see Sect. 4.2) excludes the possibility that it was caused by scattered light or other instrumental artefacts in the TESS data. Furthermore, archival imaging (see Sect. 4.3) allows us to exclude any background object as the source of TOI-4306.01’s transits.

thumbnail Fig. 7

Distribution of ages in the full GALAH DR3 sample (light grey histogram) and sources that match the UVW and metallicity of LP 890-9 to within 10 km s−1 and 1σ, respectively (dark grey histogram). The latter has a broad peak at Gyr (yellow line with uncertainties indicated by the green shaded region), which we take as an estimate for the age of LP 890-9.

4.2 Follow-up photometry

Owing to TESS’s large pixel scale of 21″ per pixel, it is not uncommon that the photometric apertures used by SPOC to extract the light curves also include some flux from contaminating stars. In this context, ground-based follow-up photometry at a higher spatial resolution is important to explore the possible contamination from nearby sources and confirm that the target star is the source of the transits. Figure 1 shows that LP 890-9 is relatively isolated. SPOC reports only a small amount of contamination from outside sources, ranging from 2.2 to 5.3% of the total flux depending on the sector and aperture used (the PDC-SAP photometry is corrected for this contamination). The closest known Gaia DR2 source (Gaia Collaboration 2018), labelled ‘2’ in Fig. 1, is 43.24″ away and about 3 magnitudes fainter (TESS-mag = 17.02, Stassun et al. 2018b). LP 890-9 is clearly resolved in the SPECULOOS-South images, which have a pixel scale of 0.35″ per pixel, thus allowing the extraction of light curves with photometric apertures of only a few arcseconds (see Sect. 2.2.1). The smallest aperture that was tested was 1.4″ (4 pixels). No transit event was observed on any nearby star, while transits with a depth matching that of TOI-4306.01 were clearly detected on the target star at the predicted times.

Our follow-up transit photometry was obtained in seven different bandpasses − g′, r′, i′, blue-blocking, I + z, z′, and 0.851.55 µm – covering a wavelength range from ~0.4 to 1.6 µm Together with the TESS transit photometry, this allowed us to check for a wavelength dependence of the transit depth, that could indicate a blend with a stellar eclipsing binary either in the background/foreground or bound to the target star. We found that the transit depths in the individual bandpasses are in very good agreement and do not show any chromatic dependence (besides the effect of stellar limb-darkening) within the accuracy of the measurements (see Sect. 5.1.2 and Table 4).

4.3 Archival imaging

The high proper motion of LP 890-9 (333 mas yr−1, Gaia Collaboration 2021) makes it possible to investigate archival images for a possible background object that could be blended with the target at its current position (see, e.g., Bakos et al. 2006). We inspected POSS I/DSS (Minkowski & Abell 1963), POSS II/DSS2 (Reid et al. 1991), and Pan-STARRS1 (Chambers et al. 2016) images spanning 64 yr with our recent observations. None of these archival images show any source at the present-day target location (see Fig. 8). Given the detection limits of these images, we can thus rule out background sources brighter than ~5 magnitudes below that of LP 890-9.

4.4 Statistical validation

To fully vet the two planet candidates, we used the open-source software package TRICERATOPS8 (Giacalone & Dressing 2020; Giacalone et al. 2021), which uses a Bayesian framework that incorporates prior knowledge of the target star, planet occurrence rates, and stellar multiplicity to calculate the probability that the transit signal is due to a planet transit or another astro-physical source. To consider a planet as statistically validated, we require its False Positive Probability (FPP), which is the sum of probabilities for all false positive scenarios, to be less than 0.01 (1%). This condition was defined by Giacalone et al. (2021) based on results obtained with TRICERATOPS for 68 TOIs (TESS Objects of Interest) that were previously identified as either confirmed planets or astrophysical false positives by TFOP using follow-up observations.

Using the phase-folded SSO I + z transit light curve (which provides the tightest photometric constraints) and the Zorro speckle imaging 5σ contrast curve at 832 nm as inputs to TRICERATOPS, we found FPPs of 6.5 × 10−6 and 3.2 × 10−7 for the inner and outer candidates, respectively. These values were obtained for each planet individually, not accounting for the fact that they are both transiting in the same system. Applying a ‘multiplicity boost’ to account for the a priori higher likelihood of a candidate in a multi-transiting system to be a real planet would allow to further reduce these FPPs, by a factor of ~20–50 (Lissauer et al. 2012; Guerrero et al. 2021). In any case, these FPPs are well below 1%, so we consider that the two transit signals correspond to validated planets. Henceforth, we denote the inner one (TOI-4306.01) as LP 890-9 b and the second outer transiting planet (SPECULOOS-2 c) as LP 890-9 c.

thumbnail Fig. 8

Imaging of the present-day position of LP 890-9 over 64 yr, to assess for a current blend with an unresolved background object. Archival images from POSS I (1957), POSS II (1994), and Pan-STARRS (2014) are presented, along with a stack image from one night of observation with SSO/Europa (2021). For each image, the position of LP 890-9 during the SSO/Europa observations is shown as a white circle.

5 Photometric analysis

5.1 Transit analysis

5.1.1 Derivation of the system parameters

We performed a joint fit of all the transit light curves described in Sect. 2 using the most recent version of the adaptive MCMC code presented in Gillon et al. (2012, 2014). This implementation makes uses of the Metropolis-Hastings algorithm (Metropolis et al. 1953; Hastings 1970) and the Gibbs sampler (Geman & Geman 1984). The data were modelled using the quadratic limb-darkening transit model of Mandel & Agol (2002) multiplied by a photometric baseline model, different for each light curve, aimed at representing the photometric variations caused by other astrophysical, instrumental, or environmental effects. Fitting the transit signals simultaneously with the possibly underlying correlated noise (instead of pre-detrending the data) ensures a better propagation of the uncertainties to the derived system parameters of interest. For each light curve, we explored a large range of baseline models consisting of low-order polynomials with respect to, e.g., time, airmass, full width at half maximum of the point spread function, target’s location on the detector, background, or any combination of these parameters. The minimal baseline model is a simple constant to account for any out-of-transit flux offset. We show in Table B.1 the baseline model selected for each light curve based on the Bayesian Information Criterion (BIC, Schwarz 1978). This table also gives for each light curve the two scaling factors, βw and βr, that were applied to the photometric error bars to account respectively for over- or under-estimated white noise and the presence of residual correlated (red) noise (see Gillon et al. 2012 for details).

The transit model parameters sampled by the MCMC for each of the two planets were: the log of the orbital period (log P), the mid-transit time (T0), the transit depth ( where Rp is the radius of the planet and R★ the stellar radius), the cosine of the orbital inclination (cos ip), and the two parameters cos ωp and sin ωp (with ep the eccentricity and ωp the argument of periastron). We also fitted the log of the stellar density (log ρ★), the log of the stellar mass (log M★), and the effective temperature (Teff). Finally, we fitted for each bandpass the combinations q1 = (u1 + u2)2 and q2 = 0.5 u1(u1 + u2)−1 of the quadratic limb-darkening coefficients (u1 and u2), following the triangular sampling scheme advocated by Kipping (2013). The coefficients of the baseline models were not sampled by the MCMC, but determined from the residuals at each step of the procedure using a singular value decomposition method (Press et al. 1992).

The priors used in our analysis are listed in the last column of Table 3. We assumed normal prior distributions for M★ and the luminosity L★ based on the values derived in Sect. 3.2 (Table 2). Normal priors for the quadratic limb-darkening coefficients (u1 and u2) were computed for each bandpass using the PyLDTk code (Parviainen & Aigrain 2015) and the library of PHOENIX high-resolution synthetic spectra of Husser et al. (2013). We increased the widths of these computed priors by a conservative factor of ten to account for model-dependent uncertainties.

We performed two different analyses, assuming either eccentric or circular orbits for the two planets (setting in this case their respective cos ωp and sin ωp values to zero). For each analysis, we ran two chains of 250 000 steps (including 20% burn-in) and checked their convergence by using the statistical test of Gelman & Rubin (1992), ensuring that the test values for all sampled parameters were <1.01. The physical parameters of the system were deduced from the above transit model parameters at each step of the MCMC. For each planet, the semi-major axis in stellar radii (ap/R★) was computed from ρ★ and P using Kepler’s third law (Seager & Mallén-Ornelas 2003). R★ was deduced from ρ★ and M★. L★ was computed from R★ and Teff. For each planet, Rp, ip, ap, ep, ωp, the stellar irradiation (Sp), and equilibrium temperature (Teq) were then easily obtained from the above parameters. We did not find evidence for eccentric orbits for any of the two planets based on the transit photometry, so we adopted the results of the circular analysis as our nominal solution (this is also consistent with simulations of the system’s tidal evolution, see Sect. 7.3).

The posterior distribution functions for the main fitted transit parameters are displayed in Fig. C.1. Table 3 gives the medians and 1σ credible intervals of the posterior distribution functions obtained for the system parameters. The corresponding best-fit transit models are shown in Fig. 9 (planet b) and 10 (planet c), together with the detrended phase-folded photometry. Most notably, we find that the two planets have similar radii of and R⊕, respectively, and that they orbit rather close to the second-order 3:1 mean motion resonance, but not exactly within (Pc/Pb = 3.098). We also note that the stellar density derived from our global transit fit is in very good agreement with the value of expected based on the stellar radius from our SED fit and the stellar mass from our evolutionary modelling (see Table 2).

Table 3

Properties of the LP 890-9 planetary system based on our global transit analysis (see Sect. 5.1.1).

thumbnail Fig. 9

Phase-folded detrended transit photometry of LP 890-9 b in each observed bandpass. The unbinned data points are shown in grey, while the black circles with error bars correspond to 10-min bins. The best-fit transit model is shown in red, with a different limb darkening for each bandpass.

5.1.2 Check for transit chromaticity

To assess the transit chromaticity, we repeated the same global analysis as in Sect. 5.1.1, but allowing this time some transit depth variations between the different bandpasses. The transit depths found for the two planets in each observed bandpass are given in Table 4. For each planet, the transit depths are all in agreement (at the ~1.1σ level) and do not show any chromatic dependence.

5.1.3 Search for transit timing variations

As mentioned in Sect. 5.1.1, the system is within 10% of the second-order 3:1 mean motion resonance. To explore our photometric dataset for possible transit timing variations (TTVs), we repeated the same global analysis as in Sect. 5.1.1, but allowing this time for each transit of the two planets a timing offset with respect to the reference transit ephemerides defined by the P and T0 values reported in Table 3. We limited our search for TTVs to the ground-based transit photometry, as TESS individual transits do not have a high enough S/N to allow precise measurements of their timings. Table 5 presents the timings and corresponding TTVs that we obtained for the 13 transits of planet b and 4 transits of planet c that were observed from the ground. Figure 11 shows the TTVs as a function of time. We do not detect any significant TTVs in the current dataset, nor any trend that could suggest the existence of additional planets in the system.

thumbnail Fig. 10

Phase-folded detrended transit photometry of LP 890-9 c in each observed bandpass. The unbinned data points are shown in grey, while the black circles with error bars correspond to 10-min bins. The best-fit transit model is shown in red, with a different limb darkening for each bandpass. We note that the MuSCAT3 r′- and i′- light curves show a possible hump around mid-transit which could be related to spot crossing.

5.2 Searches for additional transiting planets and detection limits

5.2.1 TESS photometry

As mentioned in Sect. 2.1, LP 890-9 was observed by TESS in four sectors. While the TESS Science Office issued an alert for TOI-4306.01 (LP 890-9 b), it did not report any other transit signal that could be related to LP 890-9 c or other planets in the system. However, due to the detection threshold (MES = 7.1σ) to issue an alert by the automatic SPOC and Quick Look (QLP) pipelines, some lower S/N transit signals may have gone unnoticed. We thus performed our own search for transit signals in the TESS data using our custom SHERLOCK9 pipeline (Pozuelos et al. 2020; Demory et al. 2020).

SHERLOCK is an end-to-end pipeline that combines six different modules that allow exploring the TESS data fast and robustly. These six modules consist of (1) downloading and preparing the light curves from their online repositories, (2) searching for planetary candidates, (3) performing a semiautomatic vetting of the interesting signals, (4) computing a statistical validation, (5) modelling the signals to refine their ephemerides, and (6) computing observational windows from ground-based observatories to trigger a follow-up campaign. By default, SHERLOCK works with the PDCSAP light curve and applies a multi-detrend approach employing the bi-weight algorithm provided in the wōtan package (Hippke et al. 2019) to optimise the transit search. This strategy allows the user to maximise the S/N and the signal detection efficiency (SDE) during the transit search, which is performed over the nominal PDCSAP light curve, jointly with the new detrended light curves, employing the transit least squares (TLS) package (Hippke & Heller 2019). TLS is optimised for detecting shallow periodic transits using an analytical transit model based on stellar parameters. The transit search is carried out in a loop; once a signal is found, it is stored and masked, and then the search keeps running until no more signals above a user-defined S/N threshold are found in the dataset. Each of these search-find-mask actions is called a ‘run’. Our experience points out that results found beyond five runs are less reliable due to the accumulated gaps in a given light curve after many mask-and-run iterations. Thus, in our search, we allowed a maximum of five runs. In addition, SHERLOCK allows one to apply a preliminary Savitzky-Golay filter (Savitzky & Golay 1964), which enhances the detection of shallow transits with the associated risk of obtaining more false-positive detections that the user needs to vet carefully. This filter is very useful when dealing with threshold-crossing events with low S/N like TOI-4306.01, which was found with a MES value of 8.3 by SPOC.

We first tried to recover this candidate by running two suites of searches, one with a pure TLS search and the other with the Savitzky-Golay filter applied previously. In both cases, we successfully recovered the candidate issued by TESS in the first run, with signal-to-noise ratios ranging from 3.2 to 3.9 and 6.1 to 7.1, respectively. The very low signal-to-noise ratios reported by the pure TLS search motivated us to keep using only the Savitzky-Golay filter strategy. Then, we performed two independent transit searches for extra planets by considering the four sectors simultaneously. (1) We focussed our search on orbital periods ranging from 0.5 to 40 days where a minimum of two transits was required to claim a detection. (2) We then focussed on longer orbital periods, ranging from 40 to 80 days, where single events could be recovered. None of these strategies yielded positive results, with all the signals found being attributable to systematics, noise, or variability.

Following Wells et al. (2021), several scenarios may explain the lack of extra signals in TESS data: (1) no other planets exist in the system; (2) if they do exist, they do not transit; or (3) they do exist, but the photometric precision of the data is not good enough to detect them, or they have periods longer than the ones explored in this study. If scenario (2) or (3) is true, extra planets might be detected by means of RV follow-up as discussed in Sect. 8.2. On the other hand, we know from the ground-based photometric follow-up that there is at least a second transiting planet in the system, LP 890-9 c. The fact that it was undetected by both the SPOC and SHERLOCK transit searches suggests that we are limited in this case by its detectability in the TESS data (scenario 3).

To explore the detection limits of the TESS data, we used our MATRIX ToolKit10 (Dévora-Pajares & Pozuelos 2022), which was used successfully in a number of previous studies (see, e.g., Wells et al. 2021; Schanche et al. 2022). MATRIX performs transit injection-and-recovery experiments by inserting synthetic planets over the PDCSAP light curve, combining all the sectors available. The user defines the ranges in the Rp–P parameter space to be examined. In our case, we explored the ranges of 0.5–3.R⊕ with steps of 0.28R⊕, and 0.5–30.0 days with steps of 0.33 days. In addition, each combination of Rp–P was explored using four different phases, that is, different values of T0. Hence, we explored a total of 3600 scenarios. For simplicity, the synthetic planets were injected assuming their impact parameters and eccentricities equal zero. Moreover, we detrended the light curves using a bi-weight filter with a window size of 1 day, which was found to be the optimal value during the SHERLOCK executions, and masked the transits corresponding to LP 890-9 b. We considered a synthetic planet to be recovered when its epoch was detected with 1 h accuracy and the recovered period was within 5% of the injected period (Pinjected). In addition, for each period found that did not match the injected period, MATRIX checked if it could correspond to the first harmonic (2 × Pinjected) or sub-harmonic (Pinjected/2). If this was the case, then the signal was also considered as recovered. It is worth noting that since we injected the synthetic signals in the PDCSAP light curve, then these signals were not affected by the PDC-SAP systematic corrections. In this respect, the detection limits that we find here should be considered as an optimistic scenario (see, e.g., Pozuelos et al. 2020; Eisner et al. 2020). As demonstrated previously during the SHERLOCK executions, using a preliminary Savitzky-Golay filter before the transit search yields higher S/N detections. We thus used the same strategy during the injection-and-recovery experiments. Finally, for each scenario, we performed up to three runs (as defined above for SHERLOCK) to try to recover the injected signal. Under these considerations, our transit search performed above via SHERLOCK, with a larger number of runs (five) and the multi-detrend approach, can be considered more efficient than the injection-and-recovery map obtained by MATRIX, which still offers a reasonable estimation of the detection limits.

The results, shown in Fig. 12 (upper panel), allow us to rule out planets whose sizes are > 1.8 R⊕ and orbital periods < 14 days, with recovery rates ranging from 100 to ~80%. The same planetary sizes but with orbital periods between 14 and 24 days have recovery rates ranging from ~80 to 40%. Beyond an orbital period of 24 days, the sensitivity decreases to ~20–0%. Smaller planets with sizes between 1.0 and 1.8 R⊕ have good recovery rates, between 100 and ~60%, for short orbital periods <3.2 days. LP 890-9 b lies in this region, with a recovery rate of ~60%. On the other hand, when exploring longer orbital periods >3.2 days, the detectability of such small planets with sizes between 1.0 and 1.8R⊕ decreases to ~50–0%. In this region lies LP 890-9 c, with a recovery rate of only ~6%. This result is consistent with the non-detection of LP 890-9 c in the TESS data during our transit search with SHERLOCK. For the whole range of periods explored, planets with Rp ≤ 1.0R⊕ can not be found, which hints at a clear limitation of the TESS data to discover transiting (sub-)Earth-sized planets in this system.

Table 4

Transit depths returned by our global data analysis for which we allowed transit depth variations between the different bandpasses (see Sect. 5.1.2).

Table 5

Transit timings returned by our global data analysis for which we allowed TTVs from the linear transit ephemerides defined by the orbital period P and mid-transit time T0 given in Table 3 (see Sect. 5.1.3).

thumbnail Fig. 11

TTVs of the two planets measured from the ground-based transit photometry (see Sect. 5.1.3).

5.2.2 SPECULOOS-South photometry

After the announcement of the TOI-4306.01 planet candidate by TESS in July 2021, we initiated a blind search for other transiting planets in the system with SSO and detected the first transit of planet c after about 2 months of daily monitoring. We stopped our follow-up in mid-January 2022, only when we were confident that the empirical (optimistic) habitable zone (HZ) of the star (Kopparapu et al. 2013, 2014, see also Sect. 8.3) had been reasonably well explored. To quantify this, we used a metric that we introduced in Sebastian et al. (2021), the effective phase coverage. This metric gives an estimation of how the phase of a hypothetical planet would be covered for a range of periods. In that regard, we computed the percentage of phase covered for each orbital period from P = 0.1 day to a given P = Pmax and took the effective phase coverage for periods ≤Pmax to be the integral over the period range. We stopped our photometric monitoring of LP 890-9 when we reached an effective coverage of at least 80% for Pmax = 25.823 days, which is the orbital period associated to the outer limit of the empirical HZ (Kopparapu et al. 2013, 2014, see also Sect. 8.3). Fig. 13 shows the evolution of the phase coverage as a function of the orbital period of a putative transiting planet based on the SSO photometry. The periodic drops in phase coverage that are visible for all integer periods are expected for ground-based observations from a unique site. The effective phase coverage for the outer edge of the empirical HZ is 86%. As a comparison, we detected the first transit of planet c when its effective phase coverage was only 57%, and observed its second transit when the effective coverage was 78%.

We ran an automatic transit search through the SSO I + z photometry to check for potential additional transit signals that may have been missed by our visual inspection of the data.

Transit-finding tools often require a rather ‘flat’ light curve to work optimally. Therefore, we first stitched together each night’s normalised I + z differential light curve, as these light curves are optimised for the atmospheric conditions of the night, combining light curves from all SSO telescopes. This ‘combined’ light curve was then cleaned for bad weather, as described in Murray et al. (2020), and frames were removed where the atmospheric conditions would significantly affect the photometry. We also removed all observations between 28 December 2θ2l and 3 January 2022 inclusive, as our autoguiding software, DONUTS (McCormac et al. 2013), did not run on these nights and the corresponding light curves were significantly more noisy. Then, to model remaining photometric variations (caused by intrinsic stellar variability and nightly atmospheric variations) without modelling transit features, we employed a Gaussian Process (GP) approach (Rasmussen & Williams 2006). We used a squared exponential kernel (Rasmussen & Williams 2006) implemented by the GEORGE package (Foreman-Mackey et al. 2014; Ambikasaran et al. 2015):

(1)

where i and j are two different data points in the light curve (taken at times ti and tj), A is the amplitude of correlation, and l is the length-scale over which correlations decay. To apply the GP regression, the light curve was binned every 10 min and sigma-clipped twice to lower the significance of any infrequent, short features. In essence, binning and sigma-clipping have the same effect; they remove short-duration structures to prevent over-fitting with the GP. The first sigma-clip was a basic 3σ clip of the entire light curve, to remove the most extreme outliers. We also included a second, ‘running’ sigma-clip to account for any photometric variability. Here, we binned the light curve into hours, longer than a typical transit duration, to produce a running median (a series of points representing the median of each hour bin) and a running standard deviation (defined similarly). We then computed the median running standard deviation, σ, and clipped the data points more than 3σ above or below the running median. Even though this target does not demonstrate significant photometric variability, this method can also be used to sigma-clip rapidly-varying stars.

For a transit-finding algorithm, we used the ASTROPY (Astropy Collaboration 2013, 2018) implementation of Box Least Squares, or BLS (Kovács et al. 2002). We considered using Transit Least Squares, or TLS (Hippke & Heller 2019); however, for very red stars, such as those in SPECULOOS’s target list, the detection efficiencies of BLS and TLS converge while TLS is more demanding in terms of computing time and power (Hippke & Heller 2019). We limited BLS to search for 10000 periods between 0.5 and 30 days, with a transit duration between 0.02 and 0.09 days. The highest significance period is flagged, and the corresponding transits masked before re-running the BLS transit search on the light curve. This automatic, iterative process is run until the S/N of the transits is less than 5. For LP 890-9, we first successfully detected planet b with a period of 2.730 days. Next, planet c was detected with a period of 8.457 days. All I + z transits from planets b and c in Table 1 were recovered. There were no additional, convincing detections with S/N ≥ 5.

To assess the detection efficiency of our transit-search pipeline and the detection limits of the data themselves, we performed injection and recovery tests on the SSO I + z light curve. We generated transits for 6000 artificial planets using PyTransit (Parviainen 2015), a transit light curve modelling package that implements the Mandel & Agol (2002) model. We used the limb-darkening coefficients u1 = 0.384 ± 0.028 and u2 = 0.303 ± 0.089 for the I + z filter from Sect. 5.1. For each of the planets injected, their parameters (radius Rp, orbital period P, and inclination ip) were drawn from the following distributions:

(2)

where u(a, b) represents a uniform distribution between a and b. imin and imax are the minimum and maximum inclinations for a transiting planet. The host mass and radius were assumed constant, using the values in Table 2. However, the inclination limits depend on the orbital period; therefore, when drawing the planetary parameters, we drew individual inclinations from a range set by each period. We only considered circular orbits (with e = 0) as we do not know the underlying eccentricity distribution. Close-in planets are also expected to experience significant circularisation of their orbits (Luger et al. 2017). The time at which the first transit was injected was also randomly drawn from a uniform distribution, ϕ ~ u(0,1), where ϕ is the phase of the period, such that the first transit was injected at ϕP from the start of observations.

We injected each of the 6000 artificial planets into the target’s differential light curve and ran the transit-search pipeline. The planets were injected before the GP-detrending to allow for the possibility that the GP could over-model and remove or warp the transit signals. If we were to inject planets after GP-detrending, this could artificially inflate the recovery results, as it is easier to detect a transit in a ‘flat’ light curve than in a light curve exhibiting time-dependent structure. From BLS, we only considered the most likely transit and did not look at other peaks in the periodogram. A planet was successfully recovered if at least one transit was detected by BLS above a S/N threshold of 5, with an epoch detected within an hour of the injected transit epoch. We did not include any requirements on the recovered period due to large gaps in the data from the day-night cycle and bad weather conditions. These gaps mean that we often recover period aliases. Additionally, when testing periods up to 30 days, we have a high chance of only detecting a single transit with no meaningful period information. As we only require one ‘real’ transit for recovery, it is possible that we can recover a planet even if we miss some of its transits, or incorrectly detect noise as transits. The lack of a period condition on planet recovery may thus result in an overestimation of our recovery fractions. Therefore, this injection-recovery framework is most useful for detecting individual transit signals (including single transits) which would then be followed up with manual vetting.

For all synthetic planets with radii Rp = 0.5–3.0 R⊕ and periods P = 0.5–30 days, we recovered ~72%. The lower panel of Fig. 12 shows how the recovery rates vary in the Rp – P parameter space. Comparing these results with the ones obtained in Sect. 5.2.1 for the TESS data (upper panel of Fig. 12), we see that the SSO data are sensitive to smaller planets and/or longer orbital periods. SSO has a high detection potential for short-period (P ≤ 10 days) planets in the Earth to super-Earth regime (Rp ≥ 1 R⊕) with recovery rates ranging from 86 to 100%, and an average of ~97%. This is consistent with the detection of both the b and c planets in the SSO data, with respective recovery rates of 100 and 96%. The sensitivity is still reasonable for smaller planets with 0.75 R⊕ ≤ Rp < 1 R⊕ in the same period range (P < 10 days), with an average of ~57%. However, it drops to ~15% for even smaller planets with 0.5 R⊕ ≤ Rp < 0.75 R⊕. Due to the long observation span of LP 890-9 with the SSO and our decision not to include a period criterion in our recovery, we obtain good recovery rates (above ~60%) for Earth-sized planets and super-Earths (Rp ≥ 1 R⊕) with periods longer than 10 days, up to 30 days. For the same range of periods, the sensitivity drops to ~25–50% for planets with 0.75 R⊕ ≤ Rp < 1 R⊕ and below 15% for the smallest planets with 0.5 R⊕ ≤ Rp < 0.75 R⊕. Based on these results, we conclude that it is unlikely there are any additional transiting short-period (super-)Earth-sized planets hidden in this system. However, it is still possible that there could be more transiting bodies, either too small or too far away from their host to be detected in the SSO data.

thumbnail Fig. 12

Results from the transit injection-and-recovery tests performed on the TESS (upper panel) and SSO (lower panel) data to assess the detectability of transiting planets in the LP 890-9 system. Upper panel: we explored a total of 3600 different scenarios. Each pixel evaluated 32 scenarios, that is, 32 light curves with injected planets having different P, Rp, and T0. Larger recovery rates are presented in yellow and green, while lower recovery rates are shown in blue and darker hues. LP 890-9 b (recovery rate ~60%) and c (recovery rate ~6%) are displayed as red and blue stars, respectively (see Sect. 5.2.1 for details). Lower panel: we injected 6000 artificial planets into the combined I + z-filter SSO light curve. On average, each box in this plot contains 50 injection scenarios with differing P, Rp, ip, and ϕ. Similarly to the upper panel, we include LP 890-9 b (recovery rate of 100%) and c (recovery rate of 96%) as red and blue stars, respectively (see Sect. 5.2.2 for details).
thumbnail Fig. 13

Phase coverage as a function of the orbital period using all photometric data obtained with SSO for LP 890-9. The effective coverage of 86% is the integral of the black curve. The orbital periods of planets b and c are shown as orange and red lines, respectively. The figure also shows the empirical (optimistic) HZ in green, i.e. the region between the ‘Recent Venus’ and ‘Early Mars’ limits (Kopparapu et al. 2013, 2014, and Sect. 8.3).
thumbnail Fig. 14

Global I + z-filter light curves for the SSO/Io (red) and SSO/Europa (blue) observations. The black points correspond to 20-min bins.

5.3 Stellar variability

We searched for flares and rotational variability in the SSO I + z-fìlter photometry (which has a higher precision than the TESS photometry for this late-type target) using the techniques described in Murray et al. (2022). The light curves used in this section were cleaned for bad weather and poor photometric conditions, as described in Sect. 5.2.2. However, unlike in Sect. 5.2.2, we did not use the nightly normalised light curves. Instead, we used the raw flux values from aperture photometry across all nights observed with the same telescope as one time series, and performed our differential photometry process on the entire time series, to obtain what we call a ‘global’ light curve (see also Murray et al. 2020). Doing so allows to preserve the flux relationships between nights and study long-term photometric trends, but limits our ability to optimise the photometry night-by-night (e.g., aperture size, comparison stars). The global light curves for SSO/Io and SSO/Europa are presented in Fig. 14.

The search for flares was done in two parts: a simple, automatic flare-detection algorithm was run first to extract all the flare candidates (based on the method described in Lienhard et al. 2020), followed by a manual vetting process to confirm them. We detected no conclusive flares (flare structures that could not be attributed to atmospheric effects) in any of the light curves using this method.

To search for a rotation period, we applied the Lomb-Scargle periodogram analysis (Lomb 1976; Scargle 1982) to the global light curves, binned every 20 min. The SSO light curves are not uniformly sampled as there are gaps, not only from the day/night cycle, but also from bad weather and changes to SSO’s observation strategy. Therefore, careful treatment of the resulting periodograms is essential to remove aliases. We searched for periods between the Nyquist limit (twice the bin size) and the entire observation window. We removed all peaks in the periodogram with a false alarm probability below 3σ (0.0027). We visually inspected the periodograms of LP 890-9, shown in Fig. 15, and the corresponding phase-folded light curves for all possible rotation periods (i.e., all significant peaks in this periodogram). In addition, by comparing with periodograms of the time stamps, airmass, and FWHM, we eliminated signals which arose from the non-uniform sampling and from ground-based systematics. The most likely rotation period detected for the SSO/Europa light curve (individually and combined with SSO/Io) was found to be around 50-55 days. However, as this value is approximately a third of the entire observation span of SSO/Europa, more observations would be required to confirm. The Lomb-Scargle results for SSO/Io were less clear, with the most significant peaks at 23.0, 54.0, and 4.3 days. The signal around 54 days could correspond to the one detected for SSO/Europa, but the other two are not present in SSO/Europa’s periodogram. These results were consistent regardless of whether the transits of planets b and c were masked. However, when applying the Lomb-Scargle analysis to the light curves with a longer binning of 60 min, we found the same result for SSO/Europa, but no significant peaks in the periodogram for SSO/Io. Additionally, there were unexplained jumps in the light curve on the nights of 3 December 2021 and 15 January 2022 that did not correlate with FWHM, sky background level, target’s position on the CCD, airmass, nor with the photometric aperture or any instrumental parameter (CCD temperature, exposure time, …). If these jumps are due to unknown systematics, then they will affect the reliability of our rotation analysis.

thumbnail Fig. 15

Lomb-Scargle periodograms for SSO/Io (top, red) and SSO/Europa (bottom, blue). The five highest peaks are marked by black crosses in both cases. For SSO/Europa, the most significant peak is around 50–55 days and the next highest four are at ~1 day and were determined to be aliases. The results for SSO/Io are less clear due to large gaps and much less data. The most significant peaks are found at 23.0, 54.0, and 4.3 days. The peak around 54 days could correspond to the one detected for SSO/Europa, but the other two are not present in SSO/Europa’s periodogram.
thumbnail Fig. 16

Subaru/IRD radial velocities obtained in the YJ- (grey squares) and H- (black circles) bands. The data are plotted with both their original error bars (thick solid lines) and the error bars enlarged by the best-fit jitter terms (thinner transparent lines). The top panel displays the data as a function of time. The solid red line shows the 2-planet model generated from the posterior median, while semi-transparent red lines show 25 models randomly drawn from the posteriors. The residuals around the median solution are shown in the second panel. The two lower panels show the data folded on the orbital periods of planets b and c, respectively, after removing the RV component from the other planet.

6 Radial velocity analysis

We analysed the Subaru/IRD data using the juliet package (Espinoza et al. 2019), which is built over radvel (Fulton et al. 2018) for the modelling of radial velocities and the dynesty (Speagle 2020) dynamic nested sampling algorithm for estimating Bayesian posteriors and evidences. We modelled the RVs using a sum of two Keplerians (one for each planet) and assumed circular orbits. We imposed normal priors on the orbital period P and mid-transit time T0 of each planet based on the results of our global transit analysis reported in Table 3. We assumed wide uniform priors between 0 and 500 m s−1 for the RV semi-amplitudes K. We treated the YJ-band and H-band RVs as independent datasets, thus fitting for each of them a different zero point (systemic RV) and extra jitter term. For each dataset, this jitter is added quadratically to the measurement errors of the data points, to account for any underestimation of the uncertainties or any excess noise not captured by the model. Figure 16 shows the results of our RV fit. The data are plotted with both their original error bars (thick solid lines) and the error bars enlarged by the best-fit jitter terms (thinner transparent lines). These jitters are large ( ms-1 for the YJ-band and m s−1 for the H-band) and the datasets are rather limited, so that the RV semi-amplitudes of the two planets are poorly constrained and we were only able to derive upper limits: K < 25.1 m s−1 for planet b and K < 33.0 m s−1 for planet c (2σ upper limits). The derived upper limits on the planetary masses (using the orbital inclinations from our transit analysis, Table 3) are Mp < 13.2 M⊕ for planet b and Mp < 25.3 M⊕ for planet c at 2σ. Both are well within the planetary regime.

The large extra jitters returned by our 2-planet fit suggest that there is some variability in the data that is not captured by the model. One possible explanation is that there may be an additional planet in the system which produces a detectable RV signal but was not detected in our photometry (either because it is not transiting or because it is beyond the detection limits of our photometric dataset). We thus tried adding a third planet as a third Keplerian in the model, using a wide uniform prior between 1 and 125 days (the time baseline of the RV data) for its orbital period. Comparing the Bayesian evidences Z of the two models, we found that the 3-planet model is only weakly disfavoured compared to the 2-planet one (∆ ln Z = −0.65). The orbital period found for the third Keplerian is days, with a 2σ upper limit on its semi-amplitude of 49 m s−1. Instead of adding a third Keplerian in the model, we also tried including a linear trend that could be created by a third planet with an orbital period longer than the time baseline of the RV data. We found this model to be marginally disfavoured (∆ ln Z = −3.72) compared to the 2-planet model. Since the 2-planet model is the simplest of the three models we tested and it has the highest Bayesian evidence, it appears to be the best model given the data at hand. However, additional planets in the system might be revealed with more RV data. The variability seen in the current RV dataset may also be (at least partly) related to stellar activity. In this context, it is worth noting that the period of days found for the third Keplerian in the above 3-planet fit is within the range of the possible rotation period of 50-55 days suggested by the SSO photometry (Sect. 5.3). Finally, it is also possible that the RV measurements are impacted by systematic effects. Indeed, large and currently unexplained RV variabilities have been found in Subaru/IRD data obtained for other faint and apparently quiet M dwarfs (see, e.g., Morello et al. 2022; Mori et al. 2022).

Using the probabilistic mass-radius relationship of Chen & Kipping (2017), we find predicted masses of and M⊕ for LP 890-9 b and c, respectively. Combined with the orbital parameters of the two planets, these masses would correspond to respective radial-velocity semi-amplitudes of and m s−1, significantly smaller than the measurement errors and variability of the current Subaru/IRD data. Additional RV measurements, possibly obtained with other available or next-generation instruments, are required to improve the preliminary mass constraints derived here and investigate the possible presence of other planets in the system. We discuss some prospects for RV follow-up in Sect. 8.2.

7 Dynamical analysis

7.1 Expected TTVs

As described in Sect. 5.1.3, we do not detect any significant TTVs in the current photometric dataset, although the two planets are close to the second-order 3:1 mean motion resonance (but not exactly within: Pc/Pb = 3.098). In general, the lack of observed TTVs might be caused by either their low amplitude or the poor baseline coverage provided by the observational data. In our case, the ground-based observations that we used to explore TTVs were obtained over ~150 days, which is much longer than the super-period of 2.64 days (Lithwick et al. 2012) for LP 8909 b and c. Hence, our ground-based data should show any TTVs if they were large enough. To test the expected TTV amplitudes in the system, we used the TTVFast2Furious package (Hadden 2019) following the strategy presented by Cloutier et al. (2020). We ran 103 realisations with orbital periods and times of mid-transit sampled from the posterior distributions given in Table 3. The planetary masses were drawn from their predictions from the probabilistic mass-radius relationship of Chen & Kipping (2017; see Sect. 6) and the argument of periastron following a random distribution from 0 to 2π. We found maximum TTV values for LP 890-9 b and c of only 1.5 and 1.0 s, respectively. As a comparison, the typical timing uncertainties returned by our global transit analysis are on the order of 1 minute (see Table 5). Based on the low expected TTV amplitudes, we conclude that the LP 890-9 system is not suitable for a determination of the planetary masses (and dynamical studies in general) using the TTV method. This is in contrast with the TRAPPIST-1 system, where the outer five planets are close to first-order11 resonances with adjacent planets, resulting in large-amplitude TTVs (from ~5 to 100 min) – an important feature which made it possible to derive precise planetary masses and eccentricities (Agol et al. 2021).

7.2 Stability

From the current dataset, we concluded that the LP 890-9 system is composed of at least two planets orbiting in nearly circular orbits and with similar planetary sizes. To assess the global stability of the system, we used the Mean Exponential Growth factor of Nearby Orbits (MEGNO; Cincotta & Simó 1999, 2000; Cincotta et al. 2003), Y(t) parameter. This chaos index is widely used to explore the stability of extrasolar planetary systems (e.g., Hinse et al. 2015; Jenkins et al. 2019; Delrez et al. 2021). In short, its time-averaged mean value, 〈Y(t)〉, amplifies any stochastic behaviour, which can be used to distinguish between quasi-periodic orbits, if 〈Y(t → ∞)〉 → 2, or chaotic trajectories, if 〈Y(t → ∞)〉 → ∞. We made use of the MEGNO implementation with the N-body integrator REBOUND (Rein & Liu 2012), which employs the Wisdom-Holman WHfast code (Rein & Tamayo 2015). In this context, we studied a set of scenarios, where the planetary properties of LP 890-9 b and c were drawn from the posterior distributions and derived parameters given in Table 3. In particular, semi-major axes, orbital inclinations, and planetary masses were sampled from uniform distributions using their nominal values up to 3σ of their uncertainties. In addition, the mean anomalies and arguments of periastron were sampled randomly from 0 to 2π. We ran two suites of simulations: first, we considered circular orbits, and secondly, slightly eccentric orbits using a uniform distribution of 0.10 ± 0.05 for both planets. In each experiment, we generated 1000 scenarios, and they were integrated over 106 orbits of the outer planet, with a time-step of 5% of the orbital period of the inner planet. We found that the LP 890-9 system is highly stable in both experiments with 〈Y(t)〉 clustering between 1.999–2.001 in 99.3 and 93.8% of cases for circular and eccentric orbits, respectively.

Having shown the stability of the LP 890-9 system on its nominal configuration and within 3σ of the planetary properties, we considered possible stability limitations for the planetary masses and the orbital eccentricities, the two main parameters that have a major impact on the orbital dynamics and that are poorly constrained by the current data. To this end, we extended our analyses by building stability maps following Delrez et al. (2021), that is, by exploring the Mb–Mc and eb–ec parameter spaces. We explored planetary masses ranging from 1 to 20 M⊕ and eccentricities ranging from 0.0 to 0.6 for both planets. In these sets of simulations, we took 20 values from each range, so each stability map has 20 × 20 pixels. Moreover, we set 10 different initial conditions for each realisation, randomly sampling the argument of periastron and the mean anomaly from 0 to 2π. Hence, each stability map consisted of 4000 scenarios. We then averaged the results of the 10 different initial conditions to obtain the averaged value of 〈Y〉 for each pixel of the stability map. The integration time and time-step were the same as previously used. We found that the system is fully stable over the whole range of masses explored with ∆〈Y(t)〉= 2.0 − 〈Y(t)〉 < 10−5. On the other hand, we found that the eccentricities have to satisfy the condition eb ≤ 0.39 − 1.11 × ec. This implies that the maximum eccentricity for planet b (when planet c has a circular orbit) is 0.39, and the maximum eccentricity for planet c (when planet b has a circular orbit) is 0.35. For eccentricity values that do not meet this condition, the system may become unstable. By combining the results obtained along this subsection and assuming that there are no more planets in the system, we conclude that the dynamical architecture of the LP 890-9 system resides in a stable set of parameters, allowing for long-term evolution.

7.3 Tidal evolution

The two planets of the LP 890-9 system reside in close-in orbits and, consequently, due to the short distance from the host star, are affected by tides. Indeed, tidal interactions between the star and planets in closely-packed systems have an impact on planetary orbits, altering the obliquities (e) and rotational periods (Pot) on short time scales, driving the planets into a particular configuration called tidal locking, where the tidal torques fix the rotation rate to a specific frequency, and the planetary and stellar spins are aligned. On longer time scales, the tidal interactions can also induce the circularisation of the orbits and reduce their semi-major axes (see, e.g., Bolmont et al. 2014; Barnes 2017). This scenario is interesting since a tidally-locked planet in a circular orbit rotates synchronously (i.e., its rotational and orbital periods are equal), meaning that the same hemisphere of the planet is always facing the star, a configuration that can have important implications for the planet’s global climate.

To explore the evolution of the LP 890-9 system under tidal interactions, we employed the constant time-lag model (CTL), where a planet is modelled as a weakly viscous fluid deformed by the gravitational interactions (see, e.g., Mignard 1979; Hut 1981; Eggleton et al. 1998; Leconte et al. 2010). We used the implementation of CTL in the POSIDONIUS N-body code (Blanco-Cuaresma & Bolmont 2017). The main factor in the CTL model that controls the tidal dissipation of a given planet is the product of two free parameters: the degree-2 potential Love number, k2, and the constant time-lag, ∆τ. The k2 accounts for the elastic properties and the self-gravity of the planets, with values ranging from 0 to 1.5. On the other hand, ∆τ refers to the lag between the line connecting the two centers of mass (star-planet) and the direction of the tidal bulges, which can span orders of magnitudes (Barnes 2017). Although the masses (and thus bulk densities) of the LP 890-9 planets are currently poorly constrained (Sect. 6), their small sizes imply that they are most likely terrestrial. Following Turbet et al. (2018) and Bolmont et al. (2020), who studied the tidal effects for the TRAPPIST-1 planets, we thus used Earth’s value k2,⊕∆τ⊕ = 213 s (Neron de Surgy & Laskar 1997) as a reference and explored dissipation factors ranging from 0.1 × fc2,θ∆τφ (for dry rocky planets such as Mars) to 10 × k2,⊕∆τ⊕ (for planets dominated by fluid-solid friction), to evaluate a range of possible tidal behaviours.

In a first step, we assessed the time needed by the LP 890-9 planets to evolve towards the tidal locking configuration. To this end, we conducted a suite of simulations considering different scenarios combining three different initial planetary rotational periods P0,rot (10,100, and 1000 h), three different initial obliquities ϵ0 (15°, 45°, and 75°) and three dissipation factors (0.1×, 1×, and 10×k2,⊕∆τ⊕). We found that the slowest evolution towards tidal locking happened for the particular combination of initial conditions given by ϵ0 = 75° and P0,rot = 10 h for dry terrestrial planets (0.1 × k2,⊕∆τ⊕). In such a case, we found time scales of ~1.5 × 104 and ~1.5 × 106 yr, for planets b and c, respectively. Hence, these values might be considered as upper limits on the time scale required for each planet to evolve into a tidally-locked configuration. Considering the estimated age of the system (Sect. 3.3), the two planets are likely tidally locked.

In a second step, we evaluated the time needed for the circularisation of the orbits due to tidal interactions. From the global transit analysis presented in Sect. 5.1.1, we did not find any evidence for eccentric orbits. However, the eccentricities of the planets are poorly constrained by the transit photometry, and high-precision RV measurements are needed for a proper estimation. On the other hand, when studying the dynamical architecture of the LP 890-9 system in Sect. 7.2, we found that it may tolerate a certain level of planetary eccentricities allowing for long-term dynamical evolution. Then, assuming initial slightly eccentric orbits of 0.01 for both planets, we performed three simulations considering dissipation factors of 0.1×, 1×, and 10 × k2,⊕∆τ⊕. We found that the eccentricities of both planets decreased to reach mean equilibrium values of ~1.8 × 10−4 and −1.7 × 10−5 for planets b and c, respectively. The time scales were of the order of ~109, ~108 and ~107 yr, for the different dissipation factors considered, being longer for dry rocky planets. These small residual eccentricities, inferior to 10-3, are the result of the competition between tidal damping and planet-planet excitations (Turbet et al. 2018). Then, considering that the system is likely older than 2 Gyr (Sect. 3.3), planetary orbits may have had time to evolve towards these low eccentricities due to tidal interactions. Thus, should future observations reveal larger eccentricities than the equilibrium values found here, this could suggest that other unknown planets in the system might be exciting the orbits.

8 Discussion

8.1 LP 890-9 in the context of other planetary systems around cool dwarfs

Figure 17 shows the planets that have been found so far around cool dwarfs using the transit and radial velocity techniques, as a function of the stellar flux they receive and the effective temperature of their host star. Aside from LP 890-9, only five planetary systems are known around hosts with Teff < 3000 K: TRAPPIST-1 (Gillon et al. 2016, 2017), Teegarden’s star (Zechmeister et al. 2019), GJ 1061 (Dreizler et al. 2020), LP 79118 (Crossfield et al. 2019), and LHS 1140 (Dittmann et al. 2017; Ment et al. 2019). With an effective temperature of 2850 ± 75 K (Sect. 3), LP 890-9 is the second-coolest star found to host planets after TRAPPIST-1. Interestingly, all these systems – including LP 890-9 – are multi-planet systems, in line with the observation that compact multi-planetary systems should be a relatively common outcome of planet formation around low-mass stars (e.g., Coleman et al. 2019; Miguel et al. 2020; Mulders et al. 2021).

8.2 Prospects for RV follow-up

As mentioned in Sect. 6, the expected RV semi-amplitudes of the two planets based on the mass-radius relationship of Chen & Kipping (2017) are and m s−1, respectively. Given the faintness of the star, detecting such small signals will be challenging, even for state-of-the-art instrumentation, and will require a substantial investment of observing time.

Besides its late spectral type, the star is too faint to recover the expected RV semi-amplitudes with the SPIROU infrared spectrograph at the CFHT (Donati et al. 2020) or to be measured with the upcoming NIRPS infrared spectrograph at the 3.6m ESO telescope (limiting magnitude H = 9, Wildi et al. 2017). The expected RV semi-amplitudes are also too small to be recovered with CRIRES+ at the VLT (RV precision12 ~20 m s−1, Dorn et al. 2014). Simply scaling from the existing RV constraints derived from the current Subaru/IRD data (Kb < 25.1 m s−1 and Kc < 33.0 m s−1) and the expected RV semi-amplitudes given above suggests that ~31 and ~100 times more data would be needed to constrain the masses of planets b and c, respectively, with that facility.

In the optical range, currently available high-resolution spectrographs, like ESPRESSO at the VLT (Pepe et al. 2021) and MAROON-X at Gemini-North (Seifahrt et al. 2020), would also require substantial amounts of observing time to detect the RV signals of the planets. Even assuming an inactive and slow rotating star (stellar jitter ~1.5 m s−1 which is probably optimistic, see Sect. 6), we estimate that about 60 spectra with an average exposure time of 1800s would be required with MAROON-X to recover the RV signals of both planets with 4-5σ. With the future ANDES spectrograph at the ELT (Marconi et al. 2021), a 7–9σ detection should be feasible by investing only half of this observing time.

To summarise, while a 4–5σ detection of the RV signals would require a large amount of observing time with current high-resolution spectrographs, measuring the masses of the two planets with at least 9σ accuracy will require waiting for next-generation instruments to become available in the next decade. Furthermore, such an RV follow-up will also allow to investigate the origin of the large variability seen in the Subaru/IRD data (Sect. 6) and may reveal some additional planets in the system.

thumbnail Fig. 17

Exoplanets found around cool dwarfs with an effective temperature lower than 3600 K using the transit (black circles) or radial velocity (grey triangles) techniques. Data were taken from NASA Exoplanet Archive on 20 April 2022. The planets are shown as a function of their incident stellar flux and the effective temperature of their host star. The size of the symbols is proportional to the planetary radius, estimated from the mass using the probabilistic mass-radius relationship of Chen & Kipping (2017) when it is unknown (for non-transiting RV planets). LP 890-9 b and c are shown as orange circles. The inner (Recent Venus) and outer (Early Mars) boundaries of the empirical (optimistic) HZ are shown as solid red and blue lines, respectively (Kopparapu et al. 2014 and Sect. 8.3). The dashed green line shows an alternative conservative inner edge limit from 3D models (Leconte et al. 2013). The outer edge of the HZ in 3D models agrees with incident stellar flux for Early Mars insulation. We note that the three limits shown do vary with the stellar effective temperature, but this is not obvious on this plot due to the log-scale of the x-axis.

8.3 Potential habitability

The habitable zone (HZ) is a concept that can be used to prioritise detected rocky exoplanets for follow-up observations. It is defined as the circular region around one or multiple stars, where liquid water is possible on the surface of a geologically-active rocky planet (e.g., Kasting et al. 1993; Kopparapu et al. 2013; Kaltenegger & Haghighipour 2013; Haghighipour & Kaltenegger 2013; Kane & Hinkel 2013; Ramirez & Kaltenegger 2016). The width and distance of the HZ region depends in a first approximation on the incident stellar flux and the atmospheric composition of the planet.

The limits of the conservative HZ is set by the greenhouse effect of CO2 and H2O vapour. At the outer edge, condensation and scattering by CO2 outstrips its greenhouse capacity. At the inner edge, the mean surface temperatures exceed the critical point of water, triggering a runaway greenhouse environment and rapid water loss to space on short timescales (see, e.g., Kasting et al. 1993; Kopparapu et al. 2014; Kaltenegger 2017). Modelling clouds is complex especially because no data exists that can provide comparison data sets for, e.g., fast rotating or slow rotating, very hot or very cold Earth-like planets. However, both 1D and 3D models generally find wider boundaries when considering cloud feedback than the original classical HZ but models differ in their specific results (see, e.g., the review in Kaltenegger 2017). The updated conservative HZ limits (Kopparapu et al. 2014) have been scaled to the 3D models by Leconte et al. (2013) to account for some cloud feedback and are used here.

The concept of the empirical (or optimistic) HZ relies on observations in our own Solar System (Kasting et al. 1993). It is based on the incident flux on Venus and Mars when there were no more hints of liquid water on their surface. In our Solar System, that translates to 1.76 current Earth’s incident irradiation (S⊕) for Venus around 1 billion years ago and 0.32 S⊕ for Mars about 3.8 billion years ago, providing empirical proxies for the incident flux limits of habitability. The inner edge of the empirical HZ is not well known because of the geological active history of Venus’ surface. The outer edge of the HZ in 3D models agrees with incident stellar flux for Early Mars insulation. The empirical HZ provides an empirical dataset for the search for habitable planets, complementary to modelling efforts.

Cool stars warm an Earth-like planet’s surface more effectively than hotter stars (see, e.g., Kasting et al. 1993). Thus, for LP 890-9, the inner limits of the HZ are 1.49 S⊕ for the empirical and 0.92 S⊕ for the conservative HZ respectively (Kopparapu et al. 2013). The outer limit of the HZ for LP 890-9 is found at a stellar flux of about 0.21 S⊕. While the HZ limits can change with atmospheric composition (e.g., Pierrehumbert & Gaidos 2011; Ramirez & Kaltenegger 2017, 2018), here we used the empirical HZ and the conservative HZ defined for a planet with an N2-H2O-CO2 dominated atmosphere like Earth to identify the HZ limits for the LP 890-9 system.

With an incident stellar flux of 4.09 ± 0.12 S⊕, LP 890-9 b is too close to the star to be within the HZ. However, LP 890-9 c receives 0.906 ± 0.026 S⊕, which places it within the limits of both the empirical and the conservative HZ. It orbits close to the inner edge of the conservative HZ, and in the inner third of the empirical HZ of LP 890-9 (see Fig. 17). Due to increasing stellar irradiation towards the inner edge of the HZ, LP 890-9 c should develop a water-rich atmosphere because more surface water would evaporate than on Earth due to increased temperatures. Such an increase in water vapour in the atmosphere should be detectable with JWST and upcoming Extremely Large Telescopes (see Sect. 8.4).

While LP 890-9 c orbits within the HZ of its host star, due to its position close to a cool M star, a potential slow rotation due to tidal forces could maintain habitability even for higher incident fluxes than the nominal limits of the HZ, due to an increased cloud coverage on the day-side of the planet. This makes planets that could be synchronously locked at the inner edge of the HZ, like LP 890-9 c, very interesting test-cases for this hypothesis (see, e.g., Yang et al. 2013; Wolf & Toon 2014; Kopparapu et al. 2016).

thumbnail Fig. 18

Expected transmission signals (see Sect. 8.4 for details) of currently known transiting terrestrial exoplanets (Rp < 1.6 R⊕) as a function of their incident stellar flux. Data were taken from NASA Exoplanet Archive on 29 March 2022. We only show planets with mild incident irradiations, lower than five times that received by the Earth (Sp < 5 S⊕). The symbols are colour-coded as a function of their S/N relative to TRAPPIST-1 b, obtained by scaling the signal amplitude with the hosts’ brightness in the J-band and using TRAPPIST-1 b’s S/N as a reference. The size of the symbols is proportional to the planetary radius. Transmission signals less than 5 ppm have been removed to enhance the readability of the figure.

8.4 Prospects for atmospheric characterisation

Rocky exoplanets orbiting in the HZ of nearby late-type M dwarf stars provide unique opportunities for characterising their atmospheres as well as searching for biosignature gases. To assess the potential of the system for atmospheric characterisation in the context of other known transiting terrestrial planets, we followed the same approach as for K2-315b in Niraula et al. (2020). We retrieved the currently known transiting exoplanets with a reported radius below 1.6 R© from NASA Exoplanet Archive13. We focussed here on temperate planets with mild incident irradiations, lower than five times that received by the Earth (Sp < 5 S⊕), and estimated their typical signal amplitude in transit transmission spectroscopy as:

(3)

where heff is the effective atmospheric height (that is, the extent of the atmospheric annulus probed in transmission). heff is directly proportional to the planet’s atmospheric scale height H: , where k is Boltzmann’s constant, T is the atmospheric temperature, µ is the atmospheric mean molecular mass, and g is the surface gravity. NH is a number ranging typically between 6 and 10, depending on the presence of clouds, the spectral resolution and range covered (Miller-Ricci et al. 2009; de Wit & Seager 2013). We assumed here NH = 7, µ = 20 amu (relevant for scenarios of secondary atmospheres), and T to be equal to the equilibrium temperature, calculated for a sphere of null Bond albedo and an efficient heat redistribution. For planets with unknown masses, we estimated g using the probabilistic mass-radius relationship of Chen & Kipping (2017).

Figure 18 shows the estimated transmission signals as a function of the planets’ incident stellar flux. The points are colour-coded as a function of their S/N relative to TRAPPIST-1 b, obtained by scaling the signal amplitude with the hosts’ brightness in the J-band and using TRAPPIST-1 b’s S/N as a reference. LP 890-9 b compares well to the outer TRAPPIST-1 planets in terms of potential for atmospheric characterisation, while LP 890-9 c stands out as the second-most favourable habitable-zone terrestrial planet after the TRAPPIST-1 planets. The projected signals in transmission are ~125 and ~88 ppm for LP 890-9 b and c, respectively, with corresponding signal-to-noise ratios relative to TRAPPIST-1 b of 0.38 and 0.26. This means that assuming comparable secondary atmospheres, LP 890-9b and c would require respectively 7 and 15 times the amount of in-transit time to be characterised to the same extent as TRAPPIST-1 b. Assessing the presence of a 10bar CO2 atmospherewithJWST/NIRSpecwouldthus requirearound 14 and 30 transits for planets b and c, respectively (scaled from Lustig-Yaeger et al. 2019). Repeating the same exercise to compare LP 890-9 c with TRAPPIST-1 d and e, which receive closer amounts of stellar flux, we find a S/N of 0.29 relative to TRAPPIST-1 d and 0.38 relative to TRAPPIST-1 e. Characterising LP 890-9 c to the same extent as TRAPPIST-1 d (resp. TRAPPIST-1 e) would thus require 12 (resp. 7) times the amount of in-transit time. Based on these estimations, about 90 transits would be needed to detect a 1 bar H2O atmosphere around LP 890-9 c with JWST/NIRSpec (scaled again from Lustig-Yaeger et al. 2019). This rough comparison assumes a similar atmosphere for all planets. However, since LP 890-9 c is located close to the innerlimit(runaway greenhouse) ofthe conservative HZ (see Sect. 8.3), it may have a thicker H2O atmosphere which would boost its atmospheric signals.

9 Conclusions

We have presented the discovery and initial characterisation of the LP 890-9 system, which hosts two temperate super-Earths transiting a nearby M6 dwarf. LP 890-9 is the second-coolest star found to host planets after TRAPPIST-1. The seven terrestrial planets orbiting TRAPPIST-1 have garnered the interest of a broad community across very diverse scientific disciplines: planet formation and evolution (e.g., Ormel et al. 2017; Coleman et al. 2019; Raymond et al. 2022), star-planet interactions (e.g., Bolmont et al. 2017; O’Malley-James & Kaltenegger 2017; Dong et al. 2018), multi-planet dynamics (e.g., Grimm et al. 2018; Agol et al. 2021; Teyssandier et al. 2022), interior modelling (e.g., Dorn et al. 2018; Barr et al. 2018; Barth et al. 2021), atmospheric observations (e.g., de Wit et al. 2016, 2018; Bourrier et al. 2017) and modelling (e.g., Lustig-Yaeger et al. 2019; Pidhorodetska et al. 2020; Lin et al. 2021), or climate predictions (e.g., Turbet et al. 2018, 2021; Sergeev et al. 2021), among others. The discovery of the remarkable LP 890-9 system presented in this work offers another rare opportunity to study temperate terrestrial planets around our smallest and coolest neighbours.

Acknowledgements

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. 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. This paper includes data collected by the TESS mission that are publicly available from the Mikulski Archive for Space Telescopes (MAST). The research leading to these results has received funding from the European Research Council (ERC) under the FP/2007-2013 ERC grant agreement no 336480, and under the European Union’s Horizon 2020 research and innovation programme (grants agreements no 679030 & 803193/BEBOP); from an Action de Recherche Concertée (ARC) grant, financed by the Wallonia-Brussels Federation, from the Balzan Prize Foundation, from the BELSPO/BRAIN2.0 research program (PORTAL project), from the Science and Technology Facilities Council (STFC; grants no ST/S00193X/1, ST/00305/1, and ST/W000385/1), and from F.R.S-FNRS (Research Project ID T010920F). This work was also partially supported by a grant from the Simons Foundation (PI: Queloz, grant number 327127), as well as by the MERAC foundation (PI: Triaud). TRAPPIST is funded by the Belgian Fund for Scientific Research (Fond National de la Recherche Scientifique, FNRS) under the grant PDR T.0120.21, with the participation of the Swiss National Science Fundation (SNF). This work is partly supported by MEXT/JSPS KAKENHI Grant Numbers JP15H02063, JP17H04574, JP18H05439, JP18H05442, JP19K14783, JP21H00035, JP21K13975, JP21K20376, JP22000005, Grant-in-Aid for JSPS Fellows Grant Number JP20J21872, JST CREST Grant Number JPMJCR1761, the Astrobiology Center of National Institutes of Natural Sciences (NINS) (Grant Numbers AB031010, AB031014), and Social welfare juridical person SHIYUKAI (chairman MASAYUKI KAWASHIMA). This paper is based on data collected at the Subaru Telescope, which is located atop Maunakea and operated by the National Astronomical Observatory of Japan (NAOJ). We wish to recognise and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. This paper is based on observations made with the MuSCAT3 instrument, developed by the Astrobiology Center and under financial supports by JSPS KAKENHI (JP18H05439) and JST PRESTO (JPMJPR1775), at Faulkes Telescope North on Maui, HI, operated by the Las Cumbres Observatory. Some of the observations in the paper made use of the High-Resolution Imaging instrument Zorro obtained under Gemini LLP Proposal Number: GN/S-2021A-LP-105. Zorro 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. Zorro was mounted on the Gemini North (and/or South) 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). We acknowledge funding from the European Research Council under the ERC Grant Agreement n. 3 37591-ExTrA. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. This work is based upon observations carried out at the Observatorio Astronómico Nacional at the Sierra de San Pedro Mártir (OAN-SPM), Baja California, México. We warmly thank the entire technical staff of the Observatorio Astronómico Nacional at San Pedro Mártir for their unfailing support to SAINT-EX operations. Research at Lick Observatory is partially supported by a generous gift from Google. L.D. is an F.R.S.-FNRS Postdoctoral Researcher. M.G. and E.J. are F.R.S.-FNRS Senior Research Associates. V.V.G. is an F.R.S.-FNRS Research Associate. B.V.R. thanks the Heising-Simons Foundation for support. Y.G.M.C. acknowledges support from UNAM-PAPIIT IG-101321. B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2-163967 and PP00P2-190080). M.N.G. acknowledges support from the European Space Agency (ESA) as an ESA Research Fellow. A.H.M.J.T acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no 803193/BEBOP), from the MERAC foundation, and from the Science and Technology Facilities Council (STFC; grants no ST/S00193X/1, ST/00305/1, and ST/W000385/1). E.D. acknowledges support from the innovation and research Horizon 2020 program in the context of the Marie Sklodowska-Curie subvention 945298. V.K. acknowledges support from NSF award AST2009343. This publication benefits from the support of the French Community of Belgium in the context of the FRIA Doctoral Grant awarded to M.T.

Appendix A Radial velocities

Table A.1

Subaru/IRD radial velocity measurements. We recommend treating the YJ-band and H-band RVs as independent datasets, as done in the analysis presented in Sect. 6.

Appendix B Photometric baseline models and error scaling factors

Table B.1

Photometric baseline models, residual RMS, and error scaling factors βw and βr for each transit light curve used in our global transit analysis (see Sect. 5.1.1). For the baseline function, p(αn) denotes a n-order polynomial function of the parameter α, with α that can be t = time, a = airmass, x and y = x- and y- position of the target on the detector, or f = full width at half maximum of the point spread function.

Appendix C Transit fit posterior distributions

thumbnail Fig. C.1

MCMC posterior distributions for the main fitted transit parameters (see Sect. 5.1.1): the log of the orbital period (log P), the mid-transit time (T0), the transit depth (dF), and the cosine of the orbital inclination (cos ip) for each of the two planets; and the log of the stellar density (log ρ★). This plot was made using the corner python package (Foreman-Mackey 2016).

References

  1. Agol, E., Dorn, C., Grimm, S. L., et al. 2021, Planet. Sci. J., 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
  2. 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]
  3. Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38, 252 [Google Scholar]
  4. Astropy Collaboration (Robitaille, T.P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A.M., et al.) 2018, AJ, 156, 123 [NASA ADS] [CrossRef] [Google Scholar]
  6. Bailey, S. 2012, PASP, 124, 1015 [NASA ADS] [CrossRef] [Google Scholar]
  7. Bakos, G. Á., Pál, A., Latham, D. W., Noyes, R. W., & Stefanik, R. P. 2006, ApJ, 641, L57 [NASA ADS] [CrossRef] [Google Scholar]
  8. Barclay, T., Pepper, J., & Quintana, E. V. 2018, ApJS, 239, 2 [Google Scholar]
  9. Barnes, R. 2017, Celesti. Mech. Dyn. Astron., 129, 509 [NASA ADS] [CrossRef] [Google Scholar]
  10. Barr, A. C., Dobos, V., & Kiss, L. L. 2018, A&A, 613, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Barth, P., Carone, L., Barnes, R., et al. 2021, Astrobiology, 21, 1325 [NASA ADS] [CrossRef] [Google Scholar]
  12. Bensby, T., Feltzing, S., & Lundström, I. 2003, A&A, 410, 527 [CrossRef] [EDP Sciences] [Google Scholar]
  13. Blanco-Cuaresma, S., & Bolmont, E. 2017, in IAU Symp., 325, 341 [NASA ADS] [Google Scholar]
  14. Bochanski, J. J., West, A. A., Hawley, S. L., & Covey, K. R. 2007, AJ, 133, 531 [NASA ADS] [CrossRef] [Google Scholar]
  15. Bolmont, E., Raymond, S. N., von Paris, P., et al. 2014, ApJ, 793, 3 [CrossRef] [Google Scholar]
  16. Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bolmont, E., Demory, B. O., Blanco-Cuaresma, S., et al. 2020, A&A, 635, A117 [EDP Sciences] [Google Scholar]
  18. Bonfils, X., Almenara, J. M., Jocou, L., et al. 2015, SPIE Conf. Ser., 9605, 96051L [NASA ADS] [Google Scholar]
  19. Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121 [NASA ADS] [CrossRef] [Google Scholar]
  20. Buder, S., Sharma, S., Kos, J., et al. 2021, MNRAS, 506, 150 [NASA ADS] [CrossRef] [Google Scholar]
  21. Burdanov, A., Delrez, L., Gillon, M., & Jehin, E. 2018, SPECULOOS Exoplanet Search and Its Prototype on TRAPPIST, eds. H.J. Deeg, & J.A. Belmonte (Berlin: Springer), 130 [Google Scholar]
  22. Burgasser, A. J., & Mamajek, E. E. 2017, ApJ, 845, 110 [Google Scholar]
  23. Burgasser, A. J., & Splat Development Team 2017, ASI Conf. Ser., 14, 7 [NASA ADS] [Google Scholar]
  24. Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, ArXiv e-prints [arXiv:1612.05560] [Google Scholar]
  25. Chen, J., & Kipping, D. 2017, ApJ, 834, 17 [Google Scholar]
  26. Ciardi, D. R., Beichman, C. A., Horch, E. P., & Howell, S. B. 2015, ApJ, 805, 16 [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&AS, 147, 205 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Phys. D Nonlinear Phenomena, 182, 151 [NASA ADS] [CrossRef] [Google Scholar]
  30. Cloutier, R., Eastman, J. D., Rodriguez, J. E., et al. 2020, AJ, 160, 3 [Google Scholar]
  31. Cointepas, M., Almenara, J. M., Bonfils, X., et al. 2021, A&A, 650, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  32. Coleman, G. A. L., Leleu, A., Alibert, Y., & Benz, W. 2019, A&A, 631, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Crossfield, I. J. M., Waalkes, W., Newton, E. R., et al. 2019, ApJ, 883, L16 [Google Scholar]
  34. Cruz, K. L., & Reid, I. N. 2002, AJ, 123, 2828 [NASA ADS] [CrossRef] [Google Scholar]
  35. Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, PASP, 116, 362 [NASA ADS] [CrossRef] [Google Scholar]
  36. Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog: II/328 [Google Scholar]
  37. de Wit, J., & Seager, S. 2013, Science, 342, 1473 [Google Scholar]
  38. de Wit, J., Wakeford, H. R., Gillon, M., et al. 2016, Nature, 537, 69 [NASA ADS] [CrossRef] [Google Scholar]
  39. de Wit, J., Wakeford, H. R., Lewis, N. K., et al. 2018, Nat. Astron., 2, 214 [NASA ADS] [CrossRef] [Google Scholar]
  40. Delrez, L., Gillon, M., Queloz, D., et al. 2018, SPIE Conf. Ser., 10700, 107001I [Google Scholar]
  41. Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., 5, 775 [NASA ADS] [CrossRef] [Google Scholar]
  42. 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]
  43. Dévora-Pajares, M., & Pozuelos, F. J. 2022, https://doi.org/10.5281/zenodo.6570831 [Google Scholar]
  44. Dhital, S., West, A. A., Stassun, K. G., et al. 2012, AJ, 143, 67 [NASA ADS] [CrossRef] [Google Scholar]
  45. Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333 [NASA ADS] [CrossRef] [Google Scholar]
  46. Donati, J. F., Kouach, D., Moutou, C., et al. 2020, MNRAS, 498, 5684 [Google Scholar]
  47. Dong, C., Jin, M., Lingam, M., et al. 2018, Proc. Natl. Acad. Sci., 115, 260 [NASA ADS] [CrossRef] [Google Scholar]
  48. Dorn, R. J., Anglada-Escude, G., Baade, D., et al. 2014, The Messenger, 156, 7 [NASA ADS] [Google Scholar]
  49. Dorn, C., Mosegaard, K., Grimm, S. L., & Alibert, Y. 2018, ApJ, 865, 20 [NASA ADS] [CrossRef] [Google Scholar]
  50. Douglas, S. T., Agüeros, M. A., Covey, K. R., et al. 2014, ApJ, 795, 161 [NASA ADS] [CrossRef] [Google Scholar]
  51. Dreizler, S., Jeffers, S. V., Rodríguez, E., et al. 2020, MNRAS, 493, 536 [Google Scholar]
  52. Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, ApJ, 499, 853 [NASA ADS] [CrossRef] [Google Scholar]
  53. Eisner, N. L., Barragán, O., Aigrain, S., et al. 2020, MNRAS, 494, 750 [NASA ADS] [CrossRef] [Google Scholar]
  54. Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262 [Google Scholar]
  55. Fernandes, C. S., Van Grootel, V., Salmon, S. J. A. J., et al. 2019, ApJ, 879, 94 [NASA ADS] [CrossRef] [Google Scholar]
  56. Foreman-Mackey, D. 2016, J. Open Source Softw., 1, 24 [Google Scholar]
  57. Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64 [Google Scholar]
  58. Fukui, A., Narita, N., Tristram, P. J., et al. 2011, PASJ, 63, 287 [NASA ADS] [Google Scholar]
  59. Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504 [Google Scholar]
  60. Gaia Collaboration (Brown, A.G.A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Gaia Collaboration (Brown, A.G.A., et al.) 2021, A&A, 650, C3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Gan, T., Soubkiou, A., Wang, S. X., et al. 2022, MNRAS, 514, 4120 [CrossRef] [Google Scholar]
  63. Garcia, L. J., Timmermans, M., Pozuelos, F. J., et al. 2022, MNRAS, 509, 4817 [Google Scholar]
  64. Gelman, A., & Rubin, D. B. 1992, Stat. Sci., 7, 457 [Google Scholar]
  65. Geman, S., & Geman, D. 1984, IEEE Trans. Pattern Anal. Mach. Intell., PAMI-6, 721 [CrossRef] [Google Scholar]
  66. Giacalone, S., & Dressing, C. D. 2020, Astrophysics Source Code Library [record ascl:2002.004] [Google Scholar]
  67. Giacalone, S., Dressing, C. D., Jensen, E. L. N., et al. 2021, AJ, 161, 24 [Google Scholar]
  68. Gillon, M. 2018, Nat. Astron., 2, 344 [Google Scholar]
  69. Gillon, M., Jehin, E., Magain, P., et al. 2011, Euro. Phys. J. Web Conf., 11, 06002 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Gillon, M., Triaud, A. H. M. J., Fortney, J. J., et al. 2012, A&A, 542, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Gillon, M., Jehin, E., Fumel, A., Magain, P., & Queloz, D. 2013, Euro. Phys. J. Web Conf., 47, 03001 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Gillon, M., Demory, B. O., Madhusudhan, N., et al. 2014, A&A, 563, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature, 533, 221 [Google Scholar]
  74. Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Nature, 542, 456 [NASA ADS] [CrossRef] [Google Scholar]
  75. Gillon, M., Meadows, V., Agol, E., et al. 2020, ArXiv e-prints [arXiv:2002.04798] [Google Scholar]
  76. Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39 [NASA ADS] [CrossRef] [Google Scholar]
  78. Günther, M. N., Pozuelos, F. J., Dittmann, J. A., et al. 2019, Nat. Astron., 3, 1099 [Google Scholar]
  79. Hadden, S. 2019, https://doi.org/10.5281/zenodo.3356829 [Google Scholar]
  80. Haghighipour, N., & Kaltenegger, L. 2013, ApJ, 777, 166 [Google Scholar]
  81. Hamuy, M., Walker, A. R., Suntzeff, N. B., et al. 1992, PASP, 104, 533 [NASA ADS] [CrossRef] [Google Scholar]
  82. Hamuy, M., Suntzeff, N. B., Heathcote, S. R., et al. 1994, PASP, 106, 566 [NASA ADS] [CrossRef] [Google Scholar]
  83. Hastings, W. K. 1970, Biometrika, 57, 97 [Google Scholar]
  84. Hinse, T. C., Haghighipour, N., Kostov, V. B., & Goździewski, K. 2015, ApJ, 799, 88 [NASA ADS] [CrossRef] [Google Scholar]
  85. Hippke, M., & Heller, R. 2019, A&A, 623, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Hippke, M., David, T. J., Mulders, G. D., & Heller, R. 2019, AJ, 158, 143 [Google Scholar]
  87. Hirano, T., Kuzuhara, M., Kotani, T., et al. 2020, PASJ, 72, 93 [Google Scholar]
  88. Howell, S. B., Everett, M. E., Sherry, W., Horch, E., & Ciardi, D. R. 2011, AJ, 142, 19 [Google Scholar]
  89. Howell, S. B., Everett, M. E., Horch, E. P., et al. 2016, ApJ, 829, L2 [NASA ADS] [CrossRef] [Google Scholar]
  90. Husser, T.-O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  91. Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
  92. Irwin, M. J., Lewis, J., Hodgkin, S., et al. 2004, SPIE Conf. Ser., 5493, 411 [Google Scholar]
  93. Jehin, E., Gillon, M., Queloz, D., et al. 2011, The Messenger, 145, 2 [NASA ADS] [Google Scholar]
  94. Jehin, E., Gillon, M., Queloz, D., et al. 2018, The Messenger, 174, 2 [NASA ADS] [Google Scholar]
  95. Jenkins, J. M. 2002, ApJ, 575, 493 [Google Scholar]
  96. Jenkins, J. M., Caldwell, D. A., Chandrasekaran, H., et al. 2010, ApJ, 713, L87 [Google Scholar]
  97. Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, SPIE Conf. Ser., 9913, 99133E [Google Scholar]
  98. Jenkins, J. S., Pozuelos, F. J., Tuomi, M., et al. 2019, MNRAS, 490, 5585 [Google Scholar]
  99. Jenkins, J. M., Tenenbaum, P., Seader, S., et al. 2020, Kepler Data Processing Handbook: Transiting Planet Search, Kepler Science Document KSCI-19081-003 [Google Scholar]
  100. Kaltenegger, L. 2017, ARA&A, 55, 433 [NASA ADS] [CrossRef] [Google Scholar]
  101. Kaltenegger, L., & Haghighipour, N. 2013, ApJ, 777, 165 [NASA ADS] [CrossRef] [Google Scholar]
  102. Kaltenegger, L., & Traub, W. A. 2009, ApJ, 698, 519 [Google Scholar]
  103. Kane, S. R., & Hinkel, N. R. 2013, ApJ, 762, 7 [Google Scholar]
  104. Kasting, J. F., Whitmire, D. P., & Reynolds, R.T. 1993, Icarus, 101, 108 [NASA ADS] [CrossRef] [Google Scholar]
  105. Kipping, D. M. 2013, MNRAS, 435, 2152 [Google Scholar]
  106. Kirkpatrick, J. D., Looper, D. L., Burgasser, A. J., et al. 2010, ApJS, 190, 100 [Google Scholar]
  107. Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
  108. Kopparapu, R. K., Ramirez, R. M., SchottelKotte, J., et al. 2014, ApJ, 787, L29 [Google Scholar]
  109. Kopparapu, Ravi Kumar, R., Wolf, E.T., Haqq-Misra, J., et al. 2016, ApJ, 819, 84 [NASA ADS] [CrossRef] [Google Scholar]
  110. Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32 [Google Scholar]
  111. Kotani, T., Tamura, M., Nishikawa, J., et al. 2018, SPIE Conf. Ser., 10702, 1070211 [Google Scholar]
  112. Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369 [Google Scholar]
  113. Leconte, J., Chabrier, G., Baraffe, I., & Levrard, B. 2010, A&A, 516, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  114. Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013, Nature, 504, 268 [Google Scholar]
  115. Lépine, S., Rich, R. M., & Shara, M. M. 2003, AJ, 125, 1598 [CrossRef] [Google Scholar]
  116. Lépine, S., Rich, R. M., & Shara, M. M. 2007, ApJ, 669, 1235 [CrossRef] [Google Scholar]
  117. Lépine, S., Hilton, E. J., Mann, A. W., et al. 2013, AJ, 145, 102 [Google Scholar]
  118. Lester, K. V., Matson, R. A., Howell, S. B., et al. 2021, AJ, 162, 75 [NASA ADS] [CrossRef] [Google Scholar]
  119. Li, J., Tenenbaum, P., Twicken, J. D., et al. 2019, PASP, 131, 024506 [NASA ADS] [CrossRef] [Google Scholar]
  120. Lienhard, F., Queloz, D., Gillon, M., et al. 2020, MNRAS, 497, 3790 [NASA ADS] [CrossRef] [Google Scholar]
  121. Lin, Z., MacDonald, R. J., Kaltenegger, L., & Wilson, D. J. 2021, MNRAS, 505, 3562 [NASA ADS] [CrossRef] [Google Scholar]
  122. Lissauer, J. J., Marcy, G. W., Rowe, J. F., et al. 2012, ApJ, 750, 112 [NASA ADS] [CrossRef] [Google Scholar]
  123. Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122 [Google Scholar]
  124. Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
  125. Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nat. Astron., 1, 0129 [Google Scholar]
  126. Lustig-Yaeger, J., Meadows, V. S., & Lincowski, A. P. 2019, AJ, 158, 27 [NASA ADS] [CrossRef] [Google Scholar]
  127. Luyten, W. J. 1979, New Luyten catalogue of stars with proper motions larger than two tenths of an arcsecond; and first supplement; NLTT. (Minneapolis (1979)); Label 12 = short description; Label 13 = documentation by Warren; Label 14 = catalogue [Google Scholar]
  128. Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [Google Scholar]
  129. Mann, A. W., Brewer, J. M., Gaidos, E., Lépine, S., & Hilton, E. J. 2013, AJ, 145, 52 [Google Scholar]
  130. Mann, A. W., Deacon, N. R., Gaidos, E., et al. 2014, AJ, 147, 160 [CrossRef] [Google Scholar]
  131. Marconi, A., Abreu, M., Adibekyan, V., et al. 2021, The Messenger, 182, 27 [NASA ADS] [Google Scholar]
  132. McCormac, J., Pollacco, D., Skillen, I., et al. 2013, PASP, 125, 548 [Google Scholar]
  133. McCully, C., Volgenau, N. H., Harbeck, D.-R., et al. 2018, SPIE Conf. Ser., 10707, 107070K [Google Scholar]
  134. Ment, K., Dittmann, J. A., Astudillo-Defru, N., et al. 2019, AJ, 157, 32 [NASA ADS] [CrossRef] [Google Scholar]
  135. Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, J. Chem. Phys., 21, 1087 [Google Scholar]
  136. Mignard, F. 1979, Moon Planets, 20, 301 [Google Scholar]
  137. Miguel, Y., Cridland, A., Ormel, C. W., Fortney, J. J., & Ida, S. 2020, MNRAS, 491, 1998 [NASA ADS] [Google Scholar]
  138. Miller-Ricci, E., Seager, S., & Sasselov, D. 2009, ApJ, 690, 1056 [Google Scholar]
  139. Minkowski, R. L., & Abell, G. O. 1963, The National Geographic Society- Palomar Observatory Sky Survey, ed. K.A. Strand (Chicago: University of Chicago Pres), 481 [Google Scholar]
  140. Morello, G., Parviainen, H., Murgas, F., et al. 2022, A&A, submitted [arXiv:2201.13274] [Google Scholar]
  141. Mori, M., Livingston, J. H., Leon, J.D., et al. 2022, AJ, 163, 298 [NASA ADS] [CrossRef] [Google Scholar]
  142. Muirhead, P. S., Dressing, C. D., Mann, A. W., et al. 2018, AJ, 155, 180 [NASA ADS] [CrossRef] [Google Scholar]
  143. Mulders, G. D., Drtţżkowska, J., van der Marel, N., Ciesla, F.J., & Pascucci, I. 2021, ApJ, 920, L1 [NASA ADS] [CrossRef] [Google Scholar]
  144. Murray, C. A., Delrez, L., Pedersen, P. P., et al. 2020, MNRAS, 495, 2446 [Google Scholar]
  145. Murray, C. A., Queloz, D., Gillon, M., et al. 2022, MNRAS, 513, 2615 [NASA ADS] [CrossRef] [Google Scholar]
  146. Narita, N., Fukui, A., Yamamuro, T., et al. 2020, SPIE Conf. Ser., 11447, 114475K [NASA ADS] [Google Scholar]
  147. Neron de Surgy, O., & Laskar, J. 1997, A&A, 318, 975 [NASA ADS] [Google Scholar]
  148. Newton, E. R., Irwin, J., Charbonneau, D., et al. 2017, ApJ, 834, 85 [Google Scholar]
  149. Niraula, P., Wit, J.D., Rackham, B.V., et al. 2020, AJ, 160, 172 [NASA ADS] [CrossRef] [Google Scholar]
  150. O’Malley-James, J.T., & Kaltenegger, L. 2017, MNRAS, 469, L26 [CrossRef] [Google Scholar]
  151. Ormel, C. W., Liu, B., & Schoonenberg, D. 2017, A&A, 604, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  152. Parviainen, H. 2015, MNRAS, 450, 3233 [Google Scholar]
  153. Parviainen, H., & Aigrain, S. 2015, MNRAS, 453, 3821 [Google Scholar]
  154. Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  155. Pidhorodetska, D., Fauchez, T. J., Villanueva, G. L., Domagal-Goldman, S. D., & Kopparapu, R. K. 2020, ApJ, 898, L33 [NASA ADS] [CrossRef] [Google Scholar]
  156. Pierrehumbert, R., & Gaidos, E. 2011, The Astrophysical Journal, 734, L13 [NASA ADS] [CrossRef] [Google Scholar]
  157. 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]
  158. Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes in FORTRAN. The art of Scientific Computing (Cambridge: Cambridge University Press) [Google Scholar]
  159. Ramirez, R. M., & Kaltenegger, L. 2016, ApJ, 823, 6 [NASA ADS] [CrossRef] [Google Scholar]
  160. Ramirez, R. M., & Kaltenegger, L. 2017, ApJ, 837, L4 [Google Scholar]
  161. Ramirez, R. M., & Kaltenegger, L. 2018, ApJ, 858, 72 [Google Scholar]
  162. Rasmussen, C. E., & Williams, C. K. I. 2006, Gaussian Processes for Machine Learning (Cambridge, USA: The MIT Press) [Google Scholar]
  163. Raymond, S. N., Izidoro, A., Bolmont, E., et al. 2022, Nat. Astron., 6, 80 [NASA ADS] [CrossRef] [Google Scholar]
  164. Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, PASP, 115, 362 [NASA ADS] [CrossRef] [Google Scholar]
  165. Reid, I. N., Brewer, C., Brucato, R. J., et al. 1991, PASP, 103, 661 [NASA ADS] [CrossRef] [Google Scholar]
  166. Rein, H., & Liu, S.-F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  167. Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [Google Scholar]
  168. Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Teles. Instrum. Syst., 1, 014003 [Google Scholar]
  169. Riddick, F. C., Roche, P. F., & Lucas, P. W. 2007, MNRAS, 381, 1067 [Google Scholar]
  170. Rojas-Ayala, B., Covey, K. R., Muirhead, P.S., & Lloyd, J.P. 2012, ApJ, 748, 93 [NASA ADS] [CrossRef] [Google Scholar]
  171. Sabin, L., Gómez Maqueo Chew, Y., Demory, B.-O., Petrucci, R., & Saint-Ex Consortium 2018, in 20th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun, 59 [Google Scholar]
  172. Savitzky, A., & Golay, M. J. E. 1964, Anal. Chem., 36, 1627 [Google Scholar]
  173. Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
  174. Schanche, N., Pozuelos, F. J., Günther, M. N., et al. 2022, A&A, 657, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  175. Schwarz, G. 1978, Ann. Stat., 6, 461 [Google Scholar]
  176. Scott, N. J., Howell, S. B., Gnilka, C. L., et al. 2021, Front. Astron. Space Sci., 8, 138 [NASA ADS] [CrossRef] [Google Scholar]
  177. Seager, S. & Mallén-Ornelas, G. 2003, ApJ, 585, 1038 [NASA ADS] [CrossRef] [Google Scholar]
  178. Sebastian, D., Gillon, M., Ducrot, E., et al. 2021, A&A, 645, A100 [EDP Sciences] [Google Scholar]
  179. Seifahrt, A., Bean, J. L., Stürmer, J., et al. 2020, SPIE, 11447, 305 [Google Scholar]
  180. Sergeev, D. E., Fauchez, T. J., Turbet, M., et al. 2021, Planet. Sci. J, 3, 212 [Google Scholar]
  181. Sharma, S., Stello, D., Buder, S., et al. 2018, MNRAS, 473, 2004 [NASA ADS] [CrossRef] [Google Scholar]
  182. Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [Google Scholar]
  183. Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000 [Google Scholar]
  184. Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
  185. Stassun, K. G., & Torres, G. 2016, AJ, 152, 180 [Google Scholar]
  186. Stassun, K. G., & Torres, G. 2021, ApJ, 907, L33 [NASA ADS] [CrossRef] [Google Scholar]
  187. Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136 [Google Scholar]
  188. Stassun, K. G., Corsaro, E., Pepper, J. A., & Gaudi, B. S. 2018a, AJ, 155, 22 [Google Scholar]
  189. Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018b, AJ, 156, 102 [Google Scholar]
  190. Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985 [Google Scholar]
  191. Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100 [Google Scholar]
  192. Sullivan, P. W., Winn, J. N., Berta-Thompson, Z. K., et al. 2015, ApJ, 809, 77 [Google Scholar]
  193. Tamura, M., Suto, H., Nishikawa, J., et al. 2012, SPIE Conf. Ser., 8446, 84461T [NASA ADS] [Google Scholar]
  194. Teyssandier, J., Libert, A. S., & Agol, E. 2022, A&A, 658, A170 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  195. Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  196. Turbet, M., Fauchez, T. J., Sergeev, D. E., et al. 2021, Planet. Sci. J, 3, 211 [Google Scholar]
  197. Vanderspek, R., Huang, C. X., Vanderburg, A., et al. 2019, ApJ, 871, L24 [Google Scholar]
  198. Wells, R. D., Rackham, B. V., Schanche, N., et al. 2021, A&A, 653, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  199. Wildi, F., Blind, N., Reshetov, V., et al. 2017, SPIE Conf. Ser., 10400, 1040018 [Google Scholar]
  200. Wolf, E. T., & Toon, O. B. 2014, Geophys. Res. Lett., 41, 167 [Google Scholar]
  201. Yang, J., Cowan, N. B., & Abbot, D. S. 2013, ApJ, 771, L45 [Google Scholar]
  202. Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49 [NASA ADS] [EDP Sciences] [Google Scholar]
  203. Zhang, S., Luo, A. L., Comte, G., et al. 2019, ApJS, 240, 31 [NASA ADS] [CrossRef] [Google Scholar]
  204. Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502 [Google Scholar]
  205. Van Grootel, V., Fernandes, C. S., Gillon, M., et al. 2018, ApJ, 853, 30 [Google Scholar]
1

Internally, the SPECULOOS team also refers to this discovery as SPECULOOS-1, since it was the first official detection as part of the SPECULOOS project.

6

The luminosity of very low-mass stars evolves very slowly with time once the star has turned on core H-burning and has reached the main sequence. Hence, assuming any age ≳2 Gyr will provide the same resulting mass.

7

The correlation statistic is defined as where is the whitened core or halo aperture flux time series and s is the whitened transit model light curve for the given planet candidate (Twicken et al. 2018).

9

The SHERLOCK (Searching for Hints of Exoplanets fRom Lightcurves Of spaCe-based seeKers) code is fully available on GitHub: https://github.com/franpoz/SHERLOCK

10

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

11

The lower the order of the mean motion resonance, the larger the induced TTV amplitude.

All Tables

Table 1

Ground-based transit light curves of LP 890-9 used in our global transit analysis.

In the text
Table 2

Properties of the host star.

In the text
Table 3

Properties of the LP 890-9 planetary system based on our global transit analysis (see Sect. 5.1.1).

In the text
Table 4

Transit depths returned by our global data analysis for which we allowed transit depth variations between the different bandpasses (see Sect. 5.1.2).

In the text
Table 5

Transit timings returned by our global data analysis for which we allowed TTVs from the linear transit ephemerides defined by the orbital period P and mid-transit time T0 given in Table 3 (see Sect. 5.1.3).

In the text
Table A.1

Subaru/IRD radial velocity measurements. We recommend treating the YJ-band and H-band RVs as independent datasets, as done in the analysis presented in Sect. 6.

In the text
Table B.1

Photometric baseline models, residual RMS, and error scaling factors βw and βr for each transit light curve used in our global transit analysis (see Sect. 5.1.1). For the baseline function, p(αn) denotes a n-order polynomial function of the parameter α, with α that can be t = time, a = airmass, x and y = x- and y- position of the target on the detector, or f = full width at half maximum of the point spread function.

In the text

All Figures

thumbnail Fig. 1

TESS target pixel files of LP 890-9 for sectors 4 (upper left), 5 (upper right), 31 (lower left), and 32 (lower right). The pixel scale is 21″ per pixel. The electron counts are colour-coded. The target is indicated by a white cross. Nearby sources identified in Gaia DR2 (Gaia Collaboration 2018), up to 6 magnitudes in contrast with LP 890-9, are shown as red circles. The symbol size is proportional to the magnitude contrast with the target. The apertures used by the SPOC pipeline to extract the photometry are shown as red shaded regions. These plots were created with tpfplotter (Aller et al. 2020).
In the text
thumbnail Fig. 2

TESS photometry of LP 890-9. For each of the four sectors, the 2-min data points (in grey) have been binned into 30-min intervals to produce the black points, with error bars corresponding to the root-mean-square of the uncertainties of the points in the bins. The transits of LP 890-9 b and c are indicated by red and blue dotted lines, respectively. The region marked in orange in sector 4 was impacted by thermal effects and thus excluded from our analysis.
In the text
thumbnail Fig. 3

Zorro speckle imaging 5σ contrast curves at 562 (blue) and 832 nm (red), along with the 832 nm reconstructed image.
In the text
thumbnail Fig. 4

Lick/Kast red optical spectrum of LP 890-9 (solid black lines) compared to M5 (dashed magenta line), M6 (solid magenta line), and M7 (dotted magenta line) optical spectral templates from Bochanski et al. (2007). The three spectral templates are offset by a constant to facilitate comparison. Key atomic and molecular spectral features are labelled, and the inset box shows a close-up of the region around the 6563 Å Ha emission line and 6708 Å Li I absorption line (not detected).
In the text
thumbnail Fig. 5

SpeX spectrum of LP 890-9 (blue) obtained in August 2021. The SpeX Prism spectrum of the M6.0 standard LHS 1375 (Kirkpatrick et al. 2010) is shown in grey for comparison. Prominent spectral features are highlighted, and regions of strong telluric absorption are shaded. The bottom line gives the measurement uncertainty.
In the text
thumbnail Fig. 6

Spectral energy distribution (SED) of LP 890-9. The red symbols represent the observed photometric measurements (the horizontal bars represent the effective width of the bandpass). The blue symbols are the model fluxes from the best-fit NextGen atmosphere model (black curve).
In the text
thumbnail Fig. 7

Distribution of ages in the full GALAH DR3 sample (light grey histogram) and sources that match the UVW and metallicity of LP 890-9 to within 10 km s−1 and 1σ, respectively (dark grey histogram). The latter has a broad peak at Gyr (yellow line with uncertainties indicated by the green shaded region), which we take as an estimate for the age of LP 890-9.
In the text
thumbnail Fig. 8

Imaging of the present-day position of LP 890-9 over 64 yr, to assess for a current blend with an unresolved background object. Archival images from POSS I (1957), POSS II (1994), and Pan-STARRS (2014) are presented, along with a stack image from one night of observation with SSO/Europa (2021). For each image, the position of LP 890-9 during the SSO/Europa observations is shown as a white circle.
In the text
thumbnail Fig. 9

Phase-folded detrended transit photometry of LP 890-9 b in each observed bandpass. The unbinned data points are shown in grey, while the black circles with error bars correspond to 10-min bins. The best-fit transit model is shown in red, with a different limb darkening for each bandpass.
In the text
thumbnail Fig. 10

Phase-folded detrended transit photometry of LP 890-9 c in each observed bandpass. The unbinned data points are shown in grey, while the black circles with error bars correspond to 10-min bins. The best-fit transit model is shown in red, with a different limb darkening for each bandpass. We note that the MuSCAT3 r′- and i′- light curves show a possible hump around mid-transit which could be related to spot crossing.
In the text
thumbnail Fig. 11

TTVs of the two planets measured from the ground-based transit photometry (see Sect. 5.1.3).
In the text
thumbnail Fig. 12

Results from the transit injection-and-recovery tests performed on the TESS (upper panel) and SSO (lower panel) data to assess the detectability of transiting planets in the LP 890-9 system. Upper panel: we explored a total of 3600 different scenarios. Each pixel evaluated 32 scenarios, that is, 32 light curves with injected planets having different P, Rp, and T0. Larger recovery rates are presented in yellow and green, while lower recovery rates are shown in blue and darker hues. LP 890-9 b (recovery rate ~60%) and c (recovery rate ~6%) are displayed as red and blue stars, respectively (see Sect. 5.2.1 for details). Lower panel: we injected 6000 artificial planets into the combined I + z-filter SSO light curve. On average, each box in this plot contains 50 injection scenarios with differing P, Rp, ip, and ϕ. Similarly to the upper panel, we include LP 890-9 b (recovery rate of 100%) and c (recovery rate of 96%) as red and blue stars, respectively (see Sect. 5.2.2 for details).
In the text
thumbnail Fig. 13

Phase coverage as a function of the orbital period using all photometric data obtained with SSO for LP 890-9. The effective coverage of 86% is the integral of the black curve. The orbital periods of planets b and c are shown as orange and red lines, respectively. The figure also shows the empirical (optimistic) HZ in green, i.e. the region between the ‘Recent Venus’ and ‘Early Mars’ limits (Kopparapu et al. 2013, 2014, and Sect. 8.3).
In the text
thumbnail Fig. 14

Global I + z-filter light curves for the SSO/Io (red) and SSO/Europa (blue) observations. The black points correspond to 20-min bins.
In the text
thumbnail Fig. 15

Lomb-Scargle periodograms for SSO/Io (top, red) and SSO/Europa (bottom, blue). The five highest peaks are marked by black crosses in both cases. For SSO/Europa, the most significant peak is around 50–55 days and the next highest four are at ~1 day and were determined to be aliases. The results for SSO/Io are less clear due to large gaps and much less data. The most significant peaks are found at 23.0, 54.0, and 4.3 days. The peak around 54 days could correspond to the one detected for SSO/Europa, but the other two are not present in SSO/Europa’s periodogram.
In the text
thumbnail Fig. 16

Subaru/IRD radial velocities obtained in the YJ- (grey squares) and H- (black circles) bands. The data are plotted with both their original error bars (thick solid lines) and the error bars enlarged by the best-fit jitter terms (thinner transparent lines). The top panel displays the data as a function of time. The solid red line shows the 2-planet model generated from the posterior median, while semi-transparent red lines show 25 models randomly drawn from the posteriors. The residuals around the median solution are shown in the second panel. The two lower panels show the data folded on the orbital periods of planets b and c, respectively, after removing the RV component from the other planet.
In the text
thumbnail Fig. 17

Exoplanets found around cool dwarfs with an effective temperature lower than 3600 K using the transit (black circles) or radial velocity (grey triangles) techniques. Data were taken from NASA Exoplanet Archive on 20 April 2022. The planets are shown as a function of their incident stellar flux and the effective temperature of their host star. The size of the symbols is proportional to the planetary radius, estimated from the mass using the probabilistic mass-radius relationship of Chen & Kipping (2017) when it is unknown (for non-transiting RV planets). LP 890-9 b and c are shown as orange circles. The inner (Recent Venus) and outer (Early Mars) boundaries of the empirical (optimistic) HZ are shown as solid red and blue lines, respectively (Kopparapu et al. 2014 and Sect. 8.3). The dashed green line shows an alternative conservative inner edge limit from 3D models (Leconte et al. 2013). The outer edge of the HZ in 3D models agrees with incident stellar flux for Early Mars insulation. We note that the three limits shown do vary with the stellar effective temperature, but this is not obvious on this plot due to the log-scale of the x-axis.
In the text
thumbnail Fig. 18

Expected transmission signals (see Sect. 8.4 for details) of currently known transiting terrestrial exoplanets (Rp < 1.6 R⊕) as a function of their incident stellar flux. Data were taken from NASA Exoplanet Archive on 29 March 2022. We only show planets with mild incident irradiations, lower than five times that received by the Earth (Sp < 5 S⊕). The symbols are colour-coded as a function of their S/N relative to TRAPPIST-1 b, obtained by scaling the signal amplitude with the hosts’ brightness in the J-band and using TRAPPIST-1 b’s S/N as a reference. The size of the symbols is proportional to the planetary radius. Transmission signals less than 5 ppm have been removed to enhance the readability of the figure.
In the text
thumbnail Fig. C.1

MCMC posterior distributions for the main fitted transit parameters (see Sect. 5.1.1): the log of the orbital period (log P), the mid-transit time (T0), the transit depth (dF), and the cosine of the orbital inclination (cos ip) for each of the two planets; and the log of the stellar density (log ρ★). This plot was made using the corner python package (Foreman-Mackey 2016).
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.