© The Authors 2024

1 Introduction

From the pioneering CoRoT space mission (Fridlund 2008) that discovered the first rocky super-Earth (Léger et al. 2009), through the Kepler space telescope (Borucki et al. 2010) that truly showed the diversity of exoplanets, to the now active all-sky transit survey TESS (Ricker et al. 2015) and most recently the very precise CHEOPS (CHaracterising ExOPlanet Satellite; Benz et al. 2021), the advent of space-based transit photometry has been a game changer.

The more than 5000 planets detected to date seem to be just the tip of the iceberg. While the discovery of new systems is continuing at a fast pace, we are simultaneously faced with the formidable task of obtaining the first grains of comprehension of the incredible diversity of both the planets we discover and the host stars they orbit. CHEOPS is the first small space mission in the European Space Agency (ESA) Science Program. It was launched on December 18 2019 and is a partnership between ESA and Switzerland with important contributions from other ESA member states. CHEOPS is the first mission dedicated to directly observing transits of already-identified individual exoplanets. While it has achieved a photometric precision in the transit light curves (LCs) only surpassed by the Kepler mission, and has thus enabled significant improvements in planetary parameters over previous exoplanetary missions, CHEOPS operates in a manner completely different from those. Earlier spacecraft were all carrying out survey missions, while CHEOPS has as its objective the detailed study of individual planets (e.g. Lacedelli et al. 2022; Wilson et al. 2022). It uses ultra-high precision photometry to observe the predicted transits of single, specific targets orbiting bright stars.

In recent years, low-mass M dwarf stars have been gaining traction among the exoplanet community in the search for new worlds. This is understandable for a number of reasons. Transit and radial velocity (RV) surveys are much more likely to detect potentially habitable worlds given the proximity of the habitable zone to an M-dwarf host star. Because of their small radii and low masses, compared to Sun-like stars, relatively large transit depths and RV amplitudes are induced by planets transiting M dwarfs. The probability of uncovering planetary habitability is further increased by the fact that red dwarfs have been estimated to be by far the most numerous objects in the Galaxy, comprising 75% of the stars in our Galaxy.

While impressive, the above advantages should be taken with a pinch of salt since the faintness of M dwarfs can also make them difficult to study in the context of exoplanets. Another important factor to consider in the case of RV surveys is the stellar activity-induced signal, which would be higher for the typically more magnetically active cool stars (e.g. Reiners et al. 2010; Andersen & Korhonen 2015). In some cases the activity displays stochastic behaviour and cannot be described by traditional methods, such as for example sinusoid-fitting. Such situations can further be exacerbated by gaps in the data caused by interruptions to the observing run. Such is the case with TOI-776 which is an early M dwarf hosting two small transiting planets in the super-Earth – mini-Neptune transition regime (Luque et al. 2021, hereafter L21). The mass range of a few to ~ 10 M_⊕ that these planets populate is an interesting one since such planets are not only the most numerous, but are also not represented in our Solar System, despite the apparent diversity of its members. The discovery, as first reported by L21, relied primarily on TESS photometry providing three transits of the inner planet and two transits of the outer one, with additional help from four ground-based transits for both planets, three of which were partial. This led to the precision of the radii being about 7%. The mass determinations of the two planets were based on 29 RV measurements that led to a precision of approximately 30%.

Given the estimated location of these planets with regards to the radius valley for M dwarfs (Van Eylen et al. 2021), the potential they have for gaining further insight into planet formation and evolution mechanisms is substantial. Moreover, TOI-776b and TOI-776c are particularly suitable for atmospheric characterisation studies, as highlighted by L21. Added to the fact that these two objects are already selected for observations in cycle 1 (program ID 2512, PI: Batalha) of JWST (Gialluca et al. 2021), updated planetary and orbital parameters will certainly also be very useful.

With the above justification in hand, TOI-776 was given a high priority as a target for observations with the CHEOPS space mission and the KESPRINT program utilizing the HARPS high-precision spectrograph. Our studies using CHEOPS are intended to determine facts about planets and planetary systems that could be helpful in narrowing down the correct formation models. This would further the understanding of the causes of the large diversity that exists among types of planets, as well as the different distribution of individual planets within their systems. Furthermore, so far it appears that the properties of systems depend to a large degree on the properties of the host star, but the role of the type of host star in these issues in general needs to be investigated in greater detail (Perryman 2018; Deeg & Belmonte 2018; Fridlund et al. 2020).

In order to make a precise characterisation of the sizes of the planets, our observations of TOI-776 with CHEOPS are aimed at improving the precision in planetary radii by at least a factor of two with respect to L21. We also wanted to expand on the RV work of L21 with the aim of improving the precision in the mass determinations and, thus, in the mean densities. The latter goal proved elusive given both the apparently changed stellar behaviour between the two radial velocity campaigns and the questionable approach based on sinusoid-fitting adopted by L21 in their RV analysis. To determine the planetary masses we instead used a novel approach, which is based on applying a skew normal (SN) fit onto the cross-correlation functions (CCFs) retrieved from the HARPS spectra (Simola et al. 2019). This allowed us to extract the RV measurements and stellar activity indicators, without introducing any ad hoc hypotheses to model the stellar activity during the detrending phase.

In this paper, we present the new observations (Sect. 2) and derive the stellar parameters as well as determine the level of activity in Sect. 3. Our analysis, using both the multi-Gaussian process (GP) method as well as our reference SN-fit method to extract the RV data, is described in Sect. 4. In Sect. 5 we derive the planetary parameters, and in Sects. 6 and 7 we discuss the increased stellar activity, its impact on the achievable determination of planetary parameters, and set out our conclusions about the planets.

2 Observations and data

We collected photometric and spectroscopic observations with the aim of performing a joint fit of the data (Sect. 4) and specifically retrieving new radii and masses of TOI-776 b and TOI-776 c.

2.1 TESS photometry

The NASA TESS space mission (Ricker et al. 2015) has been launched with the objective of discovering transiting exoplanets. Its wide field is intended to make possible the identification and immediate characterisation of exoplanets orbiting brighter stars than its predecessors. TOI-776 was first observed by TESS (two-minute cadence mode) in March-April 2019 in Sector 10 on camera 2, CCD 4, when two exoplanets now identified as b and c were initially flagged by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016). The target was further observed by TESS in its two-minute cadence mode in Sector 37 between 2 and 28 April 2021, once again on camera 2, CCD 4. We note that this latter data was not included in the analysis of L21.

For the transit analysis, we used the TOI-776 TESS light curves, as extracted by the SPOC pipelines (Twicken et al. 2010; Morris et al. 2017), with instrumental systematics and dilution corrected for via the presearch data conditioning simple aperture photometry (PDCSAP) algorithm (Smith et al. 2012; Stumpe et al. 2012). We performed a five median-absolute-deviation (MAD) clipping to reject flux outliers. From each observation run containing ~ 1 month of data, we extracted those temporal windows centred around each transit of TOI-776 b and c, keeping also ~4 h of out-of-transit data both before and after the transit event for detrending purposes. We ended up with four LCs from Sector 10 (where one contains two transit events as they are very close in time) and five LCs from Sector 37. Besides the time and the flux with its errors, each of the nine TESS LCs also contains the PDCSAP data product vectors to test whether any further detrending is needed. The raw TESS LCs are displayed in Fig. A.1. The phase-folded light curves of the TESS data are displayed in Fig. 1.

2.2 CHEOPS photometry

The CHEOPS spacecraft is in a low Earth orbit (LEO) with a period of 98.77 min. It is equipped with a 32 cm Ritchey-Chretien telescope feeding a back-illuminated CCD photometer operating between the 0.33 μm and 1.1 μm wavelength range. The raw LCs are shown in Fig. A.2. More details about CHEOPS are available in Benz et al. (2021).

We used CHEOPS to observe five visits of TOI-776 between March 14 and May 4, 2021. These observations yielded a total of 52.8 h of LC data. We detected the transits of TOI-776b in two visits and of TOI-776c in three visits (see Table 1 and Fig. A.2). The observations were carried out and then reduced using the CHEOPS data reduction pipeline (DRP, v.13; Hoyer et al. 2020). In brief, after downlinking the data as CCD windows of ~200″ diameter, the DRP subtracts the bias and corrects for non-linearity and dark current, while taking flat field variation into account. The DRP corrects also for the sky-background, cosmic ray impacts, and also smearing trails of stars close to the line of sight. The DRP performs automatic aperture photometry on the processed CHEOPS images using different circular masks centred around the target. In our case we carried out photometry using the DEFAULT aperture (i.e. 25″). A stable photometry was achieved by letting the mask follow the movements of TOI-776 as the spacecraft jittered and rolled around the optical axis.

The DRP also creates a set of vectors, allowing the user to maximise performance during the following stages of reduction. These vectors consist (among others) of the orbital roll angle (roll), the x and y positions on the CCD of the centre of the point spread function (PSF), the estimated background light (bg, e.g. due to zodiacal light), the level and position of the smear factor (smear), and the degree of contamination by background stars (conta). This data can be found in the CHEOPS archive at the Data & Analysis Center for Exoplanets (DACE)¹. The phase-folded light curves of the CHEOPS data are displayed in Fig. 1.

2.3 Ground-based photometry

In order to treat the available data in an optimal way, we also considered the ground-based light curve photometry obtained by L21 from the MEarth-South, and three nodes of the LCO, namely LCO-CTIO, LCO-SSO, and LCO-SAAO. We refer to L21 for details. The lightcurve data from these observations were integrated together with the old and new TESS observations as well as the CHEOPS observations, and thus contributed to the determination of the rotation period of the host star. Stellar designation, literature photometry as well as parallaxes are given in Table 2.

2.4 Spectroscopic observations

The two planetary candidates identified by TESS (Sect. 2.1) were confirmed through spectroscopic follow-up observations by the KESPRINT consortium² (see L21), also using the ground-based transit photometry quoted above (Sect. 2.3). In order to increase the precision in the mass determinations, we collected 35 additional RV data points with the HARPS spectrograph (R =115000) mounted at the ESO 3.6m telescope (La Silla observatory, Chile: LP 106.21TJ.001; P.I. D. Gandolfi). We used the second fibre of the spectrograph to monitor the sky background and set the exposure time to 1800–2520 s depending on the sky conditions and scheduling constraints, leading to a median signal-to-noise ratio of ~42 per pixel at 550 nm. We reduced the data using the dedicated data reduction software (DRS; Pepe et al. 2002; Lovis & Pepe 2007) and computed the CCFs from each Echelle spectrum using a numerical M2 mask (Baranne et al. 1996). With 35 new HARPS spectra, the resulting data set contains 64 HARPS spectra in total, which is more than twice the number of Doppler measurements presented in L21. The RV observations now cover the time³ between BJD_TDB ≈ 2 458 884.8 and BJD_TDB ≈ 2 459 430.5, but with a significant gap between BJD_TDB ≈ 2451931 and BJD_TDB ≈ 2452 358 due to the COVID-19 pandemic (see Figs. 2 and 3, and Sect. 4.2).

3 Stellar parameters

3.1 Spectroscopic, photometric, and isochronal parameters

The actual physical parameters of single red dwarf stars are notoriously difficult to determine, for example Schweitzer et al. (2019). Therefore, we used several different methods based on both stellar photometry as well as analysis of the observed high-resolution spectrum in order to retrieve the stellar basic parameters, T_eff, log g, [Fe/H], υ sin i_*, V_mic, M_*, R_* and an estimate of the stellar age.

For spectral data, we achieved a high signal-to-noise spectrum (S/N = 350 per pixel at 550 nm), by co-adding the individual 64 HARPS radial velocity spectra. This spectrum was then normalised and used as an input for our stellar analysis. The photometric data was obtained from the NASA Exoplanet Archive⁴ and references there.

Firstly, and also taking the stellar velocities (V_mic, V_mac, and υ sin i_*) into account is the IDL code Spectroscopy Made Easy (SME; Valenti & Piskunov 1996; Piskunov & Valenti 2017), which synthesises models of individual absorption lines in the observed spectrum. We followed Fridlund et al. (2020) and references therein, and found υ sin i_* to be 2.2 ± 1.0 km s⁻¹. Using SME to fit several hundred TiO lines with T_eff as the only free parameter, we then found T_eff to be 3725 ± 50 K. Fixing this T_eff, we then found [Fe/H]−0.21 ± 0.08.

Secondly, we used Specmatch-emp (Yee et al. 2017), a code⁵ that compares the observed spectrum with a library of over 400 spectra of stars, of all types and with well-determined physical parameters. A minimising and interpolation calculation provides values of T_eff, log g, and [Fe/H], as well as estimates of M_*, R_*, and age (see Fridlund et al. (2020) for more detail on this method).

Thirdly, we applied the new ODUSSEAS code. This tool (Antoniadis-Karnavas et al. 2020)⁶ is based on the measurement of the pseudo-equivalent widths from more than 4000 stellar absorption lines and using the machine-learning Python package scikit-learn to derive the T_eff and [Fe/H] parameters accurately and with high precision. This code uses a library of HARPS spectra of M stars, with interferometric calibrations. In the library spectra pseudo-equivalent widths were measured for hundreds of lines, and then the code was trained using these reference parameters. This gave us the two parameters, T_eff and metallicity, with statistical precision errors of 30 K for T_eff and 0.04 dex for [Fe/H].

We also fitted catalogue photometry, using two methods, astroARIADNE and a version of IRFM, to calculate the spectral energy distribution (SED) of TOI-776.

The Python code astroARIADNE⁷ (Vines & Jenkins 2022) uses a Bayesian model averaging a number of stellar atmospheric model grids to obtain estimates of the spectroscopic parameters and then derive M_* and R_*, following Persson et al. (2022).

We then checked this result with a modified infrared flux method (IRFM; Blackwell & Shallis 1977) using a Markov chain Monte Carlo (MCMC) approach to determine the stellar angular diameter and the effective temperature (Schanche et al. 2020). Following Wilson et al. (2022), we found essentially the same (within 1σ) results as with astroARIADNE, that is, R_* = 0.547 ± 0.017 R_⊙.

Our results with these different methods are presented in Table 3. As can be seen, they are all very much in agreement and we selected those with the most realistic errors as our adopted values.

Using the adopted values (T_eff, [Fe/H], R_*) as a reference input set, we derived a robust estimate of the stellar mass, M_* by using two different stellar evolutionary codes, namely PARSEC⁸ v1.2S (Marigo et al. 2017) and CLES⁹ (Scuflaire et al. 2008). In particular, we fitted the input set into pre-computed PARSEC grids of isochrones and tracks through the isochrone placement algorithm developed by Bonfanti et al. (2015, 2016). Here we used the SME-derived υ sin i_* value in order to improve the convergence of the interpolating routine as discussed in Bonfanti et al. (2016), and we obtained a first estimate of the stellar mass. Furthermore, starting from our reference input set, we employed the CLES code to generate the best-fit evolutionary track according to the Levenberg-Marquadt minimisation scheme (Salmon et al. 2021) and we inferred a second estimate of the stellar mass. Finally, we checked the consistency of the two outcomes through a χ²-based criterion and we merged them together, obtaining $M_{\star }={0.542}_{-0.039}^{+0.040}{M}_{\odot }$ (see Bonfanti et al. 2021a, for further details). Following a similar procedure we arrived at a stellar age, albeit with large uncertainties, of ${6.1}_{-5.5}^{+7.0}$ Gyr.

All the relevant stellar parameters as well as the final adopted ones, can be found in Table 3.

Fig. 1 Phase-folded LCs showing the transits of TOI-776b (first column) and TOI-776c (second column). Top panels: detrended CHEOPS LCs; the best-fit transit models are shown as red lines together with their corresponding residuals. Bottom panels: same as top, but for the TESS LCs. Data are shown in the nominal cadence modes (blue dots) and also binned to 10 min (black markers) with their associated error bars. We note that the transit durations and depths as measured by TESS and CHEOPS are identical within uncertainties. The transit depths found by CHEOPS are 909 ppm for planet b and 1177 ppm for planet c.

Fig. 2 Time series of the activity indicators. The complete RV data set has been used. The thick red lines mark the best-fitting line fit. We note the change in the activity level between the two seasons significantly detected in the contrast (first panel), the DLW (fourth panel), and the Ha index (seventh panel). The time is given as BJDTDB = Barycentric Julian Dates in Barycentric Dynamical Time − 2 450 000.

3.2 Stellar activity and rotation period

Red dwarf stars are generally found to be considerably more active than solar-like stars and TOI-776 is no exception. As mentioned in Sect. 2.4, there is a long gap in our RV sequence. It was found that TOI-776 was significantly more active during the second observing period.

We used the HARPS DRS (Sect. 2.4) to extract three profile activity diagnostics of the CCF, namely, the contrast, the full width at half maximum (FWHM), and the bisector inverse slope (BIS). We then followed the approach used in, for example, Fridlund et al. (2020) and extracted two independent sets of relative RV measurements from the HARPS spectra using the TERRA (Anglada-Escudé & Butler 2012) and SERVAL (Zechmeister et al. 2018)¹⁰ software packages.

While being similar in principle, the SERVAL code provides a different set of activity-related outputs. Both codes use algorithms that first create a template from the total assemblage of HARPS spectra and then compare each individual spectra with the template. As noted by, for example, L21, these methods can extract more precise RV measurements compared with the DRS package hitherto used with HARPS and similar spectrographic observations. This is true especially for M dwarfs and cool stars where the atomic and molecular spectral line lists are incomplete, with their spectra displaying strong line blending and continuum suppression. The RVs have a median internal uncertainty of 1.5 m s⁻¹ (resp. 1.5 m s⁻¹) and a root mean square of 5.2 m s⁻¹ (resp. 3.5 m s⁻¹) around the mean value for the SERVAL (resp. TERRA) extractions.

We specifically used SERVAL and TERRA to compute additional stellar activity indicators, namely, the chromatic index (CRX), the differential line width (DLW), the Ca II H & K lines Mount-Wilson S-index, and the Ha and Na D line indices.

We finally extracted a fourth set of relative RVs via a SN fit (Simola et al. 2019; see Sect. 4.2). The RV median internal uncertainty is 2.3 m s⁻¹, with a root mean square of 2.5 m s⁻¹.

We report the radial velocities and activity indicators extracted from the HARPS spectra in the data at the CDS and in Sect. 4.2. In particular, the two different values of extracted RVs (SERVAL and TERRA) with their associated uncertainties, the S-index, and the Ha and Na D line indices, together with the CRX and DLW, are all reported in a table at the CDS. The RV based on the SN fit along with the inherent activity indicators are also reported in a table at the CDS.

Figure 2 displays the time series of the activity indicators extracted with the DRS, SERVAL, and TERRA codes. The CCF contrast displays a significant decrease between the two observing periods, with a Pearson correlation coefficient of −0.87 (Fig. 2, first panel). Increasing trends are detected in the DLW and in the Ηα-index, with a Pearson correlation coefficient of 0.83 and 0.50, respectively (Fig. 2, fourth and seventh panel). As described in Kürster et al. (2003) and Zechmeister et al. (2018), the Ηα-index is defined as the ratio between the flux in the line core and the flux in two nearby reference spectral regions around the line. The increasing trend detected in the Ηα-index implies that the emission component has become stronger, increasing the filling factor of the line, and, consequently, the index of the absorption line. This translates into an increased level of stellar activity in the second observing period.

There are two aspects of this changing activity that need to be addressed. First, whether we can use the complete RV data set in order to also determine the stellar rotation period (P_rot) and, second, whether we can use the complete RV data set for the mass determination of both planets (as compared to using only the RV data of L21) without introducing more noise (and thus larger errors).

First, we computed the generalised Lomb–Scargle (GLS) periodogram (Zechmeister & Kürster 2009) of the RV measurements, the Hα, the DLW, the CRX, and the NaD1 indexes from SERVAL (see Fig. 4. After identifying the known signals from the two planets, “b” and “c” (with the periods obtained from the transits), from the RV and from the different activity indexes, we find a signal with a period of 21–22 days (the peak at ${21.129}_{-0.058}^{+0.062}$ days or 0.0473 day⁻¹ is shown by the yellow marker in Fig. 4). We analysed the data for each parameter in three ways: Epoch 1, which refers to the data first presented in L21; Epoch 2, where we used only the new RV spectra taken after BJD_TDB = 2452 358; and finally we analysed these data together as a set. We find that the DLW give discrepant results until we remove the offset seen in Fig. 2. After this removal we get very close values for the peak locations from all indexes, indicating that we have managed to remove a significant amount of the difference in the level of activity. It is therefore tempting to identify the 21.1d day signal with the rotation period of the star. This assumption is strengthened by our modelling efforts described in Sect. 4.1. A summary of the results can be seen in Fig. 4 where we also report the false alarm probability (FAP). We discuss further the second aspect, the impact of activity, on the RV analysis, below in Sect. 4.2.

Fig. 3 Whole RV time series: the break on the temporal axis avoids the big gap in the data collection and improves the visualisation. Top panel: undetrended RV data with the model of the full Keplerian signal + stellar activity superimposed in blue. Bottom panel: detrended RV data with the model of the full Keplerian signal superimposed in red.

4 Analysis

4.1 Preliminary multi-GP analysis

We began our investigation of the system by modelling the data using the code pyaneti¹¹ (Barragán et al. 2019, 2022). As discussed in Sect. 3.2, the activity of TOI-776 has changed significantly since the first RV observations published by L21. As a result, we find that modelling all RVs as a single data set and using a wide uniform prior around the ~34-day value (as reported by L21) on the third sinusoid accounting for the activity leads to deteriorated estimates of the previously presented planet parameters (~2-sigma detection of planet c) and a multi-modal posterior for the period of the third sinusoid. While unfortunate, this is not so surprising given that the two sets of observations are separated by over a year so the processes describing the star’s behaviour have evidently changed in this time due to the dynamic configuration of the activity regions across the stellar surface (see e.g., Barragán et al. 2021).

As can be seen in our periodogram analysis (Sect. 3.2), adding the new RV data suggests a stellar rotation period of ~21 days. Initially, we approached this result sceptically as we considered it probable that it is affected by both nuisance signals dominating the second half of the data, and the two time gaps in data acquisition. For this reason, and to try and understand the new RV data set better, we used a pyaneti-implemented multidimensional Gaussian process (multi-GP) regression following Rajpaul et al. (2015), to analyse all data, as well as the first and second observing seasons separately. Given the periodicity in the behaviour of the stellar activity evident in Fig. 4, we consider the quasi-periodic (QP) covariance function to be a good choice to model the stochastic part of the signal. The capability of GP regression with a QP kernel has been demonstrated both theoretically (e.g. Stock et al. 2023; Nicholson & Aigrain 2022) and practically (e.g. Rajpaul et al. 2016; Zicher et al. 2022; Deeg et al. 2023) in the literature, while the superiority of the multi-GP over 1D GP regression modelling was recently demonstrated by Barragán et al. (2023).

We tested modelling the RVs alongside different activity indicators and found that the dLW calculated by SERVAL behaves most closely like the RVs. Furthermore, while similar to the CCF-derived FWHM, the dLW is often a better choice for cool active stars as it does not suffer from negative effects in the absence of a good CCF (Zechmeister et al. 2018; Zicher et al. 2022). We thus proceeded to work with the dLW and the TERRA RVs. We started by placing a wide uniform prior around ~30 days for the activity signal in a set-up similar to, for example, Georgieva et al. (2021) for all three cases, which led to a non-detection of one or both planets. We thus finally turned to the ~21-day signal and, placing a Gaussian prior of ±1 day, used it to model the stellar rotation period in the multi-GP setup described above. This gave < 3-sigma results for both planets.

An important reason for the above-described outcomes is related to one key assumption of the multi-GP approach – that the same underlying GP can describe both the RVs and the activity indicator(s) time series, as well as the different observing seasons. In other words, GP regression relies on the fact that the data points are correlated. In a case such as the present one where the data sets are not well correlated, it is not surprising that this modelling approach does not perform well. Perhaps a future implementation that allows the GP hyperparameters to be set up in a more flexible way so as to describe stellar behaviour that has changed between observing seasons could lead to more successful results in cases such as this.

Fortunately, we find our solution in our novel approach – the SN-fit extraction of RV data, Sect. 4.2 – and its consequent results (Sect. 5.1), which gives realistic error bars, avoiding any bias issues, as discussed in Sect. 6.

Fig. 4 Generalised Lomb-Scargle periodogram of the complete HARPS RV data set and different activity indexes derived with SERVAL. The horizontal dashed lines mark the FAP level at 0.1%. The vertical yellow band marks the rotation frequency of TOI-776 at 0.0473 day−1 (Prot ≈ 21.1 days) and its 3σ uncertainty. The vertical dashed red and blue lines mark the orbital frequencies of TOI-776 b and c, respectively.

4.2 SN-fit-based radial velocity extraction

We re-reduced all of the HARPS data by using its dedicated data reduction software (DRS; Lovis & Pepe 2007). The DRS first cross-correlates the extracted Echelle spectra with a numerical mask (selected to be the closest to the stellar spectral type, M2 in our case; see Baranne et al. 1996; Pepe et al. 2002) to get one CCF per observation.

Unlike the standard approach of fitting a normal (N) function onto each CCF, following Simola et al. (2019) we performed the CCF-fitting using a SN function (Azzalini 1985). In addition to a location and a scale parameter (the counterparts of the Gaussian mean and standard deviation, respectively), we recall that the SN has a further degree of freedom (hereafter denoted with γ), which quantifies its skewness. In this way we can directly account for the intrinsic asymmetry of the CCF within the fit, without implementing a separate procedure to estimate it (e.g. the Bisector Span computation, Queloz et al. 2001). After performing the SN fit, our routine outputs the median of the best-fit function (i.e. the radial velocity measure, $\overline{RV}$), its full width at half maximum (FWHM_SN), its contrast (A), and its skewness (γ); these data are listed in a table at the CDS. We refer the reader to Simola et al. (2019) for a broad discussion about the specific output choices and the advantages of an SN fit over an N fit.

As the width and asymmetry of the CCF and their temporal evolution are typical tracers of the stellar activity (see e.g. Hatzes 1996; Queloz et al. 2001, 2009; Figueira et al. 2013; Simola et al. 2019; Bonfanti et al. 2023), the set of hyperparameters (FWHM_SN, A, γ) together with the time, t, constitute the basis vector against which to detrend the $\overline{RV}$ measurements to remove the activity component, RV_*. The RV detrending was performed within the LC+RV joint fit (see below, Sect. 5.1) using a polynomial baseline of the following form:

$$R{V}_{\star }={\beta }_0+\underset{k=1}{\overset{k_t}{{{\displaystyle \sum }}^{\text{}}}}{\beta }_{k,t}{t}^k+\underset{k=1}{\overset{k_F}{{{\displaystyle \sum }}^{\text{}}}}{\beta }_{k,F}{\text{FWHM}}_{\text{SN}}{}^k+\underset{k=1}{\overset{k_A}{{{\displaystyle \sum }}^{\text{}}}}{\beta }_{k,A}{A}^k+\underset{k=1}{\overset{k_{\gamma }}{{{\displaystyle \sum }}^{\text{}}}}{\beta }_{k,\gamma }{\gamma }^k,$$(1)

where the ß parameters are the polynomial coefficients, while (k_t, k_F, k_A, kγ) defines the polynomial order of the regression versus time, t, FWHM_SN, A, and γ, respectively.

Aware that the stellar activity of an M dwarf may dramatically impact the quality of our RV data, we also explored the possibility of a chunk-wise-based interpolation (e.g. the breakpoint method of Simola et al. 2022). The Bayesian Information Criterion (BIC; Schwarz 1978) did not support this however. Hence, we applied a unique polynomial detrending baseline in the form of Eq. (1) to the entire time series. The final phase-folded and detrended RV time series of both planets are displayed in Fig. 5.

5 Results

5.1 LC and RV joint analysis

For the joint TESS +CHEOPS transit and SN-fit-based RV modelling we used the MCMCI code (Bonfanti & Gillon 2020).

To save computational time, we switched off the simultaneous interaction with the stellar isochrones and tracks as the star was already well characterised using different approaches, as detailed in Sect. 3.1.

We imposed normal Gaussian priors on the stellar T_eff [Fe/H], M_*, and R_*, which has a double goal. First, the mean stellar density, ρ_*, inferred from M_* and R_*, constrains the transit parameters via Kepler’s third law. Second, stellar parameters are needed to get the limb darkening (LD) parameters for both the TESS (TE) and CHEOPS (CH) bandpasses following interpolation in the tables derived from ATLAS9 models by the code of Espinoza & Jordán (2015). Assuming a quadratic LD law (e.g. Claret 2000, and references therein), we estimated u_1,CH = 0.247 ± 0.041, u₂,_CH = 0.420 ± 0.030, u_1,TE = 0.177 ± 0.033, and u_2,TE = 0.421 ± 0.023, which define the normal priors of the LD coefficients that are jump parameters within the MCMC scheme. The LD posterior values are then reported in Table 4.

The other jump parameters, namely the transit depth, ${\left(\frac{R_p}{R_{\star }}\right)}^2$, the impact parameter, b, the orbital period, P, the transit timing, T₀, and the RV semi-amplitude, K, were subject to wide uniform priors (bounded by physical bounds only). Instead, modelling the eccentricity, e, with the beta distribution (Kipping 2013), with a wide uninformative prior leads to highly eccentric orbits for TOI-776 c (see Figs. A.3 and A.4). This would play against the system stability as reported by L21. Nonetheless, the orbits of planets in multi-transiting systems are not necessarily circular, even if they are expected to have low eccentricity values (Van Eylen & Albrecht 2015; Xie et al. 2016; Hadden & Lithwick 2017). Therefore, we imposed uniform priors on both $\sqrt{e}\cos \omega$ and $\sqrt{e}\sin \omega$, but bounded so as to imply e ≲ 0.2. This upper limit of e was set according to the analysis by L21, who followed similar considerations and found a 3σ upper limit on the eccentricities of both exoplanets equal to 0.18.

Both the LC- and RV-detrending were done simultaneously within the MCMC scheme and were based on polynomials, as shown for example in Eq. (1) for the RV side. To choose the best polynomial baseline (that is, the set of polynomial orders to be attributed to the vectors of additional parameters complementing the LC and RV time series) we performed several MCMCI runs, varying each time the polynomial order of the hyperparameters to establish the set of detrending orders favoured by the BIC. The adopted baseline is summarised in Table A.1.

Once the setup was completed, we launched a first preliminary MCMCI run to properly re-scale the photometric errors, as detailed in Bonfanti & Gillon (2020). After that, we performed a final MCMCI analysis comprising three independent runs of 200 000 steps each (burn-in phase equal to 20%) to check the convergence through the Gelman-Rubin (GR) test (Gelman & Rubin 1992). All the jump parameters converge nicely according to the GR test, with the posterior outcomes that are listed in Table 4.

5.2 Interior and atmospheric modelling

We attempted to perform an analysis of the internal structure of the two planets in the TOI-776 system. Here, we followed closely the global Bayesian modelling described in Leleu et al. (2021).

Briefly, this model is based on fitting the observed properties of the planets. Here, the planet-star radius ratio, the RV semi-amplitude, the orbital period, the stellar mass, radius, age, T_eff, and the photospheric abundances, [Si/Fe] and [Mg/Fe], are priors. It should be noted that we assumed that [Mg/Fe] = [Si/Fe] since we could not determine the Si abundance. The Bayesian analysis then relied on a forward model that computed the expected planetary radius and bulk internal structure as a function of the hidden parameters. The hidden parameters were, for each planet, the masses of solids (everything except the H or He gas), the mass fractions of the core, mantle, and water, the mass of the gas envelope, the Si/Fe and Mg/Fe mole ratios in the planetary mantle, the S/Fe mole ratio in the core, and the equilibrium temperature. We also assumed, in the forward model, that the planets are fully differentiated and consist of a core (Fe and S), a mantle (Si, Mg, Fe, and O), a pure water layer, and a H and He layer. For further details we refer to Leleu et al. (2021).

The posterior distributions of the main planetary hidden parameters indicate that both planets have a very small fraction (in terms of planetary mass) of H and He gas, and, albeit with large errors, that the mass of gas in the innermost planet is larger than that in the outermost one ($\log M_{\text{gas }}(b)=-{4.34}_{-0.44}^{+1.95}$, $\text{log\hspace{0.17em}}{M}_{\text{gas}}(c)=-{5.81}_{-5.54}^{+3.49}$. The fraction of water, on the other hand, is essentially unconstrained in our model.

We also attempted to model the atmospheric escape as a function of age, using the PASTA planetary atmospheric evolution code (Bonfanti et al. 2021b), which is an updated version of the original code presented by Kubyshkina et al. (2019b,a). We do find, however, that the planetary parameters (especially age) are not constrained well enough (with the exception of the radii and P_rot) to be able to reach any sensible conclusion about the atmospheric escape.

Fig. 5 Phase-folded and detrended RV time series of TOI-776 b (top panel) and TOI-776 c (bottom panel). Data are represented through black markers, while the red line is the best-fit model; the residuals are also shown.

6 Discussion

Combining TESS LCs, the photometry from CHEOPS, and the stellar host characterisation (Sect. 3.1), we reach a precision for the radii of planets b and c of 4.3% and 4.0%, respectively. This marks an improvement of almost a factor of two with respect to the results of L21. We note that we have thus approached the specific precision that is planned for planetary radii measurements with ESA’s PLATO mission later this decade. This was done with great effort in this one case, whereas PLATO will reach the same precision in tens of thousands of cases, simultaneously. We also note that the planetary radii approach the precision with which the stellar radii is known (3.0% – itself a very good value).

This implies that the current data are very close to what is currently (before PLATO’s combined asteroseismology and transit observations) possible.

Using all available RV data, on the dynamical side, we detected the RV semi-amplitudes of planets b and c at the 3.2σ and 2.7σ levels (see Table 4), respectively, which are slightly lower than what was obtained by L21. However, we find a few weaknesses in the approach followed by L21.

First of all, besides the two Keplerian models of TOI-776b and c, L21 further introduced a sinusoid based on the stellar rotation period for catching the activity-related signal of the host within the RV time series. This approach is discouraged in the literature as it is hard to physically justify that activity-induced RV signals are strictly sinusoidal (e.g. Lanza et al. 2001; Brinkworth et al. 2005). In general, those signals may also be quasi-periodic or aperiodic and subtracting sinusoids from the RV timeseries may introduce spurious harmonics that would bias the results (e.g. Pont et al. 2011; Tuomi et al. 2014; Rajpaul et al. 2015).

Secondly, as detailed in Sect. 3.2, the analysis, taking into account all the activity indices, leads us to infer a stellar rotation period, $P_{\text{rot}}={21.129}_{-00058}^{+0.062}$ days. The significant difference with respect to the rotation period derived by L21 ($P_{\text{rot},\text{L}21}={34.4}_{-2.0}^{+1.4}\text{\hspace{0.17em}\hspace{0.17em}d}$) suggests that subtracting a sinusoid with period P_rot,L21 from the RV time series makes the previously described scenario even worse.

Our novel SN-extraction of RV data provides the basic functions for detrending the RV data without introducing ad-hoc hypotheses. The theoretical model we employed is the simplest, being composed of the two Keplerian signals of TOI-776 b and c. Therefore, the data we provide to the literature are likely not affected by any bias and their error bars appear genuine.

In order to put our results for the planets of TOI-776 into context, we downloaded the parameters of all the known exoplanets with high-precision RV and radii data from the NASA Exoplanet Archive¹². After filtering those exoplanets with masses from RV observations and radii from transits that have a precision of at least 45 % and 15 %, respectively¹³, out of a total of 5502 exoplanets (August 28 2023), we were left with 814 planets in 708 systems. In particular, out of these, only 61 planets in 43 systems orbit M dwarf stars (here defined as having a T_eff in the interval 2380–3850 K). By additionally requiring that the planets should be smaller than 5 R_⊕, these numbers dropped to 56 planets in 38 systems, which constituted our reference sample of M dwarf exoplanets.

The bulk (mass, radius, and density) physical parameters of our reference sample of planets, colour-coded with instellation, are displayed in Figs. 6 and 7. In detail, in Fig. 6 we show the density-radius diagram of the M dwarf exoplanets along with the interior models from Zeng et al. (2019) overlaid.

In this figure we see hints of four groups of planets, separated by their intrinsic density. What we designate as Group I have densities larger than pure rock (i.e. above the solid brown line), and radii lower than ~ 1.4 R_⊕. The next group (II) is mainly found to have densities lower than that represented by a mix of 50% Si and 50% H₂O (i.e. below the dashed blue line) and radii larger than about 2 R_⊕. Within Group II we have progressively more and more water as we go towards higher radii.

Between the solid brown and the dashed blue line (densities between pure rock and the 50–50 mix of Si+H₂O), there is a small number of planets with radii of larger than 1.6 R_⊕: Group III. Our TOI-776 planets belong to this latter group, having densities of about 4 g cm⁻³.

Looking at Fig. 6 in more detail, Group I appears to be separated from Groups II and III by a radius gap, which can be described by a slanted line going between 1.3 R_⊕ and 1.5 R_⊕ in radius and between 3.5 gc m⁻³and 8g cm⁻³ in density. This feature is similar to that reported by Fulton et al. (2017) for Kepler planets.

Finally, separated by another gap, in the bottom right corner of Fig. 6, we have a handful of planets (including our own Uranus and Neptune; square markers) whose densities indicate a higher H + He component (Group IV).

Based on a similar sample (but with different precision requirements in terms of mass and radius determination), Luque & Pallé (2022) also identified three of the same groups (I, II and IV – as can be seen by a comparison with that paper with our Figs. 6 and 7). In our density radius plot, we see a clear division between planets with densities between pure rock and 50% Si and 50% H₂O compositions.

Among all the M-dwarf exoplanets known up to now (regardless of the mass knowledge), only 36 in 23 systems are characterised by a radius precision ≲4% like our TOI-776 b and c. Within this subgroup, only 17 planets have radii ≲2 R_⊕, including TOI-776 b and c, and the seven planets orbiting Trappist-1. Out of these 17 planets, a total of nine planets have masses determined via transit timing calculations. Our calculations here thus increase the sample of planets orbiting M dwarfs in preparation for further studies.

It is clear that many more planets orbiting all types of host stars are required to be characterised with a precision higher than the current one (especially as what concerns masses) in order to properly assess planetary demographics, formation, and evolution. This makes the scientific case for space missions like CHEOPS and PLATO extremely compelling.

Fig. 6 Radius-density plot of planets with radii < 5 R⊕ orbiting M dwarfs (here defined as having Teff = 2380–3850 K) with 45% and 15% uncertainties or lower in mass and radius, respectively (57 planets in 38 systems). All planets except nine (the Trappist-1 and K2-146 systems) have masses from RV measurements. The density gap referred to in the text is clearly seen between 1.3 R⊕ and 1.5 R⊕ in radius and between 3.5 g cm−3and 8 g cm−3in density. The planets are colour-coded with instellation. The TOI-776 planets are marked with star symbols (planet b to the left and planet c to the right). The squares are (from left to right) the Solar System planets Mars, Venus, Earth, Neptune, and Uranus. Interior models from Zeng et al. (2019) are plotted as listed in the legend.

Fig. 7 Mass-radius plot for the same planets and uncertainties as in Fig. 6, colour-coded with instellation. Notation as in Fig. 6.

7 Conclusions

1. We investigated TOI-776, merging old and new TESS data with dedicated CHEOPS photometry. We achieve a precision of 4% for the radii of the two planets, equal to a factor of two greater than previously;

2. This level was reached with two and three CHEOPS visits, respectively, for planets b and c, using a 32-cm telescope on a very active 11.5 magnitude star. This clearly demonstrates that PLATO with its significantly larger collecting area, its long times on targets, and consequently many observed transits, is likely to supersede its required precision – especially for a large number of red dwarfs;

3. The refined radii, together with the masses, allow us to definitively classify the planets as sub-Neptunes or (because of the relatively high density) as something in between such objects and super-Earths;

4. Placing the two planets, TOI-776 b and c, together with other planets orbiting red dwarfs and with a high radius and RV accuracy, into a density versus radius diagram (Fig. 6), we see a strong indication of the radius gap between planets denser than pure rock models and those with a lesser density. This gap is found at smaller radii than for solar-like stars;

5. We have characterised the increase in activity of TOI-776 since the data of L21 was obtained. Using our results we determined the rotation period of TOI-776 and find it to be significantly (~ 2/3) shorter than the previous value quoted by L21;

6. We have demonstrated a possible method of treating RV and LC data from exoplanets transiting stars that are very active and/or changing the level of activity during the period of observation, in order to retrieve the best planetary parameters;

7. Given that the planets orbiting TOI-776 appear to be very interesting and potentially important targets for further analyses, we consider the present investigation, with its significantly improved precision in the planetary radii, to be contributing necessary data. Taking into account also the changed behaviour of the host star over a relatively short time span, we demonstrate that our results will be important when considering future investigations.

Acknowledgements

CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledges the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission. Funding for the TESS mission is provided by NASA’s Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. This work is done under the framework of the KESPRINT collaboration (http://www.kesprint.science). KESPRINT is an international consortium devoted to the characterisation and research of exoplanets discovered with space-based missions. It is based in part on observations made with ESO Telescopes at the La Silla Observatory (Chile) under programme 106.21TJ.001 (P.I. D. Gandolfi). This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. M.F., I.Y.G., C.M.P., A.Br., and J.K. gratefully acknowledge the support of the Swedish National Space Agency (DNR 177/19, 174/18, 65/19, 2020-00104). A.J.M. gratefully acknowledges the support of the SNSA funding (Career grant 120/19C). T.G.W., A.C.C., and K.H. acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant ST/R003203/1. Y.A. and M.J.H. acknowledge the support of the Swiss National Fund under grant 200020.172746. D.G. and L.M.S. gratefully acknowledge financial support from the CRT foundation under Grant 2018.2323 “Gaseous or rocky? Unveiling the nature of small worlds”. D.G., M.F., X.B., S.C., and J.L. acknowledge their roles as ESA-appointed CHEOPS science team members. A.De. and D.E. acknowledge 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). A.De., A.Le., acknowledge support from the Swiss National Centre for Competence in Research “PlanetS” and the Swiss National Science Foundation (SNSF). A.H., Sz.C., M.E., and H.R. acknowledge support by DFG grants HA 3279/12-1 and RA 714/14-1 within the DFG Schwerpunkt SPP 1992, “Exploring the Diversity of Extrasolar Planets”. R.L. acknowledges funding from University of La Laguna through the Margarita Salas Fellowship from the Spanish Ministry of Universities ref. UNI/551/2021-May 26, and under the EU Next Generation funds. S.H. gratefully acknowledges CNES funding through the grant 837319. The MOC activities have been supported by the ESA contract no. 4000124370. S.C.C.B. and V.A. acknowledge support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004 and IF/00650/2015/CP1273/CT0001, respectively. 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; PTDC/FIS-AST/28987/2017, POCI-01-0145-FEDER-028987. S.G.S. acknowledges support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants PID2019-107061GB-C66, ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia Maria de Maeztu-Centro de Astrobiologia (INTA-CSIC), as well as the sup- port of the Generalitat de Catalunya/CERCA programme. This project was supported by the CNES. L.D. is an F.R.S.-FNRS Postdoctoral Researcher. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. K.G.I. is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme. She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection for the programme. G.L. acknowledges support by CARIPARO Foundation, according to the agreement CARIPARO-Università degli Studi di Padova (Pratica n. 2018/0098). This work was granted access to the HPC resources of MesoPSL financed by the Région Île-de-France and the project Equip@Meso (reference ANR-10-EQPX- 29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. M.L. acknowledges support from the Swiss National Science Foundation under grant number PCEFP2.194576. P.F.L.M. acknowledges support from STFC grant number ST/M001040/1. L.D.N. thanks the Swiss National Science Foundation for support under Early Postdoc. Mobility grant P2GEP2.200044. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). I.R. acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B- C33, as well as the support of the Generalitat de Catalunya/CERCA programme. This project has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, the MTA-ELTE Lendület Milky Way Research Group and the City of Szombathely under Agreement No. 67.177-21/2016. This research received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n 803193/BEBOP), and from the Science and Technology Facilities Council (STFC.; grant ST/S00193X/1). V.V.G. is an F.R.S-FNRS Research Associate. We thank the Mathworks community and, in particular, Peter (2023, https://www.mathworks.com/matlabcentral/fileexchange/42985-break-x-axis).

Appendix A Supplementary material

Fig. A.1 Raw TESS LCs. Left: Sector 10. Right: Sector 37.

Fig. A.2 Raw CHEOPS LCs shown in chronological order of observation from top to bottom.

Fig. A.3 Analysis comparing output eccentricity for TOI-776 b (left) and TOI-776 c (right) as a function of the prior on maximum eccentricity

Fig. A.4 Analysis comparing output planetary mass for TOI-776 b (left) and TOI-776 c (right) as a function of the prior on maximum eccentricity

References

All Tables

All Figures

Fig. 1 Phase-folded LCs showing the transits of TOI-776b (first column) and TOI-776c (second column). Top panels: detrended CHEOPS LCs; the best-fit transit models are shown as red lines together with their corresponding residuals. Bottom panels: same as top, but for the TESS LCs. Data are shown in the nominal cadence modes (blue dots) and also binned to 10 min (black markers) with their associated error bars. We note that the transit durations and depths as measured by TESS and CHEOPS are identical within uncertainties. The transit depths found by CHEOPS are 909 ppm for planet b and 1177 ppm for planet c.

In the text

Fig. 2 Time series of the activity indicators. The complete RV data set has been used. The thick red lines mark the best-fitting line fit. We note the change in the activity level between the two seasons significantly detected in the contrast (first panel), the DLW (fourth panel), and the Ha index (seventh panel). The time is given as BJDTDB = Barycentric Julian Dates in Barycentric Dynamical Time − 2 450 000.

In the text

	Fig. 3 Whole RV time series: the break on the temporal axis avoids the big gap in the data collection and improves the visualisation. Top panel: undetrended RV data with the model of the full Keplerian signal + stellar activity superimposed in blue. Bottom panel: detrended RV data with the model of the full Keplerian signal superimposed in red.
In the text

Fig. 4 Generalised Lomb-Scargle periodogram of the complete HARPS RV data set and different activity indexes derived with SERVAL. The horizontal dashed lines mark the FAP level at 0.1%. The vertical yellow band marks the rotation frequency of TOI-776 at 0.0473 day−1 (Prot ≈ 21.1 days) and its 3σ uncertainty. The vertical dashed red and blue lines mark the orbital frequencies of TOI-776 b and c, respectively.

In the text

	Fig. 5 Phase-folded and detrended RV time series of TOI-776 b (top panel) and TOI-776 c (bottom panel). Data are represented through black markers, while the red line is the best-fit model; the residuals are also shown.
In the text

Fig. 6 Radius-density plot of planets with radii < 5 R⊕ orbiting M dwarfs (here defined as having Teff = 2380–3850 K) with 45% and 15% uncertainties or lower in mass and radius, respectively (57 planets in 38 systems). All planets except nine (the Trappist-1 and K2-146 systems) have masses from RV measurements. The density gap referred to in the text is clearly seen between 1.3 R⊕ and 1.5 R⊕ in radius and between 3.5 g cm−3and 8 g cm−3in density. The planets are colour-coded with instellation. The TOI-776 planets are marked with star symbols (planet b to the left and planet c to the right). The squares are (from left to right) the Solar System planets Mars, Venus, Earth, Neptune, and Uranus. Interior models from Zeng et al. (2019) are plotted as listed in the legend.

In the text

	Fig. 7 Mass-radius plot for the same planets and uncertainties as in Fig. 6, colour-coded with instellation. Notation as in Fig. 6.
In the text

	Fig. A.1 Raw TESS LCs. Left: Sector 10. Right: Sector 37.
In the text

	Fig. A.2 Raw CHEOPS LCs shown in chronological order of observation from top to bottom.
In the text

	Fig. A.3 Analysis comparing output eccentricity for TOI-776 b (left) and TOI-776 c (right) as a function of the prior on maximum eccentricity
In the text

	Fig. A.4 Analysis comparing output planetary mass for TOI-776 b (left) and TOI-776 c (right) as a function of the prior on maximum eccentricity
In the text