© ESO 2020

1 Introduction

In the past decades, the exploration of exoplanets has moved from the detection scheme to the characterization challenge. The new dedicated ground- and space-based facilities built for this purpose now offer the possibility of fully characterizing the properties of extrasolar planets with exquisite precision (e.g., Pepe et al. 2014; Suárez Mascareño et al. 2020; Damasso et al. 2020). Understanding planet formation and evolution, determining atmospheric properties, and ultimately searching for biosignatures requires a sufficiently large sample of planetary systems. The plethora of properties of the current exoplanet population (more than 4100 known so far, according to the NASA Exoplanet Archive, Akeson et al. 2013) indeed offers a large collection of targets to carry out this endeavor. The combinationof different techniques then becomes critical to provide a complete view of the system, which in turn allows the inference of its history through feeding population synthesis models and allowing subsequent atmospheric characterization campaigns. In particular, the combination of the transits and radial velocity techniques is key to understanding the bulk composition of the planets. Exquisite photometric and radial velocity precision is needed to infer its internal structure, however, and key to setting observational constraints on formation and evolution processes.

Planets around low-mass stars have drawn the attention of the community in the recent years due to the ground-based instrumental capabilities in reaching the rocky domain (e.g., Luque et al. 2019; Zechmeister et al. 2019), also in the temperate region around the star (TRAPPIST-1 Gillon et al. 2016). The LHS 1140 planetary system is one of these examples. This M4.5 dwarf is located 10.5 parsec away (Gaia Collaboration 2018) and hosts two known planets, a small short-period telluric component (LHS 1140 c, P_c ~ 3.77 days) and a temperate rocky super-Earth (LHS 1140 b, P_b ~ 24.7 days). LHS 1140 b was detected by Dittmann et al. (2017) using MEarth (Nutzman & Charbonneau 2008) photometric time series and HARPS (Mayor et al. 2003) radial velocities. Additional HARPS measurements and observations from the Spitzer Space Telescope by Ment et al. (2019) led to the detection of the inner component in the system, LHS 1140 c, which also transits its host star.

LHS 1140 b lies within the habitable zone of the star and possesses a rocky composition, thus representing a key target for further astrobiological studies. The path toward a habitability analysis involves, among others, the study of the internal structure of the planet (Shahar et al. 2019). Its composition and distribution is a direct consequence of its formation (e.g., the abundances of the different elemental building blocks) and evolution (i.e., the different heating and cooling processes, and impacts throughout the planet history), see Dorn et al. (2018) and references there in. The two planets in LHS 1140 with their different orbital (hence irradiation) properties offer a unique opportunity to understand different evolutionary paths in the same environment.

Added to this, LHS 1140 is also an ideal target for searches of co-orbital planets. Co-orbital configurations consist of planet pairs trapped in gravitationally stable regions and 1:1 mean motion resonances. The Lagrangian points L₄ /L₅ present in the gravitational field of a two-body system (like a star and a planet) are stable points of equilibrium that are located exactly on the same orbit as the planet, but ± 60° ahead and behind it. These gravitational wells have been demonstrated to be very stable once an object is trapped (Laughlin & Chambers 2002). The only condition for stability, once trapped, is that the total mass of the planet pair (planet plus trojan) must be lower than 3.8% of the mass of the star (e.g., Gascheau 1843). This relaxed constraint allows similar-mass planets to co-orbit together in a long-term stable dance around the star. The two planets can share the same orbital path in different co-orbital configurations. In the case of circular or quasi-circular orbits, the possibilities are tadpole and horseshoe orbits. In the first case, the co-orbital surrounds (or librates around) one of the Lagrangian points in the co-rotating frame with the planet, as in the case of the Jupiter trojans. In the horseshoe configuration, the co-orbital describes a horseshoe shape in the co-rotating frame, moving from L₄ to L₅ and exchanging orbits with the planet. This is the case, for instance, of the Saturnian moons Janus and Epimetheus, which have comparable sizes. In the case of eccentric orbits, other configurations come into play, such as the quasi-satellite or the anti-Lagrange configuration (Giuppone et al. 2010). In general, several studies have shown that the formation and stability of co-orbital planets makes these configurations not only possible but probable (Cresswell & Nelson 2008, 2009; Leleu et al. 2019a) and that their dynamical properties could have hidden them in the noise of already existing data (Ford & Holman 2007; Madhusudhan & Winn 2009; Janson 2013; Hippke & Angerhausen 2015). However, despite the ample room for stability estimated from theoretical studies, no co-orbital pairs have been found so far. Several attempts have been focused on these configurations and have so far set observational constraints to their existence in different regimes (see, e.g., Lillo-Box et al. 2018a,b), and one key candidate has already been found (TOI-178, Leleu et al. 2019a).

The LHS 1140 planetary system has important properties that make it a key system in which to search for these co-orbitals. First, the low mass of the host star allows reaching the sub-Earth mass domain. The transiting nature of the two known planets allows the application of co-orbital detection techniques, thus avoiding the degeneracy with the eccentricity at first order (Leleu et al. 2017). The edge-on orientation of the planetary system also permits the search for co-orbitals through the transit technique (Janson 2013; Lillo-Box et al. 2018b). The slow rotational velocity of the star allows precise radial velocity measurements. Moreover, the multi-planet nature of the system increases the likelihood that it hosts co-orbital pairs (Leleu et al. 2019b), especially when the planets are in mean motion resonances (Cresswell & Nelson 2008; Leleu et al. 2019b).

We present the joint analysis of the first observations of this system obtained with the Echelle SPectrograph for Rocky Exoplanetsand Stable Spectroscopic Observations (ESPRESSO) instrument and the precise light curve obtained with the Transiting Exoplanets Survey Satellite (TESS) mission. The data are presented in Sect. 2. In Sect. 3 we explore the new radial velocities in search for other planets in the system. The TESS light curve is also further explored in Sect. 4. In Sect. 5 we explore the possibility of co-orbitals to the two known planets in the system from different perspectives and with different techniques. The final joint data analysis including the radial velocity and the light curve is presented in Sect. 6. We discuss the result in Sect. 7 and provide the final conclusions in Sect. 8.

2 Observations

2.1 HARPS

LHS 1140 was intensively observed by with the High Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003) at the 3.6 m telescope of the European Southern Observatory (ESO) La Silla facilities under the program IDs 191.C-0873 and 198.C-0838 (PI: X. Bonfils) and 0100.C-0884 (PI: N. Astudillo-Defru). In total, this radial velocity dataset comprised 293 HARPS radial velocity measurements that were presented in Ment et al. (2019), spanning 783 days between November 2015 (JD = 2 457 349.65) and January 2018 (JD = 2 458 133.54). The averagecadence was one epoch every 1.8 days, but normally, two spectra were obtained every night with separations of a few hours. We performed a night binning of the dataset, as suggested in Dumusque et al. (2012). This simple strategy has proved to be efficient to reduce the effect of short-term correlated noises on orbital elements (e.g., Hara et al. 2019). In total, 145 individual measurements were available. The corresponding uncertainties per binned data point are distributed around 3.2 ±1.2 m s⁻¹. The periodogram of this dataset is presented in Fig. 1.

Fig. 1 Periodogram of the radial velocity (upper panel) for each dataset individually (ESPRESSO in light green and HARPS in red) and for the joint dataset of both instruments assuming a 26 m s−1 offset (see Sect. 6). The activity indicators for the individual datasets are also shown in the lower panels. The periods of the three planets are marked as dotted vertical lines, and the rotation period (Prot = 131 days) and its first harmonic (Prot∕2 = 65 days) are marked as vertical dashed lines.

2.2 ESPRESSO

ESPRESSO is the new ultra-stable high-resolution spectrograph of the Very Large Telescope at ESO’s Paranal Observatory(Pepe et al. 2020). This facility has the capability of collecting light from any of the four 8.2 m Unit Telescopes (UTs). The instrument has two arms with a total wavelength coverage of 380–788 nm. ESPRESSO is equipped with a powerful reduction pipeline that performs all steps of the basic reduction and provides high-level data products, including radial velocities. The aim is a final radial velocity precision down to the 10 cm s⁻¹ level.

We obtained 116 spectra of LHS 1140 in three consecutive semesters¹ (P103, P104, and P105). The total time span of the observations is 404 days; they were taken between October 2018 (JD = 2 458 416.71) and December 2019 (JD = 2 458 820.56), and a mean cadence of one spectrum was obtained every 2.2 days (median of one spectrum per night). Each spectrum was obtained with an exposure time of 1820 s in P102 and 1915 s in P103 and P104, and produced a median signal-to-noise ratio (S/N) S∕N = 80 at 700 nm. We processed the whole dataset using version v2.0.0 of the ESPRESSO Data reduction Software² (DRS) pipeline (Pepe et al. 2020). We selected a binary mask corresponding to an M5 spectral type for this M4.5 star to perform the cross-correlation and obtain the final radial velocities (Baranne et al. 1996; Pepe et al. 2003). The resulting radial velocities have a mean photon-noise uncertainty of 1.07 m s⁻¹. In June 2019, ESPRESSO suffered a major intervention, and the chamber and the vacuum vessel were opened to perform a fiber-link exchange. This introduced a small jump in the radial velocity datasets taken before and after the intervention, and the actual value of this jump is not constant for the different spectral types. The data taken before (labeled ESPRESSOpre in this paper) and after (ESPRESSOpost) should therefore be treated as coming from different instruments. This fiber-link exchange occurred between our first and second observation semester. The mean precision corresponding to the data obtained before the fiber-link exchange is 1.16 m s⁻¹, and the post-change data have an average uncertainty of 1.02 m s⁻¹. This improvement goes in line with the increase in transmission measured after the intervention (Pepe et al. 2020) and a slight increase in exposure time from 1820 s in P102 to 1915 s in P103 and P104. We removed two data points with uncertainties above 3 m s⁻¹ because low clouds caused a low S/N (the observations were aborted after half the requested exposure time). The final dataset including the radial velocities, activity indicators, S/N, and exposure time for each spectrum is presented in Table B.1. Figure 2 shows the time-series radial velocity for both HARPS and ESPRESSO, and Fig. 3 shows the phase-folded curves for the two known planets in the system.

The periodogram of the full dataset (including HARPS and ESPRESSO measurements, see Fig. 1) shows clear signals of the known planets LHS 1140 b and LHS 1140 c. The two signals are detectable independently in the two instrument datasets (ESPRESSO and HARPS), and their significance is boosted when they are combined by only including a 26 m s⁻¹ radial velocity offset between the two (see Sect. 3.2). In addition to the signals of these two planets, other signals appear at a significant level. The rotation period of the star at P_rot = 131 days is clearly detected in each of the two datasets and in the combined set. Interestingly, both datasets also show a peak at half the rotation period P_rot∕2 ~ 65 days, which is notas evident in the HARPS radial velocity dataset. However, the rotation period and its alias are clearly visible in the activity index corresponding to the full width at half maximum (FWHM) of the cross-correlation function (CCF) in both datasets, see Fig. 1. Added to this, additional signals stand out in the region around 70–100 days. This is further investigated in Sect. 3.2.

Taking advantage of the large number of ESPRESSO spectra, we used ODUSSEAS (Antoniadis-Karnavas et al. 2020) to estimate theeffective temperature and metallicity of this M star. This machine-learning tool measures pseudo-equivalent widths for more than 4000 spectral lines and compares them to a training dataset composed of HARPS spectra of reference M-dwarfstars. Although the code was originally tested for spectra with resolutions from 48 000 to 115 000, we used the combined ESPRESSO spectrum directly with a higher resolution of 140 000. Although the resolution is slightly higher than the highest resolution grid of the code (115 000), we can safely assume that the method that measures the pseudo-equivalent widths provides compatible values when a spectrum with higher resolution is used³. With this code we derived an effective temperature of T_eff = 2988 ± 67 K and a metallicity of [Fe/H] = −0.262 ± 0.104 dex. These errorswere estimated by considering on the one hand the machine-learning model error (precision errors of 17 K and 0.03 dex, respectively), and on the other hand the mean absolute errors of the machine-learning models for the training dataset (65 K and 0.10 dex for a resolution of 115 000 for the dataset grid, Antoniadis-Karnavas et al. 2020).

Fig. 2 Top: radial velocity of LHS 1140 from the HARPS (red) and ESPRESSO (open for ESPRESSOpre and filled for ESPRESSOpost in green) datasets. The black line shows the median radial velocity model from the joint photometric and radial velocity analysis including the two known Keplerian signals and the GP model (see Sect. 6). The gray shaded regions correspond to the 68.7% (dark gray) and 95% (light gray) confidence intervals of the model. The median GP model is shown as a dashed blue line. Middle: radial velocity dataset after removing the median GP model. The Keplerian model is shown as a solid black line. Bottom: radial velocity residuals of the full model.

Fig. 3 Phase-folded radial velocity signal of LHS 1140 b (left panels) and LHS 1140 c (right panels). For each planet, upper panel: radial velocity data (color-coded as in the legend), binned radial velocities (open black symbols) with a bin size corresponding to 15% of the phase, the median model for each planet obtained from the joint fit analysis (see Sect. 6) as the solid black line, and the 68.7 and 95% confidence intervals as shaded dark and light gray regions. Bottom panels: residuals of the median model.

Fig. 4 TPF of LHS 1140 from the TESS observations in Sector 3 (composed with tpfplotter, Aller et al. 2020). The SPOC pipeline aperture is overplotted with shaded red squares, and the Gaia DR2 catalog is also overlaid with symbol sizes proportional to the magnitude contrast with the target, marked with a white cross.

2.3 TESS photometry

The Transiting Exoplanet Sky Survey (TESS, Ricker et al. 2014) observed LHS 1140 during Sector 3 in camera 1, from 20 September 2018 (JD = 2 459 115) to 18 October 2018 (JD = 2 459 141), immediately before our ESPRESSO campaign. We used tpfplotter⁴ (Aller et al. 2020) to check for contaminant sources in the automatically selected aperture. This is shown in Fig. 4, where we display all sources from the Gaia DR2 catalog with magnitude contrast up to Δm = 8 mag. In addition to LHS 1140, two additional sources lie inside the TESS aperture with magnitude contrasts in the Gaia passband of 3.8 and 4.0 mag. This imposes an upper limit to the dilution factor of 5.2%.

We used the light curve extracted by the SPOC pipeline. For the purpose of this paper, we use the presearch data-conditioning simple aperture photometry (PSDCSAP) detrending of the data, which has a CDPP of 0.547 parts per thousand (hereafter ppt). In Fig. 5 we show the extracted photometric time series and the phase-folded light curves for the two known planets. The photometric precision of this time series makes it sufficient for an independent analysis and characterization of the two planets. We therefore did not use the previous observations from Spitzer and MEarth presented in Ment et al. (2019).

3 Exploring the radial velocity dataset

3.1 ℓ₁ periodogram

We first analyzed the radial velocity data with the ℓ₁ -periodogram, as defined in Hara et al. (2017). This tool is designed to search for a representation of the signal as a sum of a small number of sinusoids, where “small” is compared to the number of observations. It has a similar aspect to a regular periodogram, but with far fewer peaks due to aliasing. As in Hara et al. (2020), we computed the ℓ₁ -periodogram of the data with different assumptions on the noise covariance. The covariance models were then ranked by cross-validation. We considered the ESPRESSO and HARPS data without binning.

Similarly to Haywood et al. (2014), we included in the model two activity indicators per instrument, smoothed with Gaussian kernels with different timescales, as linear predictors. These indicators were chosen because they exhibited significant variations on a short (0–2 days) and longer timescale (order of days and 10–20 days). We included the FWHM, with a smoothing timescale of 20 days and 2 days for HARPS and ESPRESSO, respectively, the pipeline-derived line asymmetry for ESPRESSO (10 days) and bisector span for HARPS (1 day). The timescales were chosen after fitting the hyperparameters of the Gaussian kernel. We further added a quadratic trend to the model. We assumed a noise model with a white-noise component σ_W, a calibration error σ_C, and a Gaussian component with amplitude σ_R and exponential decay τ. We then considered all possible combinations of values for σ_R = 0.0,0.5, 1.0, 1.5, 2 m s⁻¹, σ_W = 0.5, 1 m s⁻¹ σ_C = 0.1, 0.5, 0.75 m s⁻¹, and τ = 0.0, 1.0, 3.0, 6.0 d. All these noise models were ranked with cross validation, as in Hara et al. (2020). In Fig. 6 we present the ℓ₁ periodogram obtained with the noise model with maximum cross-validation score (corresponding to σ_W = 1 m s⁻¹, σ_R = 0.5 m s⁻¹, τ = 3 d, and σ_C = 0.75 m s⁻¹.).

We then considered the 20% highest ranked model CV₂₀, and computed the number of times that a signal is included in the model (which reaches a false-alarm probability, FAP < 0.05). We find that signals might be included in the models at seven periods. In decreasing order of significance, these are at 24.7, 3.77, 78.8, 62, 424, 129, and 16.9 days (see Table 1). The last two signals are included only in 50% of the CV₂₀ models, and16.9 days does not appear in the highest ranked noise model. The other signals appear to be present in the data. The signals at 3.77 and 24.7 days correspond to the known planets LHS 1140 c and LHS 1140 b, respectively. The signals at 130 days and 62 days are likely linked to the rotation period of the star (P_rot) and half of it, P_rot∕2. The origin of the 433 days signal is less certain because its significance also depends on the way the stellar indicators are included in the model. The signal at 78.8 days (the third most significant in the ℓ1-periodogram) might be due to a planet and hence deserves additional attention.

3.2 Evidence for a planet with a period of ~ 80 days

As reported in Ment et al. (2019) and confirmed in the new ESPRESSO dataset (see Sects. 2.2 and 3.1), the periodogram of the joint dataset shows a power excess between the first and second harmonic of the rotational period of the star; this is between 65 and 131 days. These several signals have a maximum peak at around 80 days. We here explorethe evidence for this third signal by studying different scenarios (different number of planets, from one planet to three planets; and different orbital configurations, eccentric or circular) and comparing the models using Bayesian analysis to unveil the significance of this signal.

We explored each dataset separately (i.e., HARPS and ESPRESSO), together with a final full radial velocity dataset analysis (i.e., HARPS+ESPRESSO). Each planet signal was modeled as a Keplerian function, with independent values for the radial velocity semiamplitude of each planet (K_i), orbital period (P_i), eccentricity (e_i), and argument of the periastron (ω_i). Because the clear effect of the stellar activity on the radial velocity is shown in the periodogram and has been pointed out by previous studies, we used Gaussian processes (GPs) to model the correlated noise. We used the quasi-periodic kernel from the george (Ambikasaran et al. 2016) implementation, which is a combination of an exponential decay and a periodic part, (1)

Here η₁ corresponds to the amplitude of the correlated noise, η₂ can be interpreted as the timescale of the variations of the stellar features causing the correlated noise (see, e.g., Faria et al. 2016), η₃ represents the stellar rotation period (P_rot), and η₄ is a balance between the exponential and the periodic components of the kernel. As shown in Suárez Mascareño et al. (2020), we can additionally use the FWHM of the cross-correlation function as an activity indicator to further constrain the GPs (see the FWHM time series in Fig. 7). In this case, the radial velocities and the FWHM share all hyperparameters of the GP, except for the amplitude, therefore an additional parameter per instrument (η_1,FWHM,j) is included. For ESPRESSO, which has the highest resolution of the planet searchers and the most stable instrumental profile, the FWHM is expected to be the best tracker of line-deforming stellar activity. This indicator therefore provides an excellent way to determine the radial velocity variations that are caused by this stellar activity. For every instrument, we also added a jitter value for the radial velocity and the FWHM (σ_RV,j and σ_FWHM,j, respectively) to account for additional uncorrelated noise and a systemic offset (Δ_RV,j and Δ_FWHM,j).

We used the emcee⁵ (Foreman-Mackey et al. 2013) implementation of the Goodman & Weare (2010) affine invariant Markov chain Monte Carlo (MCMC) ensemble sampler to explore the parameter space and sample the posterior distribution of each parameter. For each model, we used a number of walkers equal to four times the number of parameters involved (N_p) and 20 000 steps for each walker. Following the recommendations from the george package, we first explored the parameter space with a full run using all the steps and walkers. In a second phase, we resampled the walker positions around the best walker and ran a second iteration with the same number of steps to improve convergence. A final run of the maximum probability region was produced by using 10 000 steps per walker. Convergence was checked by estimating the autocorrelation time (τ) of the chains and ensuring a minimum length of the chain, corresponding to 30 × τ. The final joint chain was then composed of 40 000 × N_p, which for the simplest model (no planets) corresponds to 6.8 × 10⁵ steps and for the most complex model (including three Keplerians and three instruments) corresponds to 1.3 × 10⁶ steps.

The priors for each parameter involved in the different models are shown in Table B.2. In brief, we assumed Gaussian priors on the planet periods and mid-transit times, centered on the values from Ment et al. (2019), but with a broad width of five times the published uncertainties. The remaining Keplerian parameters were set to broad uniform priors throughout the entire allowed regime. The prior on the GP hyperparameter η₃ was set to a Gaussian centered on the rotational period of the star, corresponding to 131 ± 5 days (Dittmann et al. 2017). Based on the observed radial velocities, we set the η₁ and η_1,FWHM parameters to a uniform distribution with a maximum amplitude of 100 m s⁻¹. The timescale hyperparameter η₂ was set to a broad range of values between one and five times the stellar rotation period. The instrument systemic velocities were set up with uniform priors between −13.4 and −13.0 km s⁻¹, and the radial velocity jitter for each instrument was allowed to have values up to 20 m s⁻¹ (based on our first model attempts).

In total,15 models (including zero to three planets and all combinations of circular and eccentric orbits for each of them) were tested for each of the three datasets (HARPS, ESPRESSO, and HARPS+ESPRESSO), that is, 45 models were tested in total. In order to statistically compare the different models and datasets, we estimated the Bayesian evidence () of each model by using the perrakis⁶ implementation. Based on this Bayesian evidence, we can estimate the Bayes factor () between two models as the ratio between the evidence of each pair of models as . The model with the strongest evidence has the highest statistical relevance, with considered as a strong evidence for one model against the other.

The three panels in Fig. A.1 display this Bayesian evidence for all 15 models considered in each dataset. We find that the two-planet model with both planets on circular orbits has the strongest evidence in all datasets. This means that the current data are unable to provide evidence for eccentric architectures of the two known planets. We can only place upper limits of e_b < 0.096 and e_c < 0.274 (at 95% confidence level) based on the noncircular models. The median and 68.7% confidence interval of the parameters for this two-planet model are presented in Table. B.2.

Although the evidence in the full dataset is weaker, all three-planet models in the combined HARPS+ESPRESSO dataset converged to a planet candidate with a period of about 80 days, corresponding to the peak in the ℓ1-periodogram described in Sect. 3.1 and already seen in the classical periodogram presented in Sect. 2.2 and Fig. 1. The three-planet model with all planets on circular orbits (labeled 3p1c2c3c) has the second strongest evidence, even above the eccentric two-planet cases. The median and 68.7% confidence interval of the parameters for this circular three-planet model are presented in the last column of Table. B.2, and the phase-folded radial velocity curve of this third planet candidate is shown in Fig. 8. It is important to note that this strong evidence is mainly due to the addition of the exquisite ESPRESSO data, which itself show clear signs for this three-planet model.

We performed a cumulative Bayesian analysis of the dataset by obtaining the Bayes factor for the two most likely models with two and three planets (i.e., the circular case for all planets, labeled 2p1c2c and 3p1c2c3c in Fig. A.1). We obtained the log-evidence using up to N_i data points in steps of five data points and starting with the first 20 measurements. The results are presented in Fig. 9 and show that the log-evidence of the three-planet model, although the evidence is weaker than for the two-planet model, increases progressively when more data points are added. This is an indication that additional data might clearly confirm this signal. The jumps in the Bayes factor correspond to the addition of new data from a different instrument, which correspondingly adds new parameters (and so complexity) to the model. The evolution of the Bayesian evidence shows that about 50 new ESPRESSO measurements would be needed to clearly confirm this candidate signal.

According to our three-planet model in the combined dataset, the third planet would have a minimum mass of M_⊕ and an orbital period of days. This describes a world in the rocky-gaseous frontier beyond the habitable zone. We note here that the detection of its radial velocity signal is at the ~ 3.8σ level (). According to the ephemeris found in this radial velocity solution, the TESS observations would have missed the transit of this small planet (in case of a coplanar orbit with the other known planets in the system, the estimated impact parameter is 0.65), which would have occurred about five days before TESS observations started in this sector. The estimated radius of this planet according to forecaster (Chen & Kipping 2017) is R_⊕, which placesthe posterior distribution for the planet radius in the middle of the radius valley. The corresponding transit depth would be about 8 ppt, which is comparable with the transit depth of LHS 1140 b and therefore easily detectable by TESS and ground-based instrumentation. TESS will revisit the system in Sector 30 (22 September 2020 to 21 October 2020, in cycle 3). However, given the derived ephemeris for the candidate planet LHS 1140 d, the probability of a transit within the TESS observations in Sector 30 is unfortunately only 0.03% (see the analysis in Sect. 6).

To verify the results, we performed an independent analysis with kima (Faria et al. 2018) on the HARPS+ESPRESSO dataset. We considered up to three Keplerian signals in addition to a GP model for the correlated noise. Broad and equal priors were assigned to the three sets of orbital parameters (i.e., we did not use information from the transits). For the orbital periods we used a log-uniform distribution between 1 and 100 days, for the semiamplitudes a uniform prior between 0 and 20 m s⁻¹, and for the eccentricities a Kumaraswamy distribution with shape parameters a = 0.867 and b = 3.03 (Kipping 2013). The priors for the GP and remaining parameters were otherwise very similar to those in Table B.2.

The results indicate a significant detection of only two signals, corresponding to the orbital periods of the two known transiting planets, and the orbital parameters are fully consistent with those found previously. The probability ratio between the three-planetand two-planet models is estimated at 1.15, leading to the marginal detection of a third signal, with an orbital period of 78.9 ± 0.5 days. This is compatible with the results explained above. As before, the data do not constrain the eccentricities in the two-planet model. Consistent results were also obtained for the GP parameters.

Fig. 5 TESS light curve extracted from the SPOC pipeline and detrended using the PDCSAP flux. Upper panel: complete light curve displaying the transits of the two known transiting planets in the system, LHS 1140 b (red, one transit) and LHS 1140 c (light green, five transits). Lower panel: phase-folded light curve centered on the planet phase for each of the two transiting planets, LHS 1140 b (bottom left) and LHS 1140 c (bottom right). The colored symbols correspond to bins of one fifth of the transit duration (i.e., 30 min for LHS 1140 b and 11 min for LHS 1140c).

Fig. 6 ℓ1 -periodogram of the LHS 1140 HARPS and ESPRESSO data with the noise model that has the best cross-validation score.

Fig. 7 Dataset of the FWHM of the cross-correlation function from the HARPS and ESPRESSO data. The systemic offset for each instrument as retrieved from the two-planet model has been removed. The median GP model is shown as the solid black line, and the 1σ and 3σ confidence intervals are shown as shaded dark and light gray regions, respectively. Bottom panel: residuals after removing the GP model.

Fig. 8 Phase-folded radial velocity signal of the candidate planet LHS 1140 d in the three-planet model. Upper panel: radial velocity data (color-coded as in the legend), binned radial velocities (open black symbols) with a bin size corresponding to 15% of the phase, the median model from the joint three-planet analysis (see Sect. 6) as the solid black line, and the 68.7 and 95% confidence intervals as shaded dark and light gray regions. Bottom panels: residuals of the median model.

Fig. 9 Difference in the log-evidence of the two-planet and three-planet models with circular orbits for a cumulative number of data points. Top panel: Bayes factor against the time span of the data set, and bottom panel: against the number of data points.

3.3 Limits on additional planets at different periods

Because no additional planets can be confirmed with the current dataset, we explored the sensitivity of the HARPS and ESPRESSO data and sampling by injecting planetary signals (we assumed circular orbits) for different periods and masses ranging from 0.3–1000 days and 0.1–300 M_⊕. In total, 10⁴ signals were injected in the GP-removed dataset⁷. Then we proceed in the same manner as described in Sect. 3.2 (but we now removed the GP part) to recover the injected planet signal. We considered the planet to be detected when the median of the posterior distribution for the radial velocity semiamplitude was more than five times the standard deviation away for zero. This analysis is shown in Fig. 10, where we also add the detection limit assuming the scatter of the radial velocity data from the two-planet model in Sect. 3.2. We conclude that we can discard planets with masses m > 1 M_⊕ for periods up to 10 days, m > 2 M_⊕ for periods below 100 days, and m > 5 M_⊕ for periods shorter than one year. We note that these limits refer to the detection of the planet signal. The statistical significance of the model with respect to the two-planet model is not assessed in this step. The signal from the planet candidate LHS 1140 d is therefore above the sensitivity line in Fig. 10 (we can detect it in the data), but it cannot be confirmed with sufficient statistical significance.

Fig. 10 Detection limits of the radial velocity dataset for periods shorter than one year. The color code indicates the significance of the radial velocity semiamplitude parameter. The 5σ contour is shown as a solid line, and the limit corresponding to the scatter of the data is shown as a dashed line. The locations of the two known planets are shown as solid circles, and the location of the third planet candidate is shown as an open circle.

4 Exploring the TESS dataset

The full analysis of the transit signal of LHS 1140 b and LHS 1140 c is presented in Sect. 6. Here, we explore the possibility of additional signals in the TESS dataset.

We first analyzed the TESS dataset to search for transit-timing variations (TTVs) of LHS 1140 c, for which five transits are present in the data. We used allesfitter (Günther & Daylan 2020) to model the photometricdata alone. The priors on the parameters were set to Gaussian distributions according to the values published inMent et al. (2019), with a broad width corresponding to ten times the published uncertainties. We also setuniform priors to the individual TTVs of each transit of ± 3.5 h. The results for the inferred TTVs are shown in Table 2. All values are compatible with zero TTVs, except for the first transit, which occurs slightly earlier. However, the TTV is still compatible with zero within 2σ. Additionally, the TESS transit times are also compatible within one minute with the Spitzer transit time obtained six months before by Ment et al. (2019).

We used the transit least-squares (TLS) software (tls, see Heller et al. 2019) to search for additional signals in the TESS dataset. We masked out the times of transit corresponding to planets LHS 1140 b and LHS 1140 c and performed a smooth detrending of the light curve using wōtan (Hippke et al. 2019) with a window length of 0.5 and using the biweight method. We then computed the TLS periodogram on this masked and detrended light curve over the period range calculated by the tls software based on the data sampling with a range spanning from 0.6 to10 days. The result is shown in Fig. 11. The highest peak in the TLS periodogram corresponds to a periodicity of 1.26 days at a signal detection efficiency (SDE) of 6. Other intriguing signals are also present: a signal close to the orbital period of LHS 1140 c. Although the signal does not reach the statistical significance level of SDE > 7 required to accept a periodicity as statistically valid (see Heller et al. 2019), the fact that it appears at a similar orbital period as LHS 1140 c is intriguing and is further analyzed in Sect. 5 and in particular in Sect. 5.2.

Fig. 11 TLS periodogram of the TESS dataset after masking the light curve from the transit times corresponding to the two known planets LHS 1140 b and LHS 1140 c. The red line indicates the strongest power peak, and the blue line indicates the period of LHS 1140 c.

5 Co-orbital analysis

The long time span and high precision of the radial velocity dataset obtained for this target allows us to perform a detailed exploration for the possible presence of co-orbital planets (exotrojans). We explored this scenario for both planets in the system by applying the technique described in Leleu et al. (2017) (hereafter the α-test), a generalization of the technique proposed by Ford & Gaudi (2006). The α-test has been applied to different planetary systems, for instance, by Lillo-Box et al. (2018a,b), Armstrong et al. (2020), or Toledo-Padrón et al. (2020). In this technique, the radial velocity induced by the co-orbital system corresponds the sum of two Keplerians with the same orbital period, with at least one of the components transiting its host star (assumed to be the planet). When we also assume that the mass of the planet pair (m_p + m_t) is far lower than the stellar mass (M_⋆) and that in the near-circular case (typically, e < 0.1), the radial velocity can be approximated to first order in eccentricity and mass by (2)

where n = 2π∕P_orb, γ is the systemic velocity, and is the radial velocity semiamplitude of the co-orbital pair. The most relevant parameter in this equation is α, which to first order in eccentricity and trojan-to-planet mass ratio is (3)

where ζ is the difference in the mean longitude of the two components of the co-orbital pair, and t is the time, with the origin set to the mid-transit time of the main planet. This α parameter is therefore zero when no co-orbital is present (as m_t = 0) and can be qualitatively described as the mass ratio between the trojan and the planet. If the radial velocity data are compatible with a nonzero value for this parameter, we can infer a potential mass imbalance between the planet and the locations of the Lagrangian points. A negative value corresponds to L₄ and a positive value to L₅.

5.1 Co-orbital analysis for LHS 1140 b

The main interest in this planet resides on its location (in the middle of the habitable zone of this M dwarf) and its physical properties, including a rocky composition and low eccentricity (if any). These properties make this planet ideal for co-orbital searches and add the interest that it might lie in the habitable zone. Dvorak et al. (2004) have demonstrated that rocky trojans that co-orbit with planets in the habitable zone can also be habitable.

According to the results presented in Sect. 3.2, we here assumed only planets LHS 1140 b and LHS 1140 c in the system. We therefore added another Keplerian signal to Eq. (2) to account for the radial velocity contribution of LHS 1140 c. We followed the same modeling procedure as in Sect. 3.2. We again tested all possible scenarios corresponding to the two-planet models in that section (i.e., we assumed circular or eccentric orbits for the two planets).

We find the strongest evidence for the circular model for both planets LHS 1140 b and LHS 1140 c. In this model the inferred α parameter is α = 0.015 ± 0.065. The estimated value is fully compatible with zero, clearly stating that no co-orbital is detectable down to our sensitivity limits. Instead, using the posterior distribution of this parameter, we can set limits to the presence of co-orbitals to this planet at both Lagrangian points. Assuming the 95% confidence level as an upper limit, we can discard co-orbital planets to LHS 1140 b more massive than 1.8 M_⊕ at L₄ and 1.0 M_⊕ at L₅. Using only the HARPS measurements, we obtained a broader distribution of , which provided upper limits of 2.1 M_⊕ at L₄ and 1.8 M_⊕ at L₅. By adding the 116 ESPRESSO data points, we can decrease this limits by > 30%.

We also inspected the TESS light curve around the location of the Lagrangian points L₄ /L₅ of LHS 1140 b. Unfortunately, the L₄ location falls into a gap in the middle of the TESS sector. The transit of L₅ region occurs at BTJD = 1404.05, about eight hours before a transit of the inner planet LHS 1140 c. At the exact location of the Lagrangian point passage, no transit is found within the photometric precision. The scatter of this region in the light curve is 2.1 ppt (parts per thousand), which would correspond to a 3σ upper limit on any object corresponding to 3 R_⊕.

5.2 Co-orbital analysis for LHS 1140 c

We proceeded in the same manner as for LHS 1140 b, but now included the α parameter for planet c. In this case, we obtain a value of . The 95% confidence interval is between α ∈ [−0.27, 0.02]. This means that the 95% upper limit (corresponding to 2σ) is only marginally compatible with zero. A further exploration is therefore valuable. Because of its negative value, the mass imbalance producing this radial velocity signal can be interpreted as an existing body located at L₄. Assuming the mass of LHS 1140 c as calculated in Sect. 3.2 and a coplanar orbit, the mass of the trojan⁸ would correspond to 0.26± 0.18 M_⊕. Despite this ~2σ signal, the Bayesian evidence of the co-orbital model is still lower than the two-planet model. For the HARPS-only dataset, the evidence is stable throughout the entire dataset. However, in the case of the ESPRESSO-only dataset, the evidence for the co-orbital model progressively (although slowly) increases toward positive values as more data points are added. It therefore deserves additional attention.

We also performed an independent test to confirm this signal by following the technique proposed by Ford & Gaudi (2006) to compare the time of conjunction measured from the transit signal (T0_c,LC) and from the radial velocity signal (T0_c,RV). In the case of a single planet in circular orbit, it is easy to see that the lag between both times is zero, Δt = T0_c,LC − T0_c,RV = 0. However, in the presence of a co-orbital body at one of the Lagrangian points, the gravitational pull of the co-orbital will change T0_c,RV (this occurs earlier if in L₄ and later if in L₅), while leaving the time of transit (T0_c,LC) unchanged. We performed this test by analyzing both datasets independently. We modeled the TESS light curve with two planets and the radial velocity dataset with two Keplerian signals (details of the modeling of the photometric dataset are provided in Sect. 6). As a result, we obtain a time lag between the times of conjunction for LHS 1140 c measured from both techniques of Δt = 2.04 ± 1.27 h. This value is again 1.6σ away from zero and provides consistent results with the α-test, an indication for a co-orbital body in L₄ of LHS 1140 c. We note that this time lag can also be explained by a single planet with e cosω ~ 0.15, which is within the confidence interval of the eccentric two-planet model we showed in Sect. 3.2.

We further investigated this co-orbital candidate by inspecting the TESS light curve in this Lagrangian point. Five transits are observed with TESS of LHS 1140 c, but only four passages of its L₄ in front of the star because the third Lagrangian point lies inside the mid-sector gap. We explored the TESS light curve in the region of the Lagrangian point after removing the data points inside the transit of planet LHS 1140 b which occurs close to the L₄ location of transit 4 of LHS 1140 c. The phase-folded light curve in this regime shows a shallow dimming at the location of the L₄ Lagrangian point (slightly shifted by 0.006 in phase, or equivalently, ~32 min, but perfectly compatible with a dynamically stable libration orbit), see the middle panel in Fig. 12. Interestingly, the duration of the dimming is perfectly compatible with the transit duration of the planet LHS 1140 c, expected for a co-orbital in the coplanar case. The depth of this dimming is at the same level (~ 370 ppm) as the photometric scatter in the out-of-transit regions, however. We can compare this with the L₅ region. The photometric scatter in L₅ is significantly smaller (~280 ppm), and no signs of dimming are seen at the L₄ depth. This may indicate that the origin of the dimming is related to the larger systematics around the L₄ region. However, the different indications (including the radial velocity offset, the duration of the dimming corresponding to the planet transit duration, and its location close to the Lagrangian point) encouraged us to further investigate this signal. We therefore modeled it together with the signal of the two known planets to retrieve the size of the potential co-orbital candidate. The result provides a radius for this candidate of 0.44 ± 0.08 R_⊕. This result increases the certainty in the co-orbital scenario. However, the shallow depth (compatible with the photometric noise in the out-of-transit regime around L₄) prevents us from definitively confirming the nature of this dimming. If confirmed, it might also explain the peak in the TLS periodogram close to the period of LHS 1140 c described in Sect. 4 (see also Fig. 11). Additional data from the TESS extended mission will shed more light on this signal.

Fig. 12 TESS phase-folded light curve around the Lagrangian points L5 (top panel) and L4 (middle panel) of LHS 1140 c (see Sect. 5.2). The expected location of the transit is marked by vertical dashed lines. The median transit model inferred from the analysis of the L4 dimming is shown as a solid red line (middle panel) and a dashed line in the upper panel to guide the eye on the photometric scatter of the L5 region. The 1σ confidence interval is shown in the middle panel as a shaded red region. Bottom panel: residuals of the model for the L4 region. All panels show bin data points corresponding to 6.6 min (light blue, ten data points inside the transit duration) and 13.2 min (black, five data points inside the transit duration).

6 Joint photometric and radial velocity model

Based on the analysis performed in the previous section, we proceeded with a joint simultaneous analysis of the HARPS and ESPRESSO radial velocities and the TESS photometric time series to estimate the physical and orbital parameters of the system. To this end,we modeled the radial velocities as explained in Sect. 3.2, but with a slightly different parameterization that now explicitly included the individual contributions of the planet mass and orbital inclination for each planet (instead of using the radial velocity semiamplitude). We also included the GPs in the modeling and constrained them by simultaneously modeling the FWHM of the CCF of the ESPRESSO and HARPS data (see Sect. 3.2).

The TESS photometry was modeled with the batman code (Kreidberg 2015) to retrieve the transit models. We assumed a quadratic limb-darkening model with Gaussian priors around the values calculated by using the limb-darkening code⁹ (Espinoza & Jordán 2015) for the stellar parameters of LHS 1140 published in Ment et al. (2019). We also added photometric jitter to account for underestimated white noise, a mean level parameter to account for imperfect normalization of the light curve, and a dilution factor to account for the contamination in the TESS light curve due to the additional sources in the aperture (see Fig. 4). The stellar mass and radius were also included in the fit to properly account for their uncertainties by using a Gaussian prior around the published values by Ment et al. (2019).

To sample the posterior distributions, we followed the same principles as in Sect. 3.2. We used an MCMC sampler (emcee) with 112 walkers (four times the number of parameters) and 30 000 steps per walker. We finally removed the first half of each chain and combined all chains (1.68 × 10⁶ steps in total) to compute the final posterior distributions. Table B.3 shows the prior and posterior distributions for each parameter. The final models and confidence intervals for the planet transits are shown in Fig. 5, and the corresponding modeling of the radial velocity dataset is shown in Fig. 2 (full time series) and Fig. 3 (phase-folded dataset). In Fig. 7 we additionally show the joint modeling of the FWHM of the CCF sampling the stellar activity. The variations seen in this parameter mimic the variations seen in the radial velocity dataset, showing the correspondence between both measurements coupled by the activity of the star (see also Suárez Mascareño et al. 2020 for a similar analysis on Proxima b).

The inferred masses are compatible with the minimum masses obtained in the radial velocity analysis (Sect. 3.2). We find masses for the two planets that correspond to m_b = 6.38 ± 0.45 M_⊕ and m_c = 1.76 ± 0.17 M_⊕, slightly lower than the previously reported values. We also find slightly smaller radii for the two planets with the TESS data (R_b = 1.635 ± 0.46 R_⊕, R_c = 1.169 ± 0.038 R_⊕) than the values found by Ment et al. (2019) using Spitzer data ( R_⊕, R_c = 1.282 ± 0.024 R_⊕); they are different by about 1.5σ for LHS 1140 b and by 2σ for LHS 1140 c.

We also performed a three-planet modeling of the full radial velocity and TESS dataset to extract better constraints on the ephemeris of the third planet candidate. The resulting phase-folded radial velocity curve of the third planet is shown in Fig. 8. We can better constrain the time of conjunction for this planet to .

Fig. 13 Mass-radius diagram for the known exoplanets with the lowest masses with measured mass precisions better than30% (gray symbols). The new locations of the LHS 1140 planets are shown as filled blue (LHS 1140 b) and orange (LHS 1140 c) circles. Their previous locations using the values from Ment et al. (2019) are shown as open squares with dashed error bars and the same color code. The candidate planet at an orbital period of ~ 78 days is shown as the shaded magenta region. Its radius is not known, therefore its vertical location comes from the simple mass-radiusrelations from Chen & Kipping (2017). The bulk density lines corresponding to different compositions from Zeng et al. (2019) are shown as solid traces, and the dashed lines correspond to iso-densities of 1.33, 5.3, and 10 g cm−3 (from top tobottom).

7 Discussion

7.1 Mass-radius diagram

We estimated the physical and orbital properties of the two known planets with high precision thanks to the precise radial velocity measurements from ESPRESSO and the high-cadence and high-precision photometry from TESS. In Fig. 13 we show a mass-radius diagram including all planets with masses lower than 10 M_⊕ and radii smaller than2.2 R_⊕. The mass and radius precisions are better than 30%. We include the two confirmed planets in the LHS 1140 system with their derived properties estimated in this work and the previously known properties from Ment et al. (2019). The newly derived properties of both planets place them at the top of the Earth-like bulk density line, while previous measurement suggested a slightly lower density for these planets. The planet radius of LHS 1140 b places it close to the radius gap (Fulton et al. 2017), and the newly derived mass confirms the rocky nature of the planet. LHS 1140 c is one of the few planets with a derived mass lower than 2 M_⊕ in this diagram. The other two planets are GJ 357 b (Luque et al. 2019; Jenkins et al. 2019) and TRAPPIST-1 f (Gillon et al. 2016).

7.2 Transit of the planet candidate LHS 1140 d

The mass-radius diagram in Fig. 13 also shows the location of the third planet candidate (P_d~ 78 days) in the system, with an estimated mass of 3.9 M_⊕. The uncertainty on the planet radius is taken from the forecaster estimation. The measured mass and estimated radius do not allow us to discern between rocky or gaseous compositions. Based on its mass, the expected planet radius might be in the range 1.2–2.6 R_⊕. This impliesa potential transit depth of 3.5 ppt, which is suitable for ground-based instrumentation. Unfortunately, the precision of the time of conjunction from the radial velocity analysis is only about 3.5 days, which means that the ephemeris is too uncertain toplan ground-based observations. The next transits of this planet candidate will occur on Julian dates 2 459 167.8 ± 5.9 (2020-11-14 ± 5.9 days), 2 459 246.5 ± 6.4 (2021-02-01 ± 6.4 days), 2 459 325.2 ± 6.8 (2021-04-20 ± 6.8 days), 2 459 403.9 ± 7.3 (2021-07-08 ± 7.3 days), and 2 459 482.7 ± 7.6 (2021-09-25 ± 7.6 days). When a circular orbit coplanar with the other two known transiting planets is assumed, the impact parameter of this planet candidate would be b∕R_⋆ ~ 0.65, and it is therefore expected to transit its host star. LHS 1140 will be reobserved by TESS in its extended mission in Sector 30. However, despite the large uncertainty on the time of mid-transit, the TESS observations fall in between two transits of the planet and will therefore most likely not be detected by TESS (see Fig. 14). As a consequence, ground-based observations, and especially high-precision photometry, are needed to confirm this planet and further characterize its properties. Alternatively, LHS 1140 is observable with the Cheops space-based telescope (Broeg et al. 2013), although the large uncertainty on the ephemeris of the transit would require a large effort from the mission.

7.3 Internal structures of LHS 1140 b and LHS 1140 c

To investigate the interior of the confirmed planets, we used the internal structure model developed by Brugger et al. (2017) for the study of terrestrial planets. The input variables of the interior structure model are the total planetary mass, the core mass fraction (CMF), and the water mass fraction (WMF). In order to explore the parameter space, we performed a complete Bayesian analysis to obtain the probability density distributions of the parameters. The Bayesian analysis was carried out by implementing an MCMC algorithm, following the method proposed by Dorn et al. (2015). The initial values of the three input parameters were randomly drawn from their prior distributions, and we used a Gaussian distribution for the mass and two different uniform distributions for the CMF and the WMF. The uniform distribution for the CMF spans from 0 to 0.65, which is a constraint derived from the Fe/Si ratio in the protoSun (Lodders et al. 2009). Including this upper limit, we assumed that these planets have not undergone dramatic processes during or after their formation, such as mantle evaporation or giant impacts. The WMF ranges from 0 to 0.5, the upper value being derived from the composition of the Saturn moon Titan, which is one of the most hydrated Solar System bodies (Tobie et al. 2006). Along with these constraints, we took into account the masses and the radii obtained in this study for LHS1140 b and c to derive the posterior probability distributions of their CMF and WMF. We assumed a surface pressure and temperatureequal to those prevalent on Earth’s surface. This is valid for LHS 1140 b, but we note that LHS 1140 c is closer to the star, and the WMF estimated in this study therefore only represents an upper limit to its water content (which assumes a layer of liquid water, which is not viable in this case). Furthermore, for LHS 1140 c, we extended the interior structure model to include the modeling of water in supercritical phase, whose implementation is described in Mousis et al. (2020). We considered an atmosphere with a composition of 97% water and 3% carbon dioxide, and a surface pressure at its bottom of 300 bar. The atmospheric mass, thickness, albedo, and surface temperature are provided by a grid generated by the atmospheric model described in Marcq et al. (2017).

Figure 15 shows the results of this internal structure analysis by means of ternary diagrams for both planets. In the case of LHS1140 b, the CMF and WMF are estimated to be 0.49 ± 0.07 and 0.03 ± 0.07, respectively.In the case of LHS1140 c, assuming a liquid water layer, the CMF and WMF are estimated to be 0.45 ± 0.10 and 0.05 ± 0.07, respectively.This implies that the bulk composition of planet b is expected to be at least 10% richer in Fe than the Earth’s bulk composition. Planet c might be compatible with a CMF of 0.35 within its 1σ confidence interval, which is very close to the Earth’s 0.32 value. Furthermore, the WMF of LHS 1140 b implies a range of water content from no liquid water to 100 times more water than Earth. The posterior distribution peaks at WMF = 0.04, around 80 times the water content on Earth (0.0005–0.5% – Sotin et al. 2007). The 1% confidence level corresponds to WMF_1% = 0.007, still 1.5 times higher than on Earth. When the supercritical phase of water is considered for LHS 1140 c, the CMF is calculated as 0.59 ± 0.05, and the WMF is constrained as . If we were to assume a CMF similar to that of Earth for LHS 1140 c, its WMF would be lower than 0.05%, which means that the inner planetis very poor in water (see Fig. 16).

Fig. 14 Expected posterior distribution of the mid-transit time of planet candidate LHS 1140 d (magenta histograms) and the TESS coverage of Sectors 3 and 30 on the field including the host star (shaded gray region). The median of each posterior distribution is marked by a vertical dotted line, and the corresponding calendar date in Universal Time Coordinate system is annotated.

Fig. 15 Sampled 2D marginal posterior distribution for the CMF and WMF of LHS1140 b. The color code displays the probability density function (PDF). The solid and dashed black lines represent the isoradius curves for the central value of the radius and its 1σ confidence interval limits, respectively. The gray lines delimit the areas we excluded from our sampling (see text).

Fig. 16 Mass-radius diagram for different compositions considering water in supercritical phase under the surface temperatureand pressure conditions of LHS 1140 c. The black dot and error bars indicate the position of LHS 1140 c in the mass-radius diagram.

8 Conclusions

We have revisited the planetary system LHS 1140 with its two known Earth-like transiting planets by analyzing an intensive campaign taken with the new ultra-stable high-resolution spectrograph ESPRESSO and observations from Sector 3 of the TESS mission. We have searched for additional signals in the radial velocity data including HARPS and ESPRESSO measurements. The results show additional evidence for the signal of a third planet candidate at ~ 78 days. Our analysis of the joint radial velocity dataset, including a modeling of the stellar activity using the FWHM of the CCF as an activity indicator, shows positive evidence that this signal is of planetary nature (i.e., not caused by activity). We find a ~ 4σ significance in the semiamplitude of the radial velocity of this signal. This is supported by our statistical analysis with kima and the l1-periodogram (see Sect. 3.2). The Bayesian evidence of the two-planet model against the three-planet model, however, is still not enough to claim a confirmation, but its evolution when the dataset was increased points toward the planetary scenario. This evolution shows that ~ 50 more ESPRESSO measurements are required to unambiguously confirm the planet signal based on radial velocity data alone. The corresponding planet mass of the candidate planet is 3.9 ± 1.1 M_⊕. Unfortunately, no transit signal is present in Sector 3 of the TESS data, as expected from the ephemeris of our joint LC+RV analysis. This ephemeris also suggest that the extended TESS mission will not catch the planet transit when the mission revisits this field in Sector 30 between 22 September 2020 and 21 October 2020. If the system is coplanar (as the two known planet inclinations strongly suggest), the impact parameter of this planet candidate would certainly mean that it transits the host star. The expected total transit duration (assuming circular orbit) is about three hours. Further ground-based photometry to search for this transit would therefore be extremely useful and efficient to confirm this signal. We caution about the relatively large uncertainties of the ephemeris, however, which would imply a dedicated campaign of about five consecutive nights on the same target.

We also explored the possibility of co-orbital planets in the two known transiting planets through a dedicated study of the radial velocity data and light curve. For LHS 1140 b, the radial velocity data allowed us to constrain the mass of any tadpole and horseshoe co-orbitals up to 2.1 M_⊕ in L₄ and 1.8 M_⊕ in L₅ (95% confidence levels). Similarly, the TESS data discard nearly coplanar co-orbitals around the Lagrangian point regions larger than ~ 3 R_⊕. On the other hand, for LHS 1140 c, the radial velocity analysis unveils a 2σ signal in its L₄ region that would correspond to a sub-Earth co-orbital planet with a 0.26 ± 0.18 M_⊕. The Bayesian evidence of this model still does not support this scenario against the null hypothesis, however, and several hundred additional measurements with ESPRESSO-like precision are required to confirm the signal with radial velocity data alone. However, the analysis of the TESS light curve indicates a shallow dimming very close to this Lagrangian point with the same duration as the transit of LHS 1140 c (as expected in the coplanar case). If real, this dimming would correspond to a 0.44± 0.08 R_⊕ sub-Earth size body. Unfortunately, the significance of the signal is not strong enough to ensure that this dimming is not caused by systematics of the TESS light curve, and only additional data will shed more light on this co-orbital candidate.

Finally, we performed a joint analysis including the radial velocity and light curve datasets assuming the two-planet scenario to derive precise properties of the known transiting planets. We used the largest evidence model from the radial velocity-only analysis, corresponding to the two-planet case with circular orbits for both planets (this is favored by the Bayesian evidence with a Bayes factor larger than compared to the eccentric case). The results of this analysis show slightly less massive and smaller planets than previously reported values. For LHS 1140 b, we find a planet mass and radius of m_b = 6.48 ± 0.46 M_⊕ and R_c = 1.641 ± 0.048 R_⊕, resulting in a bulk density of g cm⁻³. For LHS 1140 c, we obtain m_b = 1.77 ± 0.17 M_⊕ and R_c = 1.185 ± 0.044 R_⊕, giving g cm⁻³. This analysis provides unprecedented mass and density precisions for such small rocky planets (9 and 7% for the outer and inner components, respectively). With these new estimates, both planets lie exactly on the Earth-like bulk density line in the mass-radius diagram (see Fig. 13). Taking advantage of the precise mass and radius measurements, we performed an internal structure analysis considering liquid water conditions and derived a core mass fraction for LHS 1140 b and LHS 1140 c of 49 ± 7 and 45 ± 10%, respectively. It is remarkable that these values are close to the Earth core mass fraction (32%). Additionally, we found a remarkable water content on LHS 1140 b, with the water mass fraction peaking at 4% (80 times more than on Earth). The posterior distribution for this water mass fraction is truncated at zero, however, thus not allowing a definitive confirmation of a water layer on this planet, but instead providing strong indications that LHS 1140 b is a large true water world. This water mass fraction would imply a huge ocean with a depth of 779 ± 650 km according to our internal structure model. On the other hand, if water in supercritical phase is considered for LHS 1140 c, its maximum WMF is estimated to be 0.06%, which shows that planet c is likely to be dry or very poor in water. Thus we can conclude that the system of LHS 1140 presents a gradient of water mass fraction with irradiation, with an inner dry planet and a wet outer planet. LHS 1140 b and LHS 1140 c might have acquired their water content because they formed beyond the snow line and then migrated inward. In the case of planet c, a water-dominated atmosphere could have undergone atmospheric escape due to the high irradiation it received from its host star, which caused it to loose most of its water content. The analysis performed in Sect. 3.3, however, rejects the presence of additional planets more massive than 1 M_⊕ in between LHS 1140 c and LHS 1140 b that could form a resonant chain that could explain the inward migration of these planets from orbits beyond the snow line (Delisle 2017). Although several mechanisms can explain the disruption of these chains with hot super-Earths (Cossou et al. 2014), additional theoretical studies should focus on reconciling the interior structure of these planets (requiring that they formed beyond the snow line) with the fact that the planets in the system are currently not in resonance (usual outcome of inward migration in multiplanet systems, see, e.g., Delisle 2017), and the system hosts no additional components in orbits between the two known planets that can create a resonant chain (see Sect. 3.3).

With our unprecedentedly precise mass and density measurements for the two transiting rocky worlds, LHS 1140 becomes one of the prime systems for astrobiological studies. The habitable-zone planet LHS 1140 b is now known to contain a significant amount of water with a relatively high probability. Previous works have studied the actual state of water layers on the surface of this planet by analyzing the possible eyeball, lobster, or snowball scenarios (Yang et al. 2020). Future atmospheric studies with the James Webb Space Telescope and extremely large telescopes will be able to distinguish among these scenarios. These observations will greatly benefit from the precise measurements presented here. Added to this, the possible presence of a third rocky planet in the system that may be transiting the host star at periods beyond the habitable zone makes this system a key target for understanding atmospheric properties of rocky worlds at different stellar irradiations.

Acknowledgements

This research has been funded by the Spanish State Research Agency (AEI) Projects No.ESP2017-87676-C5-1-R and No. MDM-2017-0737 Unidad de Excelencia “María de Maeztu” – Centro de Astrobiología (INTA-CSIC). This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019; UIDB/04434/2020; UIDP/04434/2020; PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113; PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953. N.A.-D. acknowledges the support of FONDECYT project 3180063. X.B. and J.M.A. acknowledge funding from the European Research Council under the ERC Grant Agreement n. 337591-ExTrA. A.C. acknowledges support by CFisUC projects (UIDB/04564/2020 and UIDP/04564/2020), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), funded by COMPETE 2020 and FCT, Portugal. J.P.F. is supported in the form of a work contract funded by national funds through FCT with reference DL57/2016/CP1364/CT0005.

Appendix A: Additional figure

Fig. A.1 Bayesian evidence of different models and datasets (from top to bottom: HARPS-only, ESPRESSO-only, and HARPS+ESPRESSO). Each panel contains 15 different models that include a different number of planets: no planets (GP-only, yellow), one-planet models (pink), two-planet models (red), and three-planet model (green). Each has different orbital configuration assumptions (planets with either circular or eccentric orbits). As an example, model 3p1c2c corresponds to a three-planet model in which planets 1 (LHS 1140 b) and 2 (LHS 1140 c) have assumed circular orbits, but planet 3 (LHS 1140 d) has free eccentricity. The bars corresponding to each model include this information at the bottom of the bar. At the top of each bar we show the Bayes factor of each model compared to the stronggest evidence model for each dataset. The highest bar corresponds to the strongest evidence model.

Appendix B: Additional tables

References

All Tables

All Figures

Fig. 1 Periodogram of the radial velocity (upper panel) for each dataset individually (ESPRESSO in light green and HARPS in red) and for the joint dataset of both instruments assuming a 26 m s−1 offset (see Sect. 6). The activity indicators for the individual datasets are also shown in the lower panels. The periods of the three planets are marked as dotted vertical lines, and the rotation period (Prot = 131 days) and its first harmonic (Prot∕2 = 65 days) are marked as vertical dashed lines.

In the text

Fig. 2 Top: radial velocity of LHS 1140 from the HARPS (red) and ESPRESSO (open for ESPRESSOpre and filled for ESPRESSOpost in green) datasets. The black line shows the median radial velocity model from the joint photometric and radial velocity analysis including the two known Keplerian signals and the GP model (see Sect. 6). The gray shaded regions correspond to the 68.7% (dark gray) and 95% (light gray) confidence intervals of the model. The median GP model is shown as a dashed blue line. Middle: radial velocity dataset after removing the median GP model. The Keplerian model is shown as a solid black line. Bottom: radial velocity residuals of the full model.

In the text

Fig. 3 Phase-folded radial velocity signal of LHS 1140 b (left panels) and LHS 1140 c (right panels). For each planet, upper panel: radial velocity data (color-coded as in the legend), binned radial velocities (open black symbols) with a bin size corresponding to 15% of the phase, the median model for each planet obtained from the joint fit analysis (see Sect. 6) as the solid black line, and the 68.7 and 95% confidence intervals as shaded dark and light gray regions. Bottom panels: residuals of the median model.

In the text

	Fig. 4 TPF of LHS 1140 from the TESS observations in Sector 3 (composed with tpfplotter, Aller et al. 2020). The SPOC pipeline aperture is overplotted with shaded red squares, and the Gaia DR2 catalog is also overlaid with symbol sizes proportional to the magnitude contrast with the target, marked with a white cross.
In the text

Fig. 5 TESS light curve extracted from the SPOC pipeline and detrended using the PDCSAP flux. Upper panel: complete light curve displaying the transits of the two known transiting planets in the system, LHS 1140 b (red, one transit) and LHS 1140 c (light green, five transits). Lower panel: phase-folded light curve centered on the planet phase for each of the two transiting planets, LHS 1140 b (bottom left) and LHS 1140 c (bottom right). The colored symbols correspond to bins of one fifth of the transit duration (i.e., 30 min for LHS 1140 b and 11 min for LHS 1140c).

In the text

	Fig. 6 ℓ1 -periodogram of the LHS 1140 HARPS and ESPRESSO data with the noise model that has the best cross-validation score.
In the text

Fig. 7 Dataset of the FWHM of the cross-correlation function from the HARPS and ESPRESSO data. The systemic offset for each instrument as retrieved from the two-planet model has been removed. The median GP model is shown as the solid black line, and the 1σ and 3σ confidence intervals are shown as shaded dark and light gray regions, respectively. Bottom panel: residuals after removing the GP model.

In the text

Fig. 8 Phase-folded radial velocity signal of the candidate planet LHS 1140 d in the three-planet model. Upper panel: radial velocity data (color-coded as in the legend), binned radial velocities (open black symbols) with a bin size corresponding to 15% of the phase, the median model from the joint three-planet analysis (see Sect. 6) as the solid black line, and the 68.7 and 95% confidence intervals as shaded dark and light gray regions. Bottom panels: residuals of the median model.

In the text

	Fig. 9 Difference in the log-evidence of the two-planet and three-planet models with circular orbits for a cumulative number of data points. Top panel: Bayes factor against the time span of the data set, and bottom panel: against the number of data points.
In the text

Fig. 10 Detection limits of the radial velocity dataset for periods shorter than one year. The color code indicates the significance of the radial velocity semiamplitude parameter. The 5σ contour is shown as a solid line, and the limit corresponding to the scatter of the data is shown as a dashed line. The locations of the two known planets are shown as solid circles, and the location of the third planet candidate is shown as an open circle.

In the text

	Fig. 11 TLS periodogram of the TESS dataset after masking the light curve from the transit times corresponding to the two known planets LHS 1140 b and LHS 1140 c. The red line indicates the strongest power peak, and the blue line indicates the period of LHS 1140 c.
In the text

Fig. 12 TESS phase-folded light curve around the Lagrangian points L5 (top panel) and L4 (middle panel) of LHS 1140 c (see Sect. 5.2). The expected location of the transit is marked by vertical dashed lines. The median transit model inferred from the analysis of the L4 dimming is shown as a solid red line (middle panel) and a dashed line in the upper panel to guide the eye on the photometric scatter of the L5 region. The 1σ confidence interval is shown in the middle panel as a shaded red region. Bottom panel: residuals of the model for the L4 region. All panels show bin data points corresponding to 6.6 min (light blue, ten data points inside the transit duration) and 13.2 min (black, five data points inside the transit duration).

In the text

Fig. 13 Mass-radius diagram for the known exoplanets with the lowest masses with measured mass precisions better than30% (gray symbols). The new locations of the LHS 1140 planets are shown as filled blue (LHS 1140 b) and orange (LHS 1140 c) circles. Their previous locations using the values from Ment et al. (2019) are shown as open squares with dashed error bars and the same color code. The candidate planet at an orbital period of ~ 78 days is shown as the shaded magenta region. Its radius is not known, therefore its vertical location comes from the simple mass-radiusrelations from Chen & Kipping (2017). The bulk density lines corresponding to different compositions from Zeng et al. (2019) are shown as solid traces, and the dashed lines correspond to iso-densities of 1.33, 5.3, and 10 g cm−3 (from top tobottom).

In the text

	Fig. 14 Expected posterior distribution of the mid-transit time of planet candidate LHS 1140 d (magenta histograms) and the TESS coverage of Sectors 3 and 30 on the field including the host star (shaded gray region). The median of each posterior distribution is marked by a vertical dotted line, and the corresponding calendar date in Universal Time Coordinate system is annotated.
In the text

	Fig. 15 Sampled 2D marginal posterior distribution for the CMF and WMF of LHS1140 b. The color code displays the probability density function (PDF). The solid and dashed black lines represent the isoradius curves for the central value of the radius and its 1σ confidence interval limits, respectively. The gray lines delimit the areas we excluded from our sampling (see text).
In the text

	Fig. 16 Mass-radius diagram for different compositions considering water in supercritical phase under the surface temperatureand pressure conditions of LHS 1140 c. The black dot and error bars indicate the position of LHS 1140 c in the mass-radius diagram.
In the text

Fig. A.1 Bayesian evidence of different models and datasets (from top to bottom: HARPS-only, ESPRESSO-only, and HARPS+ESPRESSO). Each panel contains 15 different models that include a different number of planets: no planets (GP-only, yellow), one-planet models (pink), two-planet models (red), and three-planet model (green). Each has different orbital configuration assumptions (planets with either circular or eccentric orbits). As an example, model 3p1c2c corresponds to a three-planet model in which planets 1 (LHS 1140 b) and 2 (LHS 1140 c) have assumed circular orbits, but planet 3 (LHS 1140 d) has free eccentricity. The bars corresponding to each model include this information at the bottom of the bar. At the top of each bar we show the Bayes factor of each model compared to the stronggest evidence model for each dataset. The highest bar corresponds to the strongest evidence model.

In the text