© The Authors 2023

1 Introduction

The Transiting Exoplanet Survey Satellite (TESS) NASA mission (Ricker et al. 2015) is a powerful observatory for detecting potential transiting candidates. To date, it has provided more than 4000 candidates. The Echelle Spectrograph for Rocky Exoplanet and Stable Spectroscopic Observations, ESPRESSO (Pepe et al. 2021), operating at the Very Large Telecope (VLT) in Chile, serves the community as an efficient tool for following these candidates up and allows determining radial velocity (RV) measurements with a 10 cm s⁻¹ precision. Because radial velocity observations are limited by photon noise precision and instrumental instability, it is natural to use ESPRESSO to determine the masses of small-size planets with a radius smaller than two Earth radii (R_⊕). The instrument is indeed mounted on the VLT, and its initial design was oriented towards high stability. This is necessary requirement for reaching the precision that is required to detect these planets, which is one of the objective of the ESPRESSO guaranteed time observation (GTO) consortium (Program ID 1102.C-744, PI: F. Pepe). The subprogram of the ESPRESSO GTO that focuses on determining precise masses of TESS and K2 transiting candidates has produced several results since its start on October 1, 2018: Damasso et al. (2020) derived a precise characterisation of the Π Mensae system, which is composed of a brown dwarf and a 4.3 ± 0.7 Earth-mass (M_⊕) super-Earth; the masses of the super-Earth and sub-Neptune planets around the G2 star K2-38 were determined by Toledo-Padrón et al. (2020) as M_⊕ and 8.3 ± 1.3 M_⊕, respectively; Sozzetti et al. (2021) unveiled two planets around the bright late-F dwarf HD 5278, a M_⊕ transiting sub-Neptune and a non-transiting Neptune at periods of 14 and 41 days, respectively; Mortier et al. (2020) reported two planets in near-resonance around the slightly evolved G2 star K2-111 at 5.3 and 15.7 days with a mass of M_⊕ for the inner planet that transits in front of its host star; and finally, Demangeon et al. (2021) announced the detection of the lowest-mass planet measured so far around L 98-59, a nearby stellar system with four planets.

The determination of the density for planets with measured radius and mass allows inferring their bulk composition. Current observations suggest that exoplanets follow two distinct trends: similar to an Earth-like composition for planets with a mass of up to ~25 M_⊕ and a more volatile-rich composition for planets more massive than ~10 M_⊕ (Otegi et al. 2020). The transition of these two populations and their origins are key questions in exoplanet research and require a precise mass determination for both populations. Using an updated exoplanet catalogue based on reliable and robust masses, Otegi et al. (2020) showed that the dispersion around the Earth-like composition for the rocky exoplanet population is small. Theoretical models show inherent degeneracies in the composition of planets with similar masses and radii, especially in the range of Earth masses to Neptune masses (e.g. Sotin et al. 2007; Seager et al. 2007; Howe et al. 2014; Dorn et al. 2017); however, the lack of precise measurements of the two physical variables remains the dominant source of uncertainty in this analysis. Efforts to increase the number of planets with accurate mass and radius measurements are therefore paramount. In the near future, the need for a broad sample of planets with a high-accuracy density will be necessary to identify suitable targets for an atmospheric characterisation with the next generation of observatories, such as the James Webb Space Telescope (JWST; Gardner et al. 2006; Deming et al. 2009) or the Extremely Large Telescope (ELT; Gilmozzi & Spyromilio 2007; Marconi et al. 2016, 2020, 2021), and also to serve as input in the atmospheric models that will study these targets.

M dwarfs account for roughly 75% of all the stars in our Galaxy (Henry et al. 2006). They provide stronger radial velocity signals and a high contrast ratio and transit depth for transiting Earth-like planets that favours the detection of spectral features in the atmosphere of these planets (Wunderlich et al. 2019). The probability of hosting multiple planets in compact orbits is higher for M dwarfs than other stars (Muirhead et al. 2015; Ballard 2019). Hence planets transiting M dwarfs are prime targets in our quest for accurate mass and radius measurements, based on which planets with an atmosphere can be identified. We also wish to be able to study planet composition as a function of stellar properties, and it is therefore required to first detect and characterise planets around stars with lower masses.

Stellar multiplicity is an increasing function of the primary mass (Duchêne & Kraus 2013). The multiplicity rate, that is, the percentage of stellar systems with more than one star, is therefore higher for solar-type stars (44%; Duquennoy & Mayor 1991; Raghavan et al. 2010; Tokovinin 2014) than for M dwarfs (27%; Leinert et al. 1997; Duchêne & Kraus 2013; Ward-Duong et al. 2015; Winters et al. 2019a). For these two stellar types, binaries form the bulk of the contingent and triple star systems account for 8% and 3% of the overall population of solar-type stars and M dwarfs, respectively. The study of planets within these systems provides a unique angle to help shape our global understanding of the planet formation mechanisms. The required stability over a period of time of planetary systems in binary or higher-order systems that is long enough adds complexity to the problem of habitability around M dwarfs (Shields et al. 2016) and forces a deeper determination of the very notion of habitability.

The LTT 1445 ABC system is a triple stellar system composed of three M dwarfs in a hierarchical configuration. A compact binary orbits the main component at a minimum separation of ~21 astronomical units (AU). A TESS candidate has been confirmed by Winters et al. (2019b) with only one TESS sector. The system is at a distance of only 6.9 parsecs (pc) from Earth, which makes it the second closest transiting exoplanet system to date after HD 219134 (Motalebi et al. 2015). The planet orbits the main component LTT 1445A with an orbital period of ~5.4 days and has a radius of ~1.4 R_⊕, which makes it an ideal target for ESPRESSO follow-up.

In this paper, we present radial velocity measurements obtained with ESPRESSO of the M dwarf LTT 1445A (TOI-455) that allow us to determine the mass for the previously announced transiting planet b with an orbital period of 5.4 days with high accuracy. We also report a new transiting planet c at 3.12 days and a potential non-transiting planet d at 24.30 days. Section 2 presents the observations and Sect. 3 the stellar properties. In Sect. 4, we detail the photometric analysis of the TESS observations, including an additional sector with respect to Winters et al. (2019b). In Sect. 5, we present the radial velocity analysis with an investigation of the stellar activity of the host star. The last section reports the results from the combined photometric and spectroscopic analysis and a final discussion.

2 Observations

2.1 TESS photometry

The LTT 1445 system has been observed by TESS in sector 4 (from 19 October to 14 November 2018; a 26-day interval) and sector 31 (from 22 October to 18 November 2020; a 27-day interval). Both sectors were observed with CCD 4 of Camera 2 with a 2-min cadence. The 20-s cadence mode is also available for sector 31. Following the method of Lillo-Box et al. (2020), we used tpfplotter¹ (Aller et al. 2020) to scan for any contaminant sources. As shown in Fig. 1, we did not find any contaminants with magnitude contrast brighter than 8 in the automatically selected aperture.

Winters et al. (2019b) detected a transiting planet with an orbital period of ~5.4 days using sector 4 around LTT 1445A. The publicly available TESS photometric data were downloaded from the Mikulski Archive for Space Telescopes (MAST).

2.2 Spectroscopic follow-up with ESPRESSO

We carried out the radial velocity monitoring of the system with the ESPRESSO spectrograph installed at the VLT in the Paranal Observatory in Chile. The observations were part of the subprogram of the GTO that focuses on determining precise masses of TESS and K2 transiting candidates. All observations were carried out with the 1UT mode, no binning, and 900-s exposure time with the calibration source (fiber B) pointed to the Fabry-Perot. The spectroscopic data were extracted using version 2.3.0 of the ESPRESSO data reduction pipeline (Lovis et al., in prep.), which is available for direct download from the ESO pipelines website². The data reduction was performed with an M4 mask. The extracted data are downloaded from the DACE platform API³. The ESPRESSO dataset consists of 84 observations taken between 20 July 2019 and 12 March 2021 (the observation with ESO OB Id 2549082 was stopped after ~500 s of exposure time and is not considered in this dataset). The gap in the time coverage of these observations is due to the COVID crisis, which caused the instrument at the Paranal observatory to be unavailable for nearly nine months. The ESPRESSO dataset is combined with 14 observations that are available in the ESO HARPS archive (Program ID 072.C-0488, PI: M. Mayor; Program IDs 183.C-0437 and 1102.C-0339, PI: X. Bonfils) that were also used in Winters et al. (2019b). The median uncertainty for the radial velocity time-series is 0.48 m s⁻¹ and 1.16 m s⁻¹ for the ESPRESSO and HARPS dataseis, respectively. The combined dataset (Fig. A.5) shows a clear long-term drift due to the interaction of LTT 1445A with the ВС component of the triple system.

Fig. 1 TPF of LTT 1445 from the TESS observations in Sector 4 (composed with tpfplotter; Aller et al. 2020). The SPOC pipeline aperture is overplotted with shaded red squares, and the Gaia DR3 catalogue is also overlaid with symbol sizes proportional to the magnitude contrast with the target, marked with a white cross. We display all sources from the Gaia DR3 catalogue with a magnitude contrast up to Δm = 8 mag.

3 Stellar properties

3.1 Spectral energy distribution

We estimated the physical parameters of LTT 1445A by fitting the photometric spectral energy distribution constructed from catalogue broadband magnitudes following the prescription of the Virtual Observatory VOSA tool (Bayo et al. 2008; Rodrigo et al. 2020). The available photometry includes optical magnitudes from Gaia, APASS, Tycho-2, and the Sloan Digital Sky Survey catalogues (Henden et al. 2015; Høg et al. 2000; Blanton et al. 2017), the JHK magnitudes from the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), and the W1 − 4 magnitudes from the Wide-field Infrared Survey Explorer mission (A11WISE; Wright et al. 2010; Mainzer et al. 2011). We employed the BT-Settl models (Allard et al. 2012), which are computed for a wide range of effective temperatures (T_eff), surface gravities (log g), and atmospheric metallicities ([Fe/H]) with steps of 100 K, 0.5 dex, and 0.5 dex, respectively. The photometric spectral energy distribution of LTT 1445A covers the wavelength interval 0.36−22.09 μm. We assumed no interstellar extinction due to the proximity of the star (d = 6.864 ± 0.001 pc) and no obvious reddening in the stellar colours at the long wavelengths of AllWISE data. From this analysis, we obtained T_eff = 3300 ± 100 K, log g = 6.0 ± 0.50 dex, and [Fe/H] = −1 dex. The high surface gravity is consistent with the star being older than ~ 1 Gyr. The observed spectral energy distribution of LTT 1445A and its best fit are shown in Fig. 2. By integrating the observed spectral energy distribution, we derived the stellar luminosity to be log L/L_⊙ = −2.099 ± 0.010 dex.

Fig. 2 Observed photometric spectral energy distribution of LTT 1445A (red dots) from the и filter through the W4 band. The best BT-Settl model (Allard et al. 2012), corresponding to Teff = 3300 K, low metallicity, and log g = 6.0 (cgs), is shown with the blue line.

3.2 Mass and radius of LTT 1445A

The location of LTT 1445A in the Gaia colour-magnitude diagram is displayed in Fig. 3, where the position of the fainter binary LTT 1445BC is also indicated for completeness. The magnitudes and colours of LTT 1445A are consistent with M4.25 ± 0.25. All stellar sequences shown in the figure were built using Gaia photometry and parallaxes (see Luhman 2018 for the young isochrones and Cifuentes et al. 2020 for the main-sequence track), thus minimising any possible systematics in the comparison exercise. The absolute G-band magnitude of LTT 1445A is fainter than that of the main sequence for its colour, which is consistent with a slightly sub-solar metallicity, as was inferred from the analysis of the stellar photometric spectral energy distribution. The position of LTT 1445BC appears to lie closer to the main sequence of stars, but the two components of the subarcsecond binary were not resolved by Gaia. When the photometry was deblended into individual magnitudes following the delta-magnitudes provided by Winters et al. (2019b), the positions of LTT 1445B and C lie at fainter magnitudes in Fig. 3 and away from the solar-metallicity track. The colours and absolute magnitudes of the stars are compatible with their being main-sequence M-type stars.

At low stellar masses, results derived from stellar models (either evolutionary or atmosphere models) strongly depend on the assumed physics, which are not perfectly controlled at low temperatures. To avoid this problem as much as we can, we therefore decided to base our mass and radius deteminations of LTT 1445A on relations constructed from compilations of detached double-lined double-eclipsing main-sequence M-dwarf binaries (Mann et al. 2019). Using the most recent Gaia distance, the 2MASS K magnitude, and the empirical mass–M_K (Mann et al. 2019) and radius–M_K (Mann et al. 2015) relations, we determined the mass and size of LTT 1445A as M = 0.249 ± 0.023 M_⊙ and R = 0.276 ± 0.010 R_⊙ for solar metallicity. For a sub-solar metallicity of [Fe/H] = −1. 0 dex, the radius would slightly increase to R = 0.289 ± 0.010 R_⊙, while the stellar mass remains unchanged within the quoted uncertainty for metallicities around solar, as discussed in Mann et al. (2019). The errors in mass and radius were derived from the uncertainties in distance, K-band magnitude, and the dispersion of the empirical relations. These values are compatible with those of Winters et al. (2019b).

Fig. 3 Location of LTT 1445A and BC in the Gaia colour-magnitude diagram. The binary BC is not resolved by Gaia. The young stellar isochrones of the β Pic (blue) and Tucana-Horologium (green) moving groups and the Pleiades cluster (purple) are taken from Luhman (2018), and the main sequence of field stars (black) is taken from Cifuentes et al. (2020). The grey area represents the dispersion observed among main-sequence stars.

4 Photometry analysis with TESS

4.1 Identification of transit candidates

The three stars of the system LTT 1445 are included in the TESS aperture because of the large TESS pixel size of 21″. We used the lightkurve python package (Lightkurve Collaboration 2018) to download the three light curves produced by the Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014) algorithm for both sectors. The three light curves for the two sectors are relatively complex and exhibit a clear modulation signal at 1.4 days (Fig. 4). In this section, we perform a blind initial analysis to identify potential planet candidates.

4.1.1 Modulation signal at 1.4 days

The modulation signal with a period of 1.4 days is clearly visible in both sectors. Winters et al. (2019b) suspected that it originates from either the B or C component, based on the activity indicator measurements from spectra obtained with the Tillinghast Reflector Echelle Spectrograph. The amplitude of the signal varies between the two epochs with a root mean square (rms) of 2.37 and 4.13 part per thousand (ppt) during observation of sector 4 and sector 31, respectively.

Similarly to Winters et al. (2019b), we chose to model it as well as any other stellar signals with a Gaussian process (GP; Williams 2006). A GP is characterised by a covariance matrix generated by a covariance function (kernel) that determines the characteristics of the underlying model. The choice of kernel and its hyperparameters is arbitrary in the sense that it is not based on some physical first principles. However, it is conditioned by our understanding of the underlying physics (periodicities and timescales of the expected signals). Throughout this paper, we use a kernel with three hyperparameters: the amplitude B, the length L, and the period P. An arbitrary element, C_ij, of the covariance matrix is (1)

This kernel is a simpler version of the kernel described in Foreman-Mackey et al. (2017), which has an additional scale parameter C. Models using parameter C were systematically rejected when compared with models without it. We therefore did not justify the use of this additional parameter. This kernel was implemented with the celerite package (Foreman-Mackey et al. 2017).

The three PDCSAP light curves were normalised, time windows affected by flares (spotted by eye) were removed, and a 3-σ clipping was performed on the time series to remove outliers. The light curves were corrected for dilution, and we kept the dilution factor provided by the pipeline in the keyword CROWDSAP. The light curves were binned with 500 bins. All these operations were made with the lightkurve package. We then performed a Bayesian inference on the hyperparameters of the GP with the three binned light curves independently. The results are similar for the two cadences of sector 31. We also performed the inference with the two sectors combined using the 2-min cadence only. The inference was made with the nested sampling algorithm (Skilling 2006) implemented with pymultinest (Buchner et al. 2014), a Python wrapper for multinest (Feroz et al. 2019, 2009; Feroz & Hobson 2008). The nested sampling algorithm computes the Bayesian evidence (the logarithm of the Bayesian evidence is noted ln 𝒵), which allows comparing the models, but also provides the posterior distribution of the parameters as a byproduct. The three hyperparameters converged to similar values for the four inferences (see Fig. A.1).

4.1.2 Identification of two transit signals

The light curves were detrended using the mean prediction of the GP with the hyperparameters set to the maximum a posteriori (MAP) values obtained through the Bayesian inference on the combined dataset (2-min cadence for both sectors). We performed a periodogram analysis to identify potential transits within each light curve using the box least-squares (BLS) method for periods between 1 and 20 days and durations between 0.01 and 0.5 days. When a transit signal was found, it was masked to allow another iteration of the periodogram analysis. Results are shown in Appendix A.2. The strongest detected transit signal is linked to the planet that orbits LTT 1445A at 5.4 days (Winters et al. 2019b). A second signal is clearly detected in the two sectors at 3.1 days; it might originate from any of the three stars in the system. The light curves stripped of these two signals do not exhibit any significant signals. Fixing the hyperparameters to the MAP values obtained with the independent inference on each light curves does not affect the identification of these two signals.

The transit signals are indicated by b_i and c_i in Fig. 4. During the time window of the dataset for the 5.4-day signal, 9 transits took place. One transit (b₈) falls in a gap in the light curve. Transit b₃ is at the edge of a gap in sector 4 and was not included in Winters et al. (2019b). During the time window of the dataset for the 3.1-day signal, 16 transits occurred. One transit (c₃) falls in a gap in the sector 4 light curve. The transit depths are about 0.23% and 0.16% for the 5.4- and 3.1-day signals, respectively.

Fig. 4 LTT1445 ligt curves from TESS sectors 4 and 31. First row: TESS PDCSAP light curve for sector 4 with 2-min cadence (first column) and sector 31 (second column) for a 2-min cadence (black) and a 20-s cadence (grey). A clear modulation with a period of 1.4 days is clearly visible. Second row: Respective detrended light curves with the best-fit model (in red) from the analysis of Sect. 4.2.2. The vertical lines identify the transit for each planet: b1 to b10 for planet b, and c1 to c16 for planet c. The light blue dots in the 2-min cadence in the upper panels are the outliers in the dataset that were removed from the PDCSAP light curves.

4.2 TESS transit analysis

In this section, we perform a Bayesian inference with the 2-min cadence PDCSAP light curves of both sectors. The light curves were cleaned in a procedure similar to that described in Sect. 4.1.1, but without any binning (hereafter, the cleaned light curves). We assumed that the potential planet c orbits the same star as LTT 1445A b, which is confirmed by the radial velocities analysis.

4.2.1 Joint GP-transit model

Our model includes a limb-darkening transit model implemented with the Python package Batman (Kreidberg 2015) and a GP with the kernel from Sect. 4.1.1. The hyperparameters of the GP were fixed to the value obtained by the inference on the combined binned light curves. Each planet was modelled with at least five free parameters: the time of inferior conjunction T₀, the orbital period P, the planet radius in stellar radii R_p/R_*, the semi-major axis in stellar radii a/R_*, and the orbital inclination i for circular orbits. In case of eccentric orbits, planets were modelled with seven free parameters (the eccentricity e and the longitude of periastron ω are added to the five retrieved parameters described above). When the model included two or more planets, we enforced Kepler’s third law by reducing the parameter space of the additional planets by one parameter. Assuming that the masses of the planet are negligible compared to the mass of the star, we computed the semi-major axis from the period of the planet and the semi-major axis and orbital period of the first planet in the model. We used a quadratic limb-darkening law with two parameters u₁ and u₂. The parameter space was explored using the nested sampling algorithm. We combined at least three nested sampling runs with 1000 live points and a sampling efficiency of 0.8 to ensure a thorough exploration of the parameter space. The priors (similar to those in Table 1) were chosen uniform and broad, except for the periods and the time of inferior conjunction, which were uniform and centred on the values obtained in the quick analysis of Sect. 4.1.2 with a width broad enough to capture any shift of their values compared to the quick analysis.

4.2.2 Transit model selection

The Bayesian inference was first performed simultaneously with the 2-min cadence light curves. The Bayesian evidence computed by the nested sampling algorithm allows a comparison of different models with different numbers of parameters (see Trotta 2008 for insights on Bayesian inference and model selection using the Bayes factor). We tested four models to evaluate the robustness of the two transiting planets: a simple line (with one parameter: a constant), a model with only planet b (7 parameters), a model with planets b and c on circular orbits (11 parameters), and a model with planets b and c on eccentric orbits (15 parameters). The two-planet models are highly favoured, with a logarithm of the Bayes factor (ratio of the Bayesian evidence) higher than 100 between them and the two other models. The logarithm of the Bayes factor between the eccentric and the circular orbit models is only 1.7, which indicates a weak preference for the eccentric orbits because photometric observations only provide weak constraints on the orbital eccentricity. The two additional parameters (e and ω) improve the fit to the data, but not sufficiently so to compensate for the increase in parameter space volume, which results in a slight improvement of the Bayesian evidence. In a second step, we performed a Bayesian inference in which planets b and c were separately modelled with circular orbits on each light curve to probe any differences between the retrieved values of the planet parameters for the two sectors. We tested both cadences for sector 31, which yielded similar results within one sigma. Hence, we only used the results for the 2-min cadence of both sectors for the rest of our analysis.

4.2.3 Monitoring the joint GP-transit model: A linear model

In order to monitor the performance of the joint GP-transit model, we also used a more simple approach with a simple linear model and a reduced dataset covering the light curves around each transit only. Using the ephemeris obtained through the Bayesian inference with the joint GP-transit model, we identified each transit of each planet in the cleaned light curves. The data were cut at 2 h before and after each mid-transit time. We then hid the transit and performed a simple fit with a linear model on the narrow dataset. The data were then detrended with the best-fit linear model. For each sector, we performed a Bayesian inference with a circular transit model for the two planets on the detrended data. Figures A.3 and A.4 show the reduced dataset with the best fit from the Bayesian inference with the linear model and the joint GP-transit model (results from the previous section, the inference was made on the full light curves). Both models provide similar fits on the reduced dataset, showing that the GP component is almost linear around the transits, except for transit b3, which occurs in a curved part of the 1.4-day modulation signal. A cut in the data shorter than 2 h or a two-degree polynomial model did not significantly change the inference for the two planets.

4.2.4 Sensitivity of T₀ and TTV

The independent Bayesian inference on the time of inferior conjunction for planet b is significantly different between the two sectors, with values of 0 ± 0.008 and 0.81 ± 0.49 h for the light curves of sector 4 and sector 31 with a 2-min cadence, respectively (the first sector is used as reference). The value and uncertainty derived for sector 4 are similar to those derived by Winters et al. (2019b), who used a different GP kernel. The posterior distribution for the light curve of sector 31 is very broad compared to sector 4 (see Fig. 5). We explain this result by the higher amplitude of the 1.4-day modulation signal (Sect. 4.1.1) for this sector, which affects the ability of the joint GP-transit model to correctly catch the real value of the time of inferior conjunction of planet b. When the inference was made on the combined light-curve dataset, the retrieval of was driven by sector 4 light curve, which helps to retrieve the correct joint GP-transit model for the second sector. The uncertainty on is only marginally better than the retrieval with sector 4 alone, which indicates that most of the information on the time of inferior conjunction is extracted from the first sector light curve. The posterior distributions of the orbital period of planet b that were retrieved for the two sectors have a similar dispersion, with a small shift of the median value for sector 31 at less than one sigma of the value obtained for sector 4. The orbital period of planet b retrieved with the combined light-curve dataset is well constrained around the median value obtained for sector 4, but with a significant improvement because information is extracted from both sectors. A similar effect is visible for planet c. However, the individual transit signals are weaker, which results in a higher uncertainty on the retrieved value of the time of inferior conjunction that dominates any impact of the amplitude difference in the 1.4-day modulation between the two sectors. We obtain the same results when the simple linear detrending approach was used.

Fig. 5 Posterior distributions for the Bayesian inference with the joint GP-transit model (plain lines) and linear model (dashed lines) for the sector 4 light curve (black), for sector 31 with the 2-min cadence (orange), and on the combined 2-min cadence light-curve dataset (blue). From top to bottom: time of inferior conjunction, period, ratio of the planet to star radius and impact parameter for planet b (first column) and planet c (second column). Inference on sector 31 shows a broad distribution and a 1-h time transit variation and shift in the orbital period. This broad distribution seems to result from the higher amplitude of the polluting signal at 1.4 days for sector 31. The inference on the combined light curve dataset converges towards the same values as the inference on sector 4.

4.2.5 Grazing transit for planet c

The impact parameter of planet c (b_c = a_c/R_*× cos i_c) obtained with both approaches (the joint GP-transit model and the linear detrending model) is close to one. Figure 6 shows the posterior distribution of the inclination, the radius, and the semi-major axis retrieved for planet c. We calculated the grazing criterion R_p/R_* + b, which indicates a grazing transit if it is greater than one (Smalley et al. 2011). The distribution (see Fig. 6) peaks around one, with a long tail towards higher values, and we derived a probability of 85% for the transit to be grazing. The three parameters are strongly correlated. The radius distribution exhibits a long tail towards a high value and up to the boundary of our prior (i.e. 0.1) linked to the possible grazing transit, which is degenerate by nature (i.e. a highly grazing planet with a large radius may have the same transit shape as a smaller planet with a lower grazing transit).

Fig. 6 Corner plot for the photometry analysis (blue) and combined RV-photometry analysis (pink) with the parameters (radius, semi-major axis, and inclination) and the grazing criterion of planet c. A grazing criterion higher than one indicates that the transit is grazing. There is a 85% chance that the transit is grazing for planet c.

4.3 Photometric planetary parameters

The period of planet b is refined compared to Winters et al. (2019b) to 5.358764 ± 0.000004 days through the five additional transits we used in this analysis. The retrieved value for the orbital period of planet c is 3.123898 ± 0.000003. The complexity of the light curve shown by the sensitivity of the retrieved value for T₀ between the two sectors means that the ephemeris and properties derived for the planets are to be considered with caution. Transit-timing variations (TTV) cannot be entirely ruled out. However, the ratio of planet to star radii is marginally affected for planet b. The values derived for the two independent inferences and the combined light curve dataset are compatible within one sigma (Fig. 5). The uncertainty for planet c is higher because the transit is likely grazing, which increases the uncertainty on the radius.

5 Radial velocity analysis

5.1 Stellar activity

Radial velocity measurements were derived from the analysis of the host star spectroscopic absorption features. However, the shape and centroid position of these spectral lines are affected by many stellar effects that are commonly grouped under the term ‘stellar activity’. It encompasses a variety of physical processes related to the host star, which affects our capacity of unambiguously analysing the radial velocity signals arising from orbital motion. The activity signals operate on different timescales and amplitudes. They emerge as both incoherent and quasi-periodic variability modulated by the rotation period of the star (Queloz et al. 2001; Cloutier et al. 2017). With a typical stellar jitter magnitude of 1 to 3 m s⁻¹ for quiet stars (Dumusque et al. 2011) and up to 10 m s⁻¹ or higher for more active stars (Santos et al. 2014), it is paramount to investigate the activity effects of LTT 1445A to the best of our ability in order to claim reliable detections of low-mass planets. Unfortunately, these activity signals cannot always be reproduced by analytical models (Fischer et al. 2016; Wright 2018). Gaussian processes offer a flexible framework for modelling non-analytical physical processes such as stellar activity and for handling systematic and correlated noise in a robust way. It is broadly used in astrophysics in general (Lavie et al. 2017; Sale & Magorrian 2018; Barros et al. 2020) and within the radial velocity community in particular (Haywood et al. 2014; Rajpaul et al. 2015; Toledo-Padrón et al. 2020; Mortier et al. 2020; Faria et al. 2020).

5.1.1 Activity indices and Lomb-Scargle periodograms

Spectroscopic observations of LTT 1445A provide several activity indices to probe the stellar activity and rotation of the star. The rotation of LTT 1445A is not seen in the light curve, which indicates either low stellar activity or a rotational period that exceeds the time that is covered by the light-curve observations. The light curve only shows a modulation at 1.4 days, which has been explained by the rotation of one of the binary components of LTT 1445BC. We considered seven activity indices provided by the ESPRESSO pipeline: Hα, the S-index, contrast, NaD, the full width at half maximum (FWHM), Ca, and R_HK. As an initial step, the time series of each activity index were cleaned by removing outliers. Activity indicators sometimes correlate with velocities (Fischer et al. 2016). The radial velocities for LTT 1445A are correlated with the FWHM, are weakly correlated with the S-index and the NaD index, and are weakly anti-correlated with the contrast index (see Fig. A.7). The Lomb-Scargle periodograms (Lomb 1976; Scargle 1982) of all indices (see Fig. A.6) exhibit significant peaks at periods longer than 35 days. The periodograms were computed using the astropy routine (VanderPlas et al. 2012; Astropy Collaboration 2013, 2018; VanderPlas & Ivezic 2015). The Hα and R_HK indices also show a significant peak at a short period (between 1.6 and 2 days). The FWHM index displays the strongest peak in terms of significance at a period of 72 days.

5.1.2 Bayesian inference with a GP

We performed a Bayesian inference with a GP using the same kernel as in Sect. 4.1.1 independently on all the activity indices (three parameters for each activity index). As noted in the previous section, the choice of kernel was guided by our understanding of the underlying physics that we wished to model. The period hyperparameter controls the rotational period of LTT 1445A, hereafter P_rot. The length parameter L_act is related to the average lifetime of the active regions on the surface of the star. The spot patterns are stable for the duration of the rotational period for M-dwarf stars (Newton et al. 2018), therefore this parameter is expected to have a value equal or superior to the period of rotation. Initially, the inference was made on the logarithm of the hyperparameters with broad uniform priors between −10 and 10 units on all three parameters and without information provided by the periodogram analysis of Sect. 5.1.1. The hyperparameter space was explored with a nested sampling algorithm, as in Sect. 4. The GP converges towards low values for the period parameter and failed to catch the rotational period of the star as revealed with the periodogram analysis. The length parameter also converged to low values for all indices as the GP tried to catch any small variations in the dataset. In a second approach, we chose a more restrictive prior on the period hyperparameter by forcing it to be higher than 35 days. This reflects the information obtained through the periodogram analysis. We also included an inference on the combined time-series of all activity indices. The inference was made with seven amplitude parameters (one for each activity index) and the two other parameters of the kernel (L_act and P_rot). Unfortunately, none of these inferences converged to a specific value for the rotational period (Fig. A.7). Only the inference with the FWHM index exhibits a weak preference for a period of 72 days.

Finally, we tested models including the C parameter in the GP kernel (from Foreman-Mackey et al. 2017). These models are not favoured by the Bayesian analysis. We wish to stress that this result is specific to the LTT 1445A system and should not be considered as a general feature of stellar activity. The analysis of other systems may require the use of the C parameter in the GP kernel.

We conclude that it is not possible to constrain the hyperparameters of the GP with the activity indices for LTT 1445A. Nonetheless, we tried to make an educated guess with the periodogram analysis and our knowledge of stellar activity.

5.1.3 Expected stellar rotation for LTT 1445A

The rotation periods of M dwarfs span a broad range of values. Irwin et al. (2011) measured periods from 0.28 to 154 days. Using Kepler observations of M-dwarf stars between 0.3 and 0.55 M_⊙, McQuillan et al. (2013) reported two distinct groups with periods between 10–25 days and 25–80 days with peaks at ~19 and ~33 days, respectively. Giacobbe et al. (2020) reported rotation periods in the range of 0.52–191.82 days with a peak in distribution at ~30 days using the APACHE transit survey of M dwarfs in the mass range 0.15–0.70 M_⊙. The authors also reported a correlation of the rotation period of the star with its mass (when fast rotators with P < 10 days were excluded) and its age. Irwin et al. (2011) favoured the hypothesis of a rapid spin-down process as a function of age, which seems to be supported by the dearth of stars with intermediate rotation periods of between 10 and 70 days reported by Newton et al. (2016a,b, 2018). However, no clear function can help in mapping the rotation period of the star to its age. Estimating the age of stars also remains difficult, which implies systematic high uncertainties on the published values of age. Suárez Mascareño et al. (2018) proposed a simple equation (equation 1) for M dwarfs to map the rotational period of the star to the value of its . The measured for LTT 1445A is −5.88 ± 0.4, which indicates a rotational period longer than 100 days according this equation.

It is therefore difficult to estimate the rotation period of LTT 1445A. The periodogram analysis indicates that it is a slow rotator with a period longer than 70 days, which is compatible with the negligible rotational broadening detected by Winters et al. (2019b). We therefore chose to set the rotational period hyperparameter of the GP to the highest peak in the FWHM periodogram at 72 days. The length parameter L_act was fixed to 72 days as well in order to take into account that the active regions on the surface of M dwarfs are stable for one or several rotations (Newton et al. 2018). We note that the value for P_rot is compatible to what we may expect if the system is not too young (see Sect. 3). When we assume that the triple star system formed at the same time, a decrease similar to that of LTT 1445A in spin for the other two stars is inconsistent with the 1.4-day signal that results from the rotation of LTT 1445B or LTT 1445C.

5.2 Identification of radial velocity signals

In this subsection, we identify potential planets in the radial velocity time-series. The search for radial velocity signals is done without using any information from the light curve analysis performed in Sect. 4.

5.2.1 Radial velocity model for a blind search

We used only circular orbits, so that one planet was modelled by a Keplerian with only three parameters: the time of inferior conjunction T₀, the orbital period P, and the radial velocity semiamplitude K. The model also includes a mean centre-of-mass velocity for each instrument (γ_ESP and γ_HARPs) and two parameters to describe the movement of LTT 1445A due to the presence of the BC components: a slope and a curvature. The Keplerian model can be coupled with a GP with the kernel defined in the previous section. The joint GP-Keplerian models were implemented with the RADVEL package (Fulton et al. 2018), which makes use of the celerite package for the GP.

5.2.2 Iterations with periodogram analysis and Bayesian inference

We explored the different signals in the radial velocity time-series by subdividing the period parameter space into 200 bins of 0.5 days. We performed a Bayesian inference within each bin with the simple circular Keplerian (from one planet without GP to five planets with GP). One iteration corresponded to the inference made on all the period bins with one model. The prior on the orbital period was set to be uniform within each bin. The priors on the other parameters are given in Table 1. The multi-dimensional parameter space was explored with the nested sampling algorithm, which provided the Bayesian evidence that was used to classify the different periodic signals in the radial velocity time-series. In the first iteration, the strongest signals found with the one-planet model were at higher periods close to the possible rotation of the star at 72 days. For comparison with this period-subspace model ladder, we also computed the Lomb-Scargle periodogram using the ASTROPY implementation in Python (Astropy Collaboration 2013, 2018) for the ESPRESSO dataset alone, which confirmed that the main signal is compatible with the 72-day signal seen in the FWHM activity index. The second iteration was made with a joint GP-one Keplerian model. The hyperparameters L_act and P_rot of the GP were fixed to the value chosen in Sect. 5.1.2 (i.e. 72 days), only the amplitude was allowed to vary. In the next iteration, we fixed the parameters from the most significant signal in the previous iteration, added another planet, and proceeded with another inference on each bin. We iterated this process until no clear signal remained in the residue.

The subdivision of the orbital period parameter space and the strategy to fix the values of each signal before adding new Keplerian signals allowed a quick scan of the parameter space. In each iteration, the bulk of the posterior was found more easily with a reduced parameter space and fewer free parameters. The computed Bayesian evidence for each model (M1, M2, etc.) is a good proxy for the true value of these models with the larger period parameter space because the different bulks of the posterior of each significant signal are the main contribution to the evidence integral.

5.2.3 Identification of the different signals

The result of this procedure is summarised in Figs. 7 and A.8. The periodogram analysis yields results that are similar to those of the Bayesian inference for the first four iterations, but then failed to identify the 24-day signal. Figure 7 shows the significance of the signal detection within each orbital period bin for all the joint GP-Keplerian models with one to five planets. The significance is evaluated as the logarithm of the ratio of the Bayesian evidence (Bayes factor) of each pair inference (model/bin) with the best pair inference of the previous iteration. The Bayes factor for each bin with the joint GP-one Keplerian model (column M1) was computed with the lowest evidence value amongst the different period subdivision with the same model because any bin is strongly supported compared to the highest evidence obtained with the one-Keplerian model without GP.

Confirmation of the two transiting planets

The first two detected signals are at 5.36 and 3.13 days, which corresponds to the orbital periods obtained for planet b and c with the photometry analysis in Sect. 4. This is a direct confirmation that the two planets orbit the same host star LTT 1445A.

A complicated signal at 3.4 days

The third significant Keplerian signal is at a puzzling period of 3.4 days. A planet with an orbital period like this will have repercussions on the dynamical properties of the system and especially on the orbit of planet c with a potential TTV, which we do not see. The analysis of the light curve did not show any potential transits with a 3.4-day period, which implies a planet inclined compared to the 5.336- and 3.13-day planets. It therefore appears to be difficult to claim that this signal is due to the presence of another planet. We note that a 3.4-day period if an alias оf the 1.4 days in our time-series, which is the value of the modulation signal seen in the light curve that we have associated with either B or C. Therefore, we suspected that this signal was a contamination of some light from LTT 1445BC that entered ESPRESSO. However, the separation between LTT 1445A and the BC component is about 7″ (Mason et al. 2001; Dieterich et al. 2012) compared to the one-arcsecond diameter of the fiber of ESPRESSO (Pepe et al. 2014, 2021). The average seeing of Paranal is 0.77″, which will place the LTT 1445BC light source at several FWHM from the fiber. The contamination is therefore strongly unlikely. Nonetheless, we tested the hypothesis of a contamination trom the 1.4-day modulation by performing the same blind search with the two (GP kernel hyperparameters P_rot and L_act fixed to the value obtained with the photometry analysis. The long-period signals bound to the rotation period of the star are partially caught by this short-period GP because GPs with exponential decay components in their kernels destroy the covariance between data points that are separated by much more than the decay timescale, so they can model signals at lower frequencies than their exponential-decay timescales. The signal with the highest significance remains the orbital period off 5.36 days of planet b (see Fig. A.9). In the second iteration (M2 model), the signal of planet c at 3.13 days is detected only weakly. We added a Keplerian for planet c with the ephemeris obtained with the photometry analysis. The third signal is the 24-day signal presented below, which is also seen in the blind search with the 72-day GP kernel (it is the most significant signal for the inference with the M2 model). The other iterations with the more complex models M4 and M5 only show weak Keplerian signals. The 3.4 days disappeared completely. To test our hypothesis further, we fit a joint GP-Keplerian model at 1.4 days with the GP at 72 days, which also removed the 3.4-day signal without weakening the significance of the two transiting planets and the 24-day signal presented below. The origin of the 3.4-day signal remains unknown. A full dynamical analysis of the system is beyond the scope of this paper, but we stress that this analysis is necessary to reliably rule out the possibility of another small planet with an orbital period of about 3.4 days. We note that which hypothesis is chosen to include the signal in the models does not impact the parameters that are derived for the three other planets.

Fig. 7 Significance of the Keplerian signals for the blind search in the radial velocity time-series. The x-axis shows the joint GP-RV model we used for the inference from one Keplerian (M1) and up to five Keplerians (M5). The y-axis indicates the subdivision of the orbital period parameter space with a bin size of 0.5 days cut ait 30 days. One cell corresponds to a pair model or bin that indicase s the model we used for the inference and the prior on the orbital period of the plantet. The logarithm of the ratio of the Bayesian evidence (Bayes factor) between the pair model (Mx)/bin and the best pair model (Mx−1)/bin is colour-coded from weak (value of 0) to strong (value of 15 or above) to indicate the significance of the signal detected in the corresponding bin (except for the first column, in which the Bayes factor of each pair M1/bin is computed with the worst pair M1/worst bin). The orbital period obtained for the best signal of each iteration is written in black in the corresponding cell. An example: the cell corresponding to the M2 joint GP-two Keplerion model for the period subdivision of 3–3.5 days (the retrieved period is 3.13 days) indicates a strong significance of this signal compared to a GP-RV model with only one Keplerian at 5.36 days, which was the best signal detected with the previous iteration with M1

A possible planet at 24 days

The fourth significant Keplerian signal is at a 24-day orbital period with a semi-amplitude of ~1.48 m s⁻¹. It does not correspond to any peak in the periodogram analysis of the stellar activity indices (Fig. A.7). We searched for potential transits using the ephemeris provided by the inference on the radial velocity dataset, but did not find any signal. However, only one transit for each sector would be visible, with the risk that the transit for sector 4 falls within a gap in the light curve based on the radial velocity ephemeris. This low number of potential transits combined with the 1.4-day modulation signal complicates any detection of the planet in the light curves, so it is not possible to entirely rule out the possibility that a potential planet (with the ephemeris retrieved from the RV) does transit.

5.3 Radial velocity analysis

In this subsection, we perform a thorough analysis of the radial velocity time-series with different models and use the identification of the different Keplerian signals from the quick blind search. In this subsection, information from the light curves analysis are used as prior for the two inner planets detected with both observing methods.

5.3.1 Models and nomenclature

We considered a set of models that encapsulates different assumptions to model the Keplerian signals in the radial velocities.

Firstly, one model included a certain number of planets, indicate by the first number after the letter M (M1, M2 etc.). Each planet is modelled by a Keplerian with three parameter for circular orbits (the time of inferior conjunction T₀, the orbital period P, and the radial velocity semi-amplitude K) оr with five parameters for eccentric orbits (eccentricty e and longitude of periastron ω are added to the free parameters). The assumption for each planet is denoted by the letters following the notation of the model, a letter c or e per planet indicates circular or eccentric orbits, respectively.

Secondly, each model included a mean centre-of-mass velocity for each instrument (γ_ESP and γ_HARPS) and possibly a jitter term for each instrument (σ_ESP,σ_HARPS). These assumptions are denoted by the letter G when no jitter is added to the model and by a J when both the γ_S and σ_S are included in the model. The additional jitter terms aim at modelling any source of white noise in the radial velocities that will not be caught by the joint GP-model or the default uncertainties from the pipeline. The jitter terms are added in quadrature to the formal uncertainties derived by the ESPRESSO and HARPS pipelines.

Thirdly, three assumptions were considered for the motion of LTT 1445A due to the presence of the BC components: a slope, a slope and a curvature, or no movement (no additional parameters). The slope is denoted in the model name with the letter L and the curvature by the letter C. When the name contains no letter, then the model does not include any parameters for the movement of LTT 1445A. It is very unlikely that models without a slope are favoured by the model selection because the trend is clearly visible by eye in the radial velocity time-series. We tested it regardlessly to evaluate the significance of the two other assumptions.

Fourthly, the use of a joint GP-RV model is indicated by the notation gp1 in the name of the model. We used the kernel defined in Sect. 4.1 with the two parameters P_rot and L_rot fixed to 72 days, as discussed in Sect. 5.1.3. Only the amplitude of the GP is a free parameter.

Finally, we chose to model the 3.4-day signal in the radial velocities with either a Keplerian with an orbital period of 3.4 days (the third signal in models with at least three Keplerian) or 1.4 days (the third signal in models with at least three Keplerian and a superscript ^1.4 in the name of the model), or with a GP with its hyperparameters fixed to the photometry values obtained for the 1.4-day modulation signal in the light curve. Only the amplitude of the GP is a free parameter. When the model uses this assumption, the model name includes gp2 and does not possess the 72-day GP.

5.3.2 Bayesian inference and model selection

We performed the Bayesian inference with each model on the radial velocity time-series using the nested sampling algorithm with the same setup as in the previous sections (1000 live points, 0.8 sampling efficiency, and at least three multiple runs). The first two planets in all the multi-planet models were planet b and planet c, which were also detected in the photometry analysis. The priors on the time of inferior conjunction and the orbital periods of these planets were Gaussian, with the centre and the standard deviation set to the retrieved values of the inference on the light curve (Sect. 4.2.2). The constraints on the ephemeris obtained with the photometry analysis are superiors to the one obtained with the radial velocities alone. For the other signals, the prior on the orbital period and the time of inferior conjunction were uniform in a range that encompassed the period obtained with the blind search. The priors were all uniform on the other parameters and were the same as in the combined analysis (Sect. 6) given in Table 1.

The first output of our analysis is the model ladder. It classifies the significance of all models using Bayesian evidence. The full table with all the models is shown in Table A.1, and a reduced table is presented in Table 2, in which the best ten models are broken down between the two hypothesies on the 3.4-day radial velocity signal (a GP or a Keplerian). Our threshold for weak, moderate, strong, and high significance of a model are with logarithms of the Bayes factor of 2, 3, 5, and 11, respectively (Trotta 2008). A model with more complexity (i.e. more parameters) needs at least an increase in logarithm of the evidence of two compared to a more simple model (i.e. with fewer parameters) to be moderately significant and favoured. From this analysis, we can derive the following conclusion for the LTT 1445A bcd planetary system:

GP models are favoured

Models with a GP(with period of 1.4 or 72 days) are strongly favoured because they absorb the main significant signals in the radial velocities at long periods that we assume to correspond to the stellar rotation. We point out that our analysis is skewed towards planets with small periodicity (i.e. less than 40 days) because of the time span of our dataset.

Planets b and c strongly favoured

The existence of the two planets b and c is strongly favoured by the model selection. Models with eccentric orbits (MXee* with X ⩾ 2) have no or a weak significance compared to the respective model with circular orbits for both planets (MXcc*), which indicates that the two additional parameters of the eccentric models add complexity that is not required to explain the data. More observations (or more accurate data) should help to constrain the eccentricities better.

3.4 day signal

It is necessary to model the 3.4-day signal either by a Keplerian orbit (models M4 with gp1) or by the 1.4-day GP (models M3 or M4 with gp2). Models with a Keplerian at 3.4 or 1.4 days have a higher evidence than models with the 1.4 days GP because the long-term variations are better handled by the 72-day GP.

Long-term variation

The model selection favours models with jitter terms (letter J) and a slope (letter L) in the case of the 72-day GP. Jitter terms are not favoured in the case of the 1.4-day GP because this latter absorbs any short-period jitter more efficiently than the 72-day GP. The curvature parameter is not favoured. This could indicate that LTT 1445A is in the linear regime of its orbit with the components BC or that the curvature is too small and is entirely dominated by the short-term variations (stellar jitter and planets).

A new planet candidate

The 24-day signal is moderately favoured by the model selection when the 3.4-day signal is modelled with a Keplerian (the logarithm of the Bayes factor between M4ccecJLgp1 and M3cceJLgp1 is 2.09) and strongly favoured in the case of the 1.4-day GP models (ln 𝒵_M3cccGLgp2 – ln 𝒵_M2ccGLgp2 = 7.92). The detection under the other assumptions (slope, curvature, and jitter) is in approximately the same range of significance. Therefore, we claim the existence of planet d orbiting LTT1445 A with a period of 24 days with a moderate significance. We advocate for more observations to bolster the significance of the detection.

The retrieved values for the semi-amplitude of the three planets are robust as they are stable to the different assumptions, except for models that failed to explain the data (with Bayesian evidence smaller than minus 300). Figure 8 shows the different values obtained for the best models with the different thresholds of the model selection. Finally, we ran retrievals with the four Keplerian signals, the slope and jitter (model M4ddcdJL), but with a Gaussian prior on the two other hyperparameters P_rot and L_act of the GP. We tested two Gaussian priors: 72 ± 10 days on both P_rot and L_act based on the FWHM value, and 130 ± 30 days based on the derived value from Suárez Mascareño et al. (2018). We also performed a retrieval with a GP including the C parameter as a free parameter and the same setup as above. The retrieved distributions for the Keplerian parameters are not affected by these new hypotheses.

Fig. 8 Radial velocity semi-amplitude of planets b (black), c (blue), and d (green) as a function of the logarithm of the Bayesian evidence obtained with the inference on the radial velocity dataset alone. The higher the Bayesian evidence, the more favoured the model. The vertical dashed red line indicates the different thresholds in the difference with the logarithm of the Bayesian evidence of the best model. Models with a difference larger than 5 are strongly rejected.

6 Combined RV-photometry analysis

As a final step of our analysis, we performed a joint inference on the photometric and spectroscopic datasets to derive self-consistent uncertainties taking the correlations between the different parameters of the light curve and radial velocity models into account (similar to Sozzetti et al. 2021). We used the best radial velocity model for the circular (M4ccecJLgp1) and eccentric (M4eeeeJLgp1) orbits with the respective photometric model with both transiting planets. We only used the 2-min cadence for both sectors for the photometric dataset. The results are consistent with the individual radial velocity and TESS analysis. Figures 9 and 10 show the final fit for the radial velocity and light curve observations. Physical properties of the planetary system are shown in Table 3. Table 1 shows the values for LTT1445A planetary system retrieved with this joint inference. As a summary, we confirm the discovery of planet LTT 1445A b (Winters et al. 2019b) with a refined orbital period of 5.358764 ± 0.000004 days around component A of the system of three M dwarfs LTT 1445ABC at 6.9 pc. We present the discovery of LTT 1445A c with a period of 3.123898 ± 0.000003 days and the candidate planet LTT 1445A d at days. The two inner planets are transiting, but there is a high probability that the transit of planet c is grazing. We derived radii of 1.43 ± 0.09 R_⊕ and R_⊕ for planets b and c, respectively. We found no evidence that LTT 1445A d transits, but the time window covered by the two TESS sectors and the strong impact of the 1.4-day modulation signal on the light curves prevent us from ruling this possibility out. The masses derived from this analysis are 2.32 ± 0.25 M_⊕, 1.00 ± 0.19 M_⊕ for planets LTT 1445A b and c, respectively. We derived a minimum mass for the potential planet d of 2.72 ± 0.75 M_⊕. Furthermore, LTT1445A exhibits a long-term drift. We derived a radial velocity slope of −4.4 m s⁻¹ yr⁻¹, which agrees with the maximum expected value if it is due to the gravitational of BC, which according to Torres (1999), given the total mass of BC provided by Winters et al. (2019b), the system distance and present-day separation of around 7 , is 9 m s⁻¹ yr⁻¹. Furthermore, Kervella et al. (2019) reported a statistically significant (almost 20 sigma) proper motion anomaly at the mean epoch of Gaia DR2 for LTT 1445 A. Finally, Gaia EDR3 lists low RUWE values (1.07) and astrometric excess noise (158microarcsec). All the points above confirm that the motion of the primary is due to the outer BC components, and a closer-in (within about 5 au) massive gas giant or brown dwarf companion can be ruled out.

Fig. 9 Phase-folded radial velocities for planets LTT 1445A b (top panel), c (second panel), and d (third panel). The red squares are HARPS data, the orange circles are ESPRESSO data, and the red circles are the ESPRESSO data binned in 10 bins. The blue lines indicate the best fit obtained with the combined photometric and radial velocity datasets. The two last panels show the residuals for HARPS observations (second to bottom panel) and ESPRESSO (bottom panel).

Fig. 10 Phase-folded light curves with both sectors at a 2-min cadence for planet b (top panel) and planet c (bottom panel). The black dots show the flux for the 2-min cadence detrended with the GP from the best-fit model obtained with the combined photometric and radial velocity datasets. The white dots with small error bars are the same flux binned. The red line shows the transit model component of the best fit.

7 Discussion

Figure 11 shows that planet b is compatible within two sigma with an Earth-like composition (Dorn et al. 2015). The planet is less dense than Earth, with a density of 4.36 ± 0.74 g cm⁻³, which might indicate the presence of an atmosphere. Alternatively, the planet may have lower iron-mass fractions because the host star is metal poor (Adibekyan et al. 2021; Santos et al. 2017). Follow-up observations focusing on the detection and characterisation of the atmosphere of the planet are needed to explore these different hypotheses. The uncertainties on the radius of planet c are high because the transit may be grazing. The planet may be compatible with a water-rich composition, but any real statement on the nature of the planet seems premature. The planet is also compatible at two sigma with an Earth-like composition. Future radial velocity observations combined with photometric observations for example with CHEOPS (Benz et al. 2021) should reduce the uncertainties to allow a clearer picture of the density of this planet. The amount of water is expected to be low for planet b in comparison to the refractory materials (e.g. iron or silicates). Because the two planets are on compact orbits, if planet c were water-rich, it would require a mechanism to bring water to it without enhancing b or a mechanism to deplete water on planet b without depleting the water on planet c. Alternatively, planet c may have formed beyond its current orbit far away in the disc and then migrated to its current position. In that case, its water content may be very different from that of planet b. This seems unlikely, however, as any migration of planet c would have consequences on planet b as well. In any case, the presence or absence of an atmosphere cannot be unambiguously claimed for either planet. As discussed by Winters et al. (2019b), it is unlikely that we will find better targets than LTT 1445A bc based on the known occurrence rates of planets around M dwarfs. The system is therefore a prime target for an atmospheric characterisation with space telescopes such as the JWST. Winters et al. (2019b) used Kopparapu et al. (2013, 2014) to compute the conservative inner and outer habitable zone (HZ) boundaries for LTT 1445A and found 0.093 and 0.182 AU, respectively, which places planets b and c far away from these limits. With a semimajor axis of 0.09 AU, the super-Earth LTT 1445Ad may sit within the inner boundary of the habitable zone. The habitability around M dwarfs has been the subject of extensive debate, in particular, with regard to the prolonged stellar activity or the gravitational effects that result in tidally locked configuration for planets within the habitable zone (see Shields et al. 2016 and references within for an overview). The habitability of planets in multiple stellar system also brings its share of difficulties, especially in terms of orbit stability. Many studies have been devoted to the habitability in binary systems (Jones et al. 2001; Noble et al. 2002; Takeda et al. 2008; Dvorak et al. 2010). In principle, these studies can also be expanded to systems of higher order (Cuntz 2014). The interest in the habitability of LTT 1445A d is limited with the transit technique, if indeed it does not transit. However, the planet may be a prime target for next-generation instruments such as the ELT-HIRES with its adaptive optics that may allow us to characterise a possible atmosphere.

Of the more than 3500 planetary systems known thus far, 111 are planetary systems around binary stars, and 26 lie in multiple star systems (Schwarz et al. 2016). Unfortunately, we only know a handful of multi-planetary systems around triple or higher-order star systems so far. A radial velocity detection of up to seven planets has been claimed for the M dwarf GJ 667 C (Anglada-Escudé et al. 2012; Gregory 2012). The star is a companion of the K3V+K5V binary GJ 667AB at a minimum separation of 230 au. Two Jovian planets with eccentric orbits have been detected by radial velocity around the main component of two triple systems, HD 207832 A (Haghighipour et al. 2012; Ngo et al. 2017) and HD 65216 A (Mayor et al. 2004, 2011). Both hosts are G-type stars and possess two companions at projected separations of 110 and 253 AU, respectively. Finally, two planetary systems with radial velocities and transit detection have been reported for the F star K2-290 (Hjorth et al. 2019, 2021) and the K star Kepler-444 (Campante et al. 2015). The first system possesses a close-by companion at 113 AU and a far-away companion at more than 2450 AU. The two planets are a hot mini-Neptune and a warm Jupiter. Kepler-444 possesses two unresolved M-dwarf companions at a separation of 66 AU (Lillo-Box et al. 2014; Dupuy et al. 2016) with very eccentric orbits that bring these two companions very close to the main component (within 5 AU). Kepler-444 is also known to be the oldest planetary system with five compact (within 0.1 au) small planets. In this context, the report of the LTT 1445 A bcd planetary system is a good addition to this small sample, but is not sufficient to derive robust statistics regarding multi-planet systems around triple or higher-order star systems. We note that LTT 1445 is the only known system composed of three M dwarfs. The architecture is similar to that of Kepler-444, HD 207832, and HD 65216 in the sense that the planets have been detected around the main component of the star system and that the host possesses two companions that orbit each other closely (6 AU for HD 65216 BC, 4 AU for HD 207832 BC, and unresolved for Kepler-444 BC). Current literature on exoplanets in multiple star systems focuses on binaries because on one hand, they are the main systems among the multiple stellar systems, and on the other hand, many more planets have been detected around them than around higher-order systems. It appears from observations that the dynamical properties, in terms of eccentricity, mass, and period distribution of planets orbiting binaries are not different than those of planets around single stars (Ngo et al. 2017) in the case of giant planets between 0.1 and 5 AU, but extrapolation to smaller planets is not trivial. The presence of a companion or a binary companion as in the case of LTT 1445 A has strong implications for planet formation mechanisms as the protoplanetary disc may be truncated (Artymowicz & Lubow 1994) and conditions will be less favourable for planet formation (Nelson 2000; Mayer et al. 2005). Stellar companions may also eject any formed planets (Zuckerman 2014). The three planets LTT 1445A bcd are S-type planets because they orbit one component of the stellar multiple system (Dvorak 1986). The system is sufficiently compact (within 0.2 AU) to be dynamically stable when the work done on binaries is extrapolated (Holman & Wiegert 1999; Marzari & Gallina 2016; De Cesare & Capuzzo-Dolcetta 2021) to LTT 1445A -BC and a separation of 21 AU between the two components is assumed (Winters et al. 2019b). The stability of planet orbits around binaries is strongly dependent on the eccentricity of the companion. The orbit of LTT 1445 BC is unknown, and we cannot rule out the possibility that a high eccentricity will bring the two companions closer to the host star. However, stars may retain their planets if the system is compact enough even if companions are on very eccentric orbits. For example, Kepler-444 possesses a binary companion on a very eccentric orbit that regularly brings it to ~5 AU of the primary star. Nevertheless, five planets were found to orbit the primary, and recent study showed the poor dynamical perturbative impact of the stellar binary on the stability of the planetary orbits (Stalport et al. 2022). The entire LTT 1445 system may be coplanar (Winters et al. 2019b), and the stability of the planetary system may also reflect some particular aspect of the formation history of LTT 1445ABC.

To conclude, the LTT 1445 system offers great opportunities for studying a multi-planet system in a triple star system with two terrestrial planets that transit and a super-Earth in the habitable zone. We advocate for spectroscopic and photometric follow-up observations of the system to confirm the status of planet d, to prepare an atmospheric characterisation of the two inner planets with the next generation of telescopes, and to monitor the presence or absence of other planets in the system.

Fig. 11 Mass-radius diagram including the planets selected by Otegi et al. (2020). The two inner planets LTT 1445Ab and LTT 1445Ac are indicated in red.

Acknowledgements

We are grateful for a constructive review from Sarah Blunt, which allowed us to deepen our analysis. The authors acknowledge the ESPRESSO project team for its effort and dedication in building the ESPRESSO instrument. This work has been carried out in the frame of the National Centre for Competence in Research ‘PlanetS’ supported by the Swiss National Science Foundation (SNSF). B.L., A.D. acknowledge the financial support of the SNSF. A.D. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project FOUR ACES, grant agreement No. 724427). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. A.S.M., J.I.G.H., C.A.P. and R.R. acknowledge financial support from the Spanish Ministry of Science and Innovation (MICINN) project PID2020-117493GB-I00. A.S.M., J.I.G.H. and R.R. also acknowledge financial support from the Government of the Canary Islands project ProID2020010129. V.A. acknowledges the support from FCT through Investigador FCT contract no. IF/00650/2015/CP1273/CT0001. J.I.G.H. also acknowledges financial support from the Spanish MICINN under 2013 Ramón y Cajal program RYC-2013-14875. N.J.N. acknowledges support from the following projects: UIDB/04434/2020 & UIDP/04434/2020, CERN/FIS-PAR/0037/2019, PTDC/FIS-OUT/29048/2017, COMPETE2020: POCI-01-0145-FEDER-028987 & FCT: PTDC/FIS-AST/28987/2017. A.S.M. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) under the 2018 Juan de la Cierva Programme IJC2018-035229-I. IAC: This work is partly financed by the Spanish Ministry of Economics and Competitiveness through grants PGC2018-098153-B-C31.

Note added in proof

An independent analysis of this system using a different dataset was performed and published in Winters et al. (2022).

Appendix A Figures and detailed model ladder

Fig. A.1 Posterior distribution of the GP hyperparameters for the 1.4-day signal in the TESS light curves from sector 4 (black), sector 31 with 2-minute cadence (orange), sector 31 with 20-second cadence (pink), and the combined light-curves dataset (blue). Priors are uniform between [−20,20] m.s−1 for the logarithm of the amplitude, [0,10] days for the logarithm of the length, and [1,2] days for the period.

Fig. A.2 Identification of planet candidates for TESS light curves in sector 4 (first column), sector 31 with 20-second cadence (second column), and sector 31 with 2-minute cadence (third column). The first row shows the PDCSAP light curves with the GP fit (see section 4.1). The second row shows the detrended light curves. The last panels show the iterative periodogram analysis. The detected signals are at 5.35 days (third row) and 3.12 days (fourth row). The remaining time series has no signal (last row).

Fig. A.3 Light curves centred around each transit of planet b. The plain blue line indicates the best fit from the Bayesian inference done on the combined light curves of both sectors with the joint GP-transit model. The dashed blue line shows the best fit from the Bayesian inference done with a circular transit model on the reduced dataset (two hours before and after mid-transits) detrended with a linear model.

Fig. A.4 Same as Figure A.3, but for planet c

Fig. A.5 Radial velocity measurements of LTT 1445A with HARPS (red) and ESPRESSO (orange). There is a clear long-term drift due to the interaction of the star with the BC binary of the triple system.

Fig. A.6 Time series, Lomb-Scargle periodograms, and correlation with radial velocities for each activity index of the ESPRESSO dataset for LTT 1445A. For the middle column, the dotted, dashed, and plain horizontal lines indicate the levels for a false-alarm probability of 0.1%, 1%, and 10%, respectively. For the last column, the Pearson correlation coefficient is computed for each pair (activity index and RV). It is indicated as a number and with a face colour in each subplot. A coefficient of 1 (−1) indicates that the activity indey is correlated (anti-correlated) with the radial velocities. A coefficient of 0 indicates no correlation.

Fig. A.7 Posterior distributions for the hyperparameters of the GP for the inference on the seven activity indicators of ESPRESSO (first seven rows) independently and for the combined inference with the seven indices (last row). The dotted blue histograms are from the broad uniform priors, and the black histograms are from the inference with a priorlonger than 35 days for the Prot parameter.

Fig. A.8 Blind search for Keplerian signals in the radial velocity dataset. The period parameter space is subdivided into bins of 0.5 days. An iteration corresponds to a Bayesian inference for each bin with a Keplerian model (first row) or a joint GP-Keplerian model with one to five planets (rows 2 to 6). The first column shows the logarithm of the Bayes factor of each model/period bin with the best model/period bin of the iteration. The second column shows the Lomb-Scargle periodogram computed with the residue of the best model from the previous iteration (except for the first row, which is the periodogram from the data). The last column shows the best fit to the data from the best model/period of the iteration, except for the first row, which shows the GP fit to the ESPRESSO dataset. The dashed vertical lines for columns 1 and 2 indicate the three Keplerian signals at 5.4, 3.1, and 24.3 days, the dotted vertical line indicates the 72-day period associated with the stellar rotation, and the dotted-dashed vertical line indicates the polluting signal at 3.4 days. The horizontal dashed line in column 1 indicates the moderate threshold for the model selection. The horizontal dashed line in column 2 indicates the level for a false-alarm probability of 1%.

Fig. A.9 Same as Figure 7, but with the 1.4 days GP.

References

All Tables

All Figures

Fig. 1 TPF of LTT 1445 from the TESS observations in Sector 4 (composed with tpfplotter; Aller et al. 2020). The SPOC pipeline aperture is overplotted with shaded red squares, and the Gaia DR3 catalogue is also overlaid with symbol sizes proportional to the magnitude contrast with the target, marked with a white cross. We display all sources from the Gaia DR3 catalogue with a magnitude contrast up to Δm = 8 mag.

In the text

	Fig. 2 Observed photometric spectral energy distribution of LTT 1445A (red dots) from the и filter through the W4 band. The best BT-Settl model (Allard et al. 2012), corresponding to Teff = 3300 K, low metallicity, and log g = 6.0 (cgs), is shown with the blue line.
In the text

Fig. 3 Location of LTT 1445A and BC in the Gaia colour-magnitude diagram. The binary BC is not resolved by Gaia. The young stellar isochrones of the β Pic (blue) and Tucana-Horologium (green) moving groups and the Pleiades cluster (purple) are taken from Luhman (2018), and the main sequence of field stars (black) is taken from Cifuentes et al. (2020). The grey area represents the dispersion observed among main-sequence stars.

In the text

Fig. 4 LTT1445 ligt curves from TESS sectors 4 and 31. First row: TESS PDCSAP light curve for sector 4 with 2-min cadence (first column) and sector 31 (second column) for a 2-min cadence (black) and a 20-s cadence (grey). A clear modulation with a period of 1.4 days is clearly visible. Second row: Respective detrended light curves with the best-fit model (in red) from the analysis of Sect. 4.2.2. The vertical lines identify the transit for each planet: b1 to b10 for planet b, and c1 to c16 for planet c. The light blue dots in the 2-min cadence in the upper panels are the outliers in the dataset that were removed from the PDCSAP light curves.

In the text

Fig. 5 Posterior distributions for the Bayesian inference with the joint GP-transit model (plain lines) and linear model (dashed lines) for the sector 4 light curve (black), for sector 31 with the 2-min cadence (orange), and on the combined 2-min cadence light-curve dataset (blue). From top to bottom: time of inferior conjunction, period, ratio of the planet to star radius and impact parameter for planet b (first column) and planet c (second column). Inference on sector 31 shows a broad distribution and a 1-h time transit variation and shift in the orbital period. This broad distribution seems to result from the higher amplitude of the polluting signal at 1.4 days for sector 31. The inference on the combined light curve dataset converges towards the same values as the inference on sector 4.

In the text

	Fig. 6 Corner plot for the photometry analysis (blue) and combined RV-photometry analysis (pink) with the parameters (radius, semi-major axis, and inclination) and the grazing criterion of planet c. A grazing criterion higher than one indicates that the transit is grazing. There is a 85% chance that the transit is grazing for planet c.
In the text

Fig. 7 Significance of the Keplerian signals for the blind search in the radial velocity time-series. The x-axis shows the joint GP-RV model we used for the inference from one Keplerian (M1) and up to five Keplerians (M5). The y-axis indicates the subdivision of the orbital period parameter space with a bin size of 0.5 days cut ait 30 days. One cell corresponds to a pair model or bin that indicase s the model we used for the inference and the prior on the orbital period of the plantet. The logarithm of the ratio of the Bayesian evidence (Bayes factor) between the pair model (Mx)/bin and the best pair model (Mx−1)/bin is colour-coded from weak (value of 0) to strong (value of 15 or above) to indicate the significance of the signal detected in the corresponding bin (except for the first column, in which the Bayes factor of each pair M1/bin is computed with the worst pair M1/worst bin). The orbital period obtained for the best signal of each iteration is written in black in the corresponding cell. An example: the cell corresponding to the M2 joint GP-two Keplerion model for the period subdivision of 3–3.5 days (the retrieved period is 3.13 days) indicates a strong significance of this signal compared to a GP-RV model with only one Keplerian at 5.36 days, which was the best signal detected with the previous iteration with M1

In the text

Fig. 8 Radial velocity semi-amplitude of planets b (black), c (blue), and d (green) as a function of the logarithm of the Bayesian evidence obtained with the inference on the radial velocity dataset alone. The higher the Bayesian evidence, the more favoured the model. The vertical dashed red line indicates the different thresholds in the difference with the logarithm of the Bayesian evidence of the best model. Models with a difference larger than 5 are strongly rejected.

In the text

Fig. 9 Phase-folded radial velocities for planets LTT 1445A b (top panel), c (second panel), and d (third panel). The red squares are HARPS data, the orange circles are ESPRESSO data, and the red circles are the ESPRESSO data binned in 10 bins. The blue lines indicate the best fit obtained with the combined photometric and radial velocity datasets. The two last panels show the residuals for HARPS observations (second to bottom panel) and ESPRESSO (bottom panel).

In the text

Fig. 10 Phase-folded light curves with both sectors at a 2-min cadence for planet b (top panel) and planet c (bottom panel). The black dots show the flux for the 2-min cadence detrended with the GP from the best-fit model obtained with the combined photometric and radial velocity datasets. The white dots with small error bars are the same flux binned. The red line shows the transit model component of the best fit.

In the text

	Fig. 11 Mass-radius diagram including the planets selected by Otegi et al. (2020). The two inner planets LTT 1445Ab and LTT 1445Ac are indicated in red.
In the text

Fig. A.1 Posterior distribution of the GP hyperparameters for the 1.4-day signal in the TESS light curves from sector 4 (black), sector 31 with 2-minute cadence (orange), sector 31 with 20-second cadence (pink), and the combined light-curves dataset (blue). Priors are uniform between [−20,20] m.s−1 for the logarithm of the amplitude, [0,10] days for the logarithm of the length, and [1,2] days for the period.

In the text

Fig. A.2 Identification of planet candidates for TESS light curves in sector 4 (first column), sector 31 with 20-second cadence (second column), and sector 31 with 2-minute cadence (third column). The first row shows the PDCSAP light curves with the GP fit (see section 4.1). The second row shows the detrended light curves. The last panels show the iterative periodogram analysis. The detected signals are at 5.35 days (third row) and 3.12 days (fourth row). The remaining time series has no signal (last row).

In the text

Fig. A.3 Light curves centred around each transit of planet b. The plain blue line indicates the best fit from the Bayesian inference done on the combined light curves of both sectors with the joint GP-transit model. The dashed blue line shows the best fit from the Bayesian inference done with a circular transit model on the reduced dataset (two hours before and after mid-transits) detrended with a linear model.

In the text

	Fig. A.4 Same as Figure A.3, but for planet c
In the text

	Fig. A.5 Radial velocity measurements of LTT 1445A with HARPS (red) and ESPRESSO (orange). There is a clear long-term drift due to the interaction of the star with the BC binary of the triple system.
In the text

Fig. A.6 Time series, Lomb-Scargle periodograms, and correlation with radial velocities for each activity index of the ESPRESSO dataset for LTT 1445A. For the middle column, the dotted, dashed, and plain horizontal lines indicate the levels for a false-alarm probability of 0.1%, 1%, and 10%, respectively. For the last column, the Pearson correlation coefficient is computed for each pair (activity index and RV). It is indicated as a number and with a face colour in each subplot. A coefficient of 1 (−1) indicates that the activity indey is correlated (anti-correlated) with the radial velocities. A coefficient of 0 indicates no correlation.

In the text

	Fig. A.7 Posterior distributions for the hyperparameters of the GP for the inference on the seven activity indicators of ESPRESSO (first seven rows) independently and for the combined inference with the seven indices (last row). The dotted blue histograms are from the broad uniform priors, and the black histograms are from the inference with a priorlonger than 35 days for the Prot parameter.
In the text

Fig. A.8 Blind search for Keplerian signals in the radial velocity dataset. The period parameter space is subdivided into bins of 0.5 days. An iteration corresponds to a Bayesian inference for each bin with a Keplerian model (first row) or a joint GP-Keplerian model with one to five planets (rows 2 to 6). The first column shows the logarithm of the Bayes factor of each model/period bin with the best model/period bin of the iteration. The second column shows the Lomb-Scargle periodogram computed with the residue of the best model from the previous iteration (except for the first row, which is the periodogram from the data). The last column shows the best fit to the data from the best model/period of the iteration, except for the first row, which shows the GP fit to the ESPRESSO dataset. The dashed vertical lines for columns 1 and 2 indicate the three Keplerian signals at 5.4, 3.1, and 24.3 days, the dotted vertical line indicates the 72-day period associated with the stellar rotation, and the dotted-dashed vertical line indicates the polluting signal at 3.4 days. The horizontal dashed line in column 1 indicates the moderate threshold for the model selection. The horizontal dashed line in column 2 indicates the level for a false-alarm probability of 1%.

In the text

	Fig. A.9 Same as Figure 7, but with the 1.4 days GP.
In the text