key: cord-0435402-iz3fi9fq
authors: Luque, R.; Serrano, L. M.; Molaverdikhani, K.; Nixon, M. C.; Livingston, J. H.; Guenther, E. W.; Pall'e, E.; Madhusudhan, N.; Nowak, G.; Korth, J.; Cochran, W. D.; Hirano, T.; Chaturvedi, P.; Goffo, E.; Albrecht, S.; Barrag'an, O.; Briceno, C; Cabrera, J.; Charbonneau, D.; Cloutier, R.; Collins, K. A.; Collins, K. I.; Col'on, K. D.; Crossfield, I. J. M.; Csizmadia, Sz.; Dai, F.; Deeg, H. J.; Esposito, M.; Fridlund, M.; Gandolfi, D.; Georgieva, I.; Glidden, A.; Goeke, R. F.; Grziwa, S.; Hatzes, A. P.; Henze, C. E.; Howell, S. B.; Irwin, J.; Jenkins, J. M.; Jensen, E. L. N.; K'abath, P.; Kidwell, R. C.; Kielkopf, J. F.; Knudstrup, E.; Lam, K. W. F.; Latham, D. W.; Lissauer, J. J.; Mann, A. W.; Matthews, E. C.; Mireles, I.; Narita, N.; Paegert, M.; Persson, C. M.; Redfield, S.; Ricker, G. R.; Rodler, F.; Schlieder, J. E.; Scott, N. J.; Seager, S.; vSubjak, J.; Tan, T. G.; Ting, E. B.; Vanderspek, R.; Eylen, V. Van; Winn, J. N.; Ziegler, C.
title: A planetary system with two transiting mini-Neptunes near the radius valley transition around the bright M dwarf TOI-776
date: 2020-09-17
journal: nan
DOI: nan
sha: c28d6e6660d12745957f4df11eb6b2e9bbc765ec
doc_id: 435402
cord_uid: iz3fi9fq

We report the discovery and characterization of two transiting planets around the bright M1 V star LP 961-53 (TOI-776, J=8.5mag, M=0.54+-0.03Msun) detected during Sector 10 observations of the Transiting Exoplanet Survey Satellite (TESS). Combining them with HARPS radial velocities, as well as ground-based follow-up transit observations from MEarth and LCOGT telescopes, we measured for the inner planet, TOI-776b, a period of 8.25d, a radius of 1.83+-0.11Re, and a mass of 4.66+-0.97Me; and for the outer planet, TOI-776c, a period of 15.66d, a radius of 2.06+-0.13Re, and a mass of 6.1+-1.5Me. The Doppler data shows one additional signal, with a period of 35d, associated with the rotational period of the star. The analysis of fifteen years of ground-based photometric monitoring data and the inspection of different spectral line indicators confirm this assumption. The bulk densities of TOI-776b and c allow for a wide range of possible interior and atmospheric compositions. However, both planets have retained a significant atmosphere, with slightly different envelope mass fractions. Thanks to their location near the radius gap for M dwarfs, we can start to explore the mechanism(s) responsible for the radius valley emergence around low-mass stars as compared to solar-like stars. While a larger sample of well-characterized planets in this parameter space is still needed to draw firm conclusions, we tentatively estimate that the stellar mass below which thermally-driven mass loss is no longer the main formation pathway for sculpting the radius valley is between 0.63 and 0.54Msun. Due to the brightness of the star, the TOI-776 system is also an excellent target for the JWST, providing a remarkable laboratory to break the degeneracy in planetary interior models and to test formation and evolution theories of small planets around low-mass stars.

Exoplanets with masses between those of Earth and Uranus are characterised by a broad range of measured bulk densities (e.g., Hatzes & Rauer 2015) . A low density suggests the presence of an extended H/He-envelope around a solid core. On the contrary, if the density is high, the exoplanet is considered to be fully rocky or enriched in light elements (e.g., water, methane, ammonia). The absence of an envelope might be the result of two different formation mechanisms. In the first case, the planet forms in a gas-poor inner proto-planetary disk without a thick H/Heenvelope. In the second case, the envelope disappeared as a consequence of erosion processes driven by the stellar X-ray+EUV (XUV)-radiation, coupled with non-thermal processes, such as ion pick-up (Chen & Rogers 2016; Jin et al. 2014; Jin & Mordasini 2018; Kislyakova et al. 2013 Kislyakova et al. , 2014 Lammer et al. 2012 The erosion rate becomes faster if the planetary surface gravity decreases and the amount of XUV-radiation that the planet receives increases. In addition, the intensity of XUV-radiation depends on the orbital semi-major axis and on the stellar activity level. The XUV-radiation is particularly high at young ages and then it declines as a result of age, mass and stellar rotation (Walter et al. 1988; Briceno et al. 1997; Tu et al. 2015) . A star that begins its life rapidly rotating will suffer a more rapid decline in rotation than a star that was initially a slow rotator. Thus, for a several Gyr old star, understanding its original activity level is challenging. The presence, or absence, of a hydrogen-rich envelope in a system containing just one planet can thus be equally explained by assuming that the host star was either a slow or rapid rotator when it was young. Systems containing more than one planet are necessary to test the theory of atmospheric erosion, because the origin of all the planets of a system should be explained with a unique evolutionary history of the host's XUVradiation (Owen & Campos Estrada 2020) .

On the other hand, the amount of XUV-radiation also depends on the stellar type. The XUV luminosities of young Gand M-stars are similar to each other. The average X-ray luminosity of G-stars is 10 29 erg s −1 , while in the case of M dwarfs, Article number, page 1 of 24 arXiv:2009.08338v1 [astro-ph.EP] 17 Sep 2020 A&A proofs: manuscript no. main the 50 Myr stars in α-Per, for example, have luminosities of 10 28 erg s −1 . The main difference is that M dwarfs remain in the high activity phase for up to 2 Gyr , a much longer time compared to the 300 Myr of G-stars (Güdel et al. 2004 ). This makes M dwarfs preferred targets to study planetary systems that have experienced significant stellar XUV irradiation. Another advantage of M dwarfs is their small size, which makes it easier to detect smaller transiting planets. The paucity of close-in planets around mid-K to mid-M dwarfs between approximately 1.4 and 1.7 R ⊕ , known as the radius valley, marks the transition between rocky planets and sub-Neptunes orbiting low-mass stars. As such, M dwarfs multi-planetary systems which include sub-Neptunes and/or rocky planets represent an ideal benchmark for testing the theory of atmospheric erosion.

Gas-poor formation provides an alternative to explain the absence of H/He envelopes in some low-mass planets, since the erosion scenario presents some issues. For instance, if a closein 10 M ⊕ rocky planet forms while there is still a gaseous disk, its mass is high enough to undergo runaway accretion and become a Jupiter-type planet. The detection of close-in Jupitermass planets, at least in A-stars, shows that it is hard to reconstruct a mechanism which transforms a Jupiter into a rocky super-Earth since any working physical process should be able to completely strip off the H/He atmosphere. On the other hand, stars hosting hot Jupiters have high metallicities, while rocky planets are equally distributed between metal-poor and metalrich stars (Winn et al. 2017) . Thus, there are two alternative scenarios within gas-poor formation models that could explain the existence of rocky super-Earths. Either the dust-to-gas ratio of the inner disk is 20 times higher than solar, or the gas accretion is delayed until just before the disk disperses (Lee et al. 2014; Lee & Chiang 2016) . Lopez & Rice (2018) proposed a statistical test that could allow us to understand the most likely formation history for super-Earths. If a high percentage of rocky planets are the evaporated cores of sub-Neptunes, the transition radius from rocky to sub-Neptune planets should decrease for longer orbital periods. On the contrary, if the gas-poor formation scenario is correct, the transition radius should increase with orbital period. Another methodology to test the formation theory of super-Earths requires studying the position of the radius valley for stars with different masses, thus of different stellar types. If the photoevaporation scenario is correct, the radius valley shifts towards planets of smaller radii for stars of lower mass. If, on the contrary, the gas-poor formation scenario is at work, the valley position is not affected by the stellar mass (Cloutier & Menou 2020) . However, since the radius valley represents the range of radii in which the transition between rocky planets and sub-Neptunes occurs, it is necessary to accurately determine the mass and radius of the planets to calculate the mass-fraction of their envelope and unveil their nature. Therefore, the ideal test to understand which model is more realistic between the gas-poor formation and the photo-evaporation consists of measuring the masses and radii of the planets close to, or inside, the radius valley, preferably in a multi-planetary system around low-mass stars. In this way, we can also constrain these models in a much better way than through the radius distribution alone.

As of today, there is a limited number of known multiplanetary systems which orbit early M-dwarfs (3084 K < T eff < 3759 K; M0 V to M5 V) and respect the condition (M p < 10 M ⊕ ) required to test the two mentioned formation theories. In particular, we are not aware of any system with three or more transiting planets with measured dynamical masses, while we identi-fied only three systems with two transiting planets (LHS 1140 , Dittmann et al. 2017 Ment et al. 2019; LTT 3780, Nowak et al. 2020; K2-146, Lam et al. 2020; Hamann et al. 2019) . This paucity of systems is inadequate for understanding the formation and evolution of planetary systems around M dwarfs. The discovery of each new system is thus important, especially if the host star is bright and the planets are close to the radius valley.

In this paper, we present the discovery of two transiting planets orbiting an M1 V star. The inner one has a period of 8.2 d and a radius of ∼ 1.8 R ⊕ ; thus, it is close to the radius valley. The outer planet has a period of 15.7 d and a radius of 2.1 R ⊕ , in the sub-Neptune regime. By measuring their masses, we explore whether these new planets are characterised by extended H/He envelopes. Since they orbit a relatively bright, nearby M dwarf, these new objects represent ideal targets for follow-up atmospheric studies. Smith et al. 2012; Stumpe et al. 2012 Stumpe et al. , 2014 to TESS. Figure 1 the remainder of this work we make use of the latter photometric data, shown in Fig. 2 .

On June 11, 2019, two transiting candidates orbiting LP 961-53 were announced in the TESS data public website 2 under the TESS Object of Interest (TOI) number 776. TOI-776.01 is a planet candidate with a period of 15.65 d, a transit depth of 1484 ± 127 ppm, and an estimated planet radius of 2.2 ± 0.6 R ⊕ ; while TOI-776.02 is a planet candidate with a period of 8.24 d, a transit depth of 1063 ± 104 ppm, and an estimated planet radius of 1.8 ± 1.4 R ⊕ . Both candidates passed all the tests from the TCE (Threshold Crossing Event) Data Validation Report (DVR; Twicken et al. 2018; Li et al. 2019 ): even-odd transits comparison, eclipsing binary (EB) discrimination tests, ghost diagnostic tests to help rule out scattered light, or background EB, among others. However, the vetting team at the TESS Science Office proposed the possibility that TOI-776.01 could be an EB, where the secondary transit is the primary transit of TOI-776.02 candidate. The ground-based follow-up observations discussed in the next Section refuted this scenario and confirmed the two announced candidates as bona-fide planets.

We observed the TOI-776 candidates as part of the TESS Followup Observing Program (TFOP) 3 . The goals of these groundbased photometric follow-up observations were to verify that the transits observed by TESS are on target, and to refine the transit ephemeris and depth measurements. We used the TESS Transit Finder, a customised version of the Tapir software package (Jensen 2013) , to schedule photometric time-series follow-up observations. We observed two transits of TOI-776.01 and three transits of TOI-776.02, as summarised in Table 1 (Irwin et al. 2015) at Cerro Tololo Inter-American Observatory (CTIO), Chile on June 1, 2019. Seven telescopes observed continuously from evening twilight until the target star set below airmass 2, using an exposure time of 10 s with all telescopes in focus. The target star was west of the meridian throughout the observation to avoid meridian flips.

Data were reduced following the standard procedures in Irwin et al. (2007) and Berta et al. (2012) with a photometric extraction aperture radius of r = 6 pix (5 on sky given the pixel scale of 0 . 84 pix −1 ). The light curve is shown in Fig. 3 , lower right. Due to the large variation in airmass and relatively red target star compared to the available field comparison stars, we found the light curve exhibited a small amount of residual second-order (color-dependent) atmospheric extinction, so the transit model was fitted including an extinction term (linear decorrelation against airmass).

One transit of TOI-776.01 and two transits of TOI-776.02 were observed with the 1.0 m telescopes in the Las Cumbres Observatory (LCOGT) telescope network (Brown et al. 2013 ). The 4096 × 4096 pix LCOGT SINISTRO cameras have an image scale of 0 . 389 pix −1 , resulting in a 26 × 26 field of view. The images were calibrated using the standard LCOGT BANZAI pipeline, and photometric data were extracted with AstroImageJ (Collins et al. 2017 ).

An ingress of TOI-776.01 was observed from the LCOGT node at CTIO on July 1, 2019 in the i filter, simultaneous with the MEarth-South observations mentioned above (Fig. 3 Woźniak et al. 2004) , and the Catalina surveys (Drake et al. 2014 ). The telescope location, instrument configurations, and photometric bands of each public survey were summarized in Table 1 of Díez Alonso et al. (2019) . All together, the measurements span a period of 15 yr. Additionally, TOI-776 is a candidate of the Super-Wide Angle Search for Planets (SuperWASP; Pollacco et al. 2006 ). Su-perWASP acquired more than 11 000 photometric observations, using a broad-band optical filter spanning three consecutive seasons from May to July 2006, January to June 2007, and January to June 2008. In order to detect long-term photometric modulations associated with the stellar rotation, we binned the data into one day intervals, resulting in 201 epochs.

The large pixel size of TESS increase the possibility of contamination by nearby sources that are not detected in the seeinglimited photometry or in Gaia DR2. Close companions can dilute the transit depth and thus alter the measured planet radius, or lead to false positives if the companion is itself an EB (e.g. Ciardi et al. 2015) . We thus searched for companions by collecting adaptive optics (AO) and speckle images of TOI-776 using 4 and 8 m class telescopes, providing robust limits on the presence of companions and the level of photometric dilution.

Gemini/NIRI On June 15, 2019, TOI-776 was observed using the adaptive optics near-infrared imager (NIRI) mounted on the 8.1 m Gemini North telescope at Mauna Kea, Hawai'i. We collected a total of 9×1.4 s images in the Brγ filter centered on 2.166 µm. We dithered the telescope between exposures, so that 4 http://pestobservatory.com/ the sky background can be constructed from the science frames themselves. After removing bad pixels, flat-fielding, and subtracting the sky background, we aligned the stellar position between frames and co-added the images. The sensitivity of our observations was calculated as a function of radius by injecting fake companions, and scaling their brightness such that they could be detected at 5σ. The contrast curve and image are shown in Fig. 4 . Only the central 4 ×4 are shown, but no companions are seen anywhere in the field, which has a field of view ∼13 ×13 .

VLT/NaCo On July 4, 2019, TOI-776 was observed in Brγ using the NAOS-CONICA AO instrument (NaCo), mounted at the Nasmyth A port of the 8 m UT1 Very Large Telescope (VLT) in Paranal, Chile. We collected a total of 9×10 s Brγ images. Data were reduced and analyzed using the same procedures as described above for the NIRI data, and no companions were found in the reduced image. The NaCo contrast curve is shown in Fig. 4. 3.3.2. Speckle imaging SOAR/HRCam On December 12, 2019, TOI-776 was observed in I band with a pixel scale of 0.01575 pix −1 using the HRCam imager, mounted on the 4.1 m Southern Astrophysical Research (SOAR) telescope at Cerro Tololo Inter-American Observatory, Chile. The data were acquired and reduced following the procedures described in Tokovinin (2018) and Ziegler et al. (2020) . The resulting reconstructed image achieved a contrast of ∆mag = 7.1 at a separation of 3 (see top panel of Fig. 5 ).

Gemini/Zorro On March 15, 2020, TOI-776 was observed using the Zorro speckle imager (Scott 2019) , mounted on the 8.1 m Gemini South telescope in Cerro Pachón, Chile. Zorro uses high speed electron-multiplying CCDs (EMCCDs) to simultaneously acquire data in two bands centered at 562 nm and 832 nm. The data were collected and reduced following the procedures described in Howell et al. (2011) . The resulting reconstructed image achieved a contrast of ∆mag = 7.8 at a separation of 1 in the 832 nm band (see bottom panel of Fig. 5 ). We note that at the distance of TOI-776, our Zorro speckle images cover a spatial range of 0.46 to 32 au around the star with contrasts between 5 to 8 mag. 

We obtained 29 high-resolution (R≈115000) spectra of TOI-776 using the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph mounted at the ESO 3.6 m telescope of La Silla Observatory, Chile (Mayor et al. 2003) . The observations were carried out as part of our large observing pro-gram 1102.C-0923 starting on February 5 and until March 23, 2020, when ESO observatories stopped the operations due to the COVID-19 pandemic. One spectrum was acquired under the program 60.A-9709. We used the second fiber of the instrument to monitor the sky background and we reduced the data with the HARPS data reduction software (DRS; Lovis & Pepe 2007) . To compute precise radial velocities and spectral Contrast curves from NIRI (orange) and NaCo (blue), and the central 4 ×4 of the NIRI image (inset). We rule out companions 6 mag fainter than TOI-776 beyond 250 mas, and 7.5 mag fainter beyond 900 mas. The NaCo observations have a slightly tighter inner working angle, while the NIRI observations reach a deeper sensitivity beyond 0.5 . diagnostics, we applied on the reduced data the codes serval (Zechmeister et al. 2018) and TERRA (Anglada-Escudé & Butler 2012). Both programs employ a template-matching algorithm that is better suited to derive precise radial velocities for M dwarfs, if compared to the cross-correlation function (CCF) technique implemented in the DRS. In the CCF technique, the line lists of M dwarfs used to define the binary mask are incomplete and they thus produce a CCF which is often a poor match for cool stars. The RVs have a median internal uncertainty of 1.5 m s −1 (1.5 m s −1 ) and a root mean square of 5.2 m s −1 (3.5 m s −1 ) around the mean value for the serval (TERRA) extractions, respectively. We report in Tables B.1 and B.2 of the Appendix the HARPS measurements, the extracted RVs and the associated uncertainties, Na i D, Na ii D, and Hα line indices from both programs together with the chromatic index (CRX) and differential line width (dLW) computed by serval, and the Mount Wilson S-index computed by TERRA. Gaidos et al. (2014) . For an all-sky sample of approximately 3000 M-or late K-type stars, the authors provide spectroscopically determined values of the spectral type, effective temperature and metallicity, which combined with empirical relations for cool stars, allow to estimate stellar radius, luminosity and mass. In particular, they measure that TOI-776 is a relatively inactive M1 V dwarf star with the stellar properties shown in Table 2 .

We carried out an independent analysis to improve the photospheric and fundamental parameters of TOI-776. We used SpecMatch-Emp (Yee et al. 2017) to empirically estimate the effective temperature, metallicity, and stellar radius by comparing the co-added HARPS high-resolution spectrum with a spectroscopic library of well-characterized stars. The results of this analysis are in agreement with the values of Gaidos et al. (2014) within the errors. Then, we derived the stellar radius and luminosity combining Gaia G, G BP , G RP photometry and 2MASS J, H, K s magnitudes with the spectroscopic parameters from the SpecMatch-Emp analysis and the Gaia parallax. We corrected the Gaia G photometry for the magnitude dependent offset using Eq. 3 from Casagrande & VandenBerg (2018) , and adopted a minimum uncertainty of 0.01 mag for the Gaia magnitudes to account for additional systematic uncertainties. We added 0.06 mas to the nominal Gaia parallax to account for the systematic offset found by Stassun & Torres (2018) ; Riess et al. (2018) ; Zinn et al. (2019) . Our best estimate of the stellar radius is consistent with the value from Gaidos et al. (2014) and in agreement with each of the radius estimates obtained independently using only one of the magnitudes. Finally, we computed the mass using the massradius relations for M dwarfs from Schweitzer et al. (2019) .

We also applied the methods of Reddy et al. (2006) to Gaia DR2 astrometry for TOI-776 to compute galactic U, V, W velocities in the local standard of rest and the probabilities of kinematic membership in galactic stellar populations. We found that TOI-776 has a probability of 96.3% of belonging to the thin Article number, page 6 of 24 R. Luque et al.: A multi-planetary system around TOI-776 disk population, which is in excellent agreement with the galactic population probabilities for this star in the recent catalog of Carrillo et al. (2020) . Additionally, using the code isochrones (Morton 2015), we estimated the age of TOI-776 to be loosely constrained between 2 to 10 Gyr. From the metallicity, age and kinematics given in Table 2 , we can conclude that TOI-776 is a relatively old member of the galactic thin disk population.

To determine the rotational period of the star, we used the publicly available photometric data for TOI-776. Using juliet (see more details about the algorithm in Sect. 5.2) we modeled the ASAS-SN, ASAS, NSVS, Catalina, and daily binned Super-WASP data with Gaussian processes (GPs). In particular, we adopted the quasi-periodic GP kernel introduced in Foreman- 

where τ = |t i −t j | is the time-lag, B and C define the amplitude of the GP, L is a timescale for the amplitude-modulation of the GP, and P rot is the rotational period of the modulations. As in Luque et al. (2019), we considered each of the five data-sets to have different values of B and C, in order to account for the possibility that different bands could have different GP amplitudes, while we imposed the timescale of the modulation and the rotational period as common parameters for all the data sets. In addition, we fitted for an extra jitter term for each photometric time series. We considered wide uninformative priors for the jitter, B, C, L, and a uniform rotation period prior between 10 and 100 d. Figure 6 shows the posterior samples of the GP hyperparameter P rot after fitting all the long-term monitoring ground-based photometry. The distribution is bimodal with peaks at 33 ± 1 d and 38±1 d, where the samples from the first peak have the highest likelihood. From this we can estimate that the stellar rotation of TOI-776 is between 30 to 40 d over the course of 15 yr. The 38 d peak may be an alias of the true 33 d rotation due to 1 yr window function in the photometry. Alternatively, the bimodal distribution of the P rot can be explained as a consequence of the stellar differential rotation coupled with the activity cycle (Rüdiger et al. 2014; Küker et al. 2019) . For early M dwarfs with rotational periods similar to TOI-776, the expected dynamo cycle time is between 3 to 6 yr (Küker et al. 2019) , thus detectable in our data. Additionally, assuming that this star is a solar-like rotator, the rotational velocity of the star decreases as the latitude increases. The two peaks correspond to two different groups of activity features, a bigger one, closer to the equator, which generates the first peak of the posterior distribution, and a smaller one, at a higher latitude, which produces the second peak. The opposite situation, with an anti-solar like rotator, is less likely, considering that TOI-776 is an adult star, still belonging to the main sequence.

Article number, page 7 of 24 A&A proofs: manuscript no. main 

We performed a frequency analysis of the HARPS serval/TERRA extracted measurements to search for the Doppler reflex motion induced by the two transiting planets discovered in the TESS light curve and to unveil the presence of additional signals associated with the star and/or other orbiting planets. Figure 7 shows the generalized Lomb Scargle (GLS; Zechmeister et al. 2009) periodograms of the HARPS RVs and activity indicators extracted with serval (blue lines) and with TERRA (red lines). The horizontal dashed lines mark the GLS powers corresponding to the 0.1, 1, and 5% false alarm probability 5 (FAP). The vertical dashed lines mark the orbital frequencies of the two transiting planets detected in the TESS light curve (f b = 0.121 d −1 and f c = 0.064 d −1 ) and the stellar signal at ∼0.03 d −1 (see below).

The upper panel of Fig. 7 displays the GLS periodogram of the HARPS RVs in the frequency range 0 -0.42 d −1 . The highest peak is found at 0.055 d −1 (FAP ≈ 0.3 %), which is close to the orbital frequency of TOI-776 c (f c = 0.064 d −1 ). Taking into account our frequency resolution 6 of 0.021 d −1 , the two frequencies are indistinguishable. This suggests that the highest peak seen in the periodogram of the HARPS RVs is the stellar reflex motion induced by the outer transiting planet TOI-776 c. The second highest peak is found at 0.129 d −1 (Fig. 7 , upper panel), which is close to the orbital frequency of TOI-776 b. However, this signal is an alias of the signal at 0.055 d −1 . The periodogram of the window function indeed shows a peak at 0.074 d −1 (highlighted with an arrow in the bottom panel of Fig. 7) , which is equal to the frequency spacing between the two highest peaks seen in the periodogram of the HARPS RVs.

We used the code pyaneti (Barragán et al. 2019) (Sect. 5.2.3) to subtract the Doppler signal of TOI-776 c from the HARPS RVs. We assumed a circular model (see also Sect. 5.2.2), fixing period and time of first transit to the TESS ephemeris, while allowing for the systemic velocity and RV semi-amplitude to vary. The periodogram of the RV residuals shows a broad peak centered around ∼0.04 d −1 with a FAP of about 10 %. Although the Doppler signal is not significant, the GLS of the CRX, dLW, Hα, and S-index activity indicators show also peaks at ∼0.04 d −1 , suggesting that this signal is caused by the presence of active regions appearing and disappearing from the visible stellar disk as the star rotates around its axis. It is worth noting that the peak at 0.130 d −1 is not observed in the GLS periodogram of the RV residuals, corroborating the interpretation that this peak is an alias of the dominant frequency detected in the periodogram of the HARPS data.

We removed the Doppler reflex motion of TOI-776 c and the activity-induced RV signal by jointly modeling the HARPS measurements with a circular Keplerian orbit and a sine curve. For TOI-776 c we followed the same procedure described in the previous paragraph. For the stellar signal we fitted for the phase, amplitude, and frequency. The latter was allowed to vary within a wide uniform prior centered around 0.04 d −1 . The GLS Table 3 . Model comparison of RV-only fits with juliet. The prior label N represents a normal distribution. The final model used for the joint fit is marked in boldface (see Sect. 5.2.2 for details about the selection of the final model).

Prior

with a uniform prior in P rot;GP,RV ranging from 5 to 50 d.

periodogram of the RV residuals displays a peak at 0.125 d −1 (FAP ≈ 11 %), which is very close to the frequency of the inner transiting planet TOI-776 b (f b = 0.121 d −1 ). We note that the activity indicators show also peaks close to the orbital frequency of TOI-776 b. Yet, those peaks are separated by 0.074 d −1 from the stellar signal at ∼0.04 d −1 . As such they are very likely aliases of the latter.

In this section, we use juliet (Espinoza et al. 2019) 

First, to constrain the properties of the transiting planets and use them for further analyses, we modeled the TESS, LCO, and MEarth photometry with juliet. We adopted a quadratic limb darkening law for TESS, since Espinoza & Jordán (2015) showed it was appropriate as well for space-based missions. The limb darkening parameters were then parametrized with a uniform sampling prior (q 1 , q 2 ), introduced by Kipping (2013). For LCO and MEarth transits, we used a more simple linear limb darkening law, because the lower data precision with respect toTESS prevents us from adopting a more complex law. Additionally, we followed the parametrization introduced in Espinoza (2018). In particular, for each transiting planet, rather than fitting for the planet-to-star radius ratio p = R p /R * and the impact parameter of the orbit b, we sampled from the uniform priors assigned to two parameters, r 1 and r 2 , which are connected to p and b with the equations (1)-(4) in Espinoza (2018) . r 1 and r 2 were shown in Espinoza (2018) to guarantee a full exploration of the physically plausible values in the (p, b) plane. We assumed as well circular orbits and fixed the TESS dilution factor to 1, based on our anal- Table 4 . Median and the 68% credibility intervals of the posterior distributions for each fitted parameter of the final joint model obtained for the TOI-776 system using juliet. Priors and descriptions for each parameter can be found in To account for the time-correlated noise in the light curve in Fig. 2 , even using the PDC-corrected SAP, we modeled the TESS photometry with the exponential GP kernel

where T GP,TESS is a characteristic timescale and σ GP,TESS is the amplitude of this GP modulation. For the LCO photometry, on the other hand, we used a linear model to detrend the data from airmass correlations. Our photometry-only analysis increases significantly the precision of the planet parameters with respect to the TESS DVR. The uncertainties in the period decreases by two orders of magnitudes which eases up future ground-and space-based followup efforts. The radii of the planets are determined to a precision better than 5%. In addition, the ground transit observations help to constrain the short-term transit time variations (TTVs) and allowing for a robust estimate of the mid-transit epoch.

Finally, we searched for an additional planets in the system by modeling a three-planet fit with the same priors as in Table A.1 for the transiting planets, and varying the period and mid-transit time of the third hypothetical planet. Our result significantly exclude the presence of any additional transits in the light curve (∆ ln Z = ln Z 2pl − ln Z 3pl > 7).

Even though the results of the RVs extraction slightly change whether we use serval or TERRA, the GLS analysis in both cases show the evidence of a stellar signal together with the RV trends associated with the transiting planets. To adequately describe the data, we considered several RV-only models and carried out a model comparison scheme as in Luque et al. (2019) . We used juliet, a code which efficiently computes the Bayesian log-evidence of each tested model and explores the parameter space using the importance nested sampling included in MultiNest (Feroz et al. 2009 ) via the PyMultiNest package (Buchner et al. 2014) . As discussed in Nelson et al. (2020) , this method outperforms other samplers in choosing robustly the best model for those with 3 or less planets. We considered a model to be moderately favored over another if the difference in its Bayesian log-evidence ∆ ln Z is greater than two, while strongly favored if it is greater than five (Trotta 2008) . If ∆ ln Z 2, then the models are indistinguishable. In this case, the model with fewer degrees of freedom would be chosen. Table 3 summarizes the results of our analysis on both serval-or TERRA-extracted RVs. We started by testing models with eccentric orbits for the transiting planets, however, given the current RV dataset, the eccentricity is not well constrained and such models were strongly disfavored (∆ ln Z = ln Z circ − ln Z ecc > 5). For this reason, we decided to assume simple circular orbits for both of the planets. As seen in Table 3 , including the two transiting planets in the model is favored against the fiducial model (0pl). On the other hand, we tested different types of two planet models. First, we considered just the two transiting planets (2pl), without accounting for additional noise sources. Then, we accounted for the stellar noise, modeling it in three different ways: with an exponential GP kernel (2pl+GP1), with an exponential sine-squared GP kernel (2pl+GP2) and with a simple sinusoid (2pl+sinusoid). All the tested two-planet models are statistically indistinguishable, with their Bayesian log-evidences within ∆ ln Z < 2.

However, for both serval and TERRA-extracted RVs, the nominal best model accounts for two circular orbits and an additional sinusoidal curve, whose period is equal to the stellar period of rotation we estimated through the long-term groundbased photometric data. For this test we imposed a normal prior on P rot , with a wide standard deviation (10). We additionally tried wide, uninformative priors for the period of the sinusoidal signal and we retrieved the same posterior distributions and log-evidences (Fig. 6 ) as for the test with a gaussian prior. With the RV analysis, we estimated a stellar period of rotation P rot = 35.2 +3.7 −3.0 d, consistent with the rotational period estimated from the ground-based long-term photometry in Section 4.2. Additionally, all models presented in Table 3 derive the same RV semi-amplitude for TOI-776 b and TOI-776 c, well within their 1σ uncertainties. This proves the robustness of the mass determination for the transiting planets, independently of the stellar noise distribution.

Leveraging the prior information on the stellar rotation from photometry discussed in Sect. 4.2 with the presence of a significant periodicity in the RV residuals of a two-planet model (Fig. 7b) , we decided to choose the 2pl+sinusoid as our final model for the joint fit. With respect to the RV extraction, although we achieved the nominal highest log-evidence using TERRA RVs, a slightly lower jitter is obtained with serval RVs (σ serval = 1.29 +0.38 −0.19 m s −1 ) compared to TERRA ones (σ TERRA = 1.37 +0.38 −0.25 m s −1 ). For this reason, we preferred to use the serval extracted RVs in the final joint fit.

We performed a joint fit using juliet of the TESS, LCO, and MEarth photometry and HARPS serval extracted RVs, using the 2pl+sinusoid model we selected after the RV-only analysis in Sect. 5.2.2. Table A .1 and 4 shows the priors and posteriors of all the fitted parameters. The data, residuals, and joint fit best Article number, page 10 of 24 model are shown in Figs. 3 and 8 for the photometry and the RVs, respectively. Table 5 lists the transit and physical parameters, derived using the stellar parameters in Table 2 .

As a sanity check, we performed an independent joint analysis of the transit photometry and Doppler measurements using the code pyaneti (Barragán et al. 2019) , which estimates the parameters of planetary systems in a Bayesian framework, combined with an MCMC sampling. We imposed uniform priors for all the fitted parameters. Following Winn (2010), we sampled for the mean stellar density ρ and recovered the scaled semi-major axis (R p /R ) for each planet using Kepler's third law. We found that the modeling of the transit light curves provides a mean stellar density of ρ = 5203 +1782 −1228 kg m −3 , which agrees with the density of 4834 +651 −559 kg m −3 derived from the stellar mass and radius presented in Sect. 4. As for the remaining parameters, the analysis provides consistent parameter estimates with those derived with juliet, corroborating our results.

The TOI-776 system consists of two transiting planets. The inner planet, TOI-776 b, has a period of 8.25 d, a radius of 1.83 ± 0.11 R ⊕ , a mass of 4.66 ± 0.97 M ⊕ , and a bulk density of 4.2 +1.3 −1.0 g cm −3 . The outer planet, TOI-776 c, has a period of 15.66 d, a radius of 2.06 ± 0.13 R ⊕ , a mass of 6.1 ± 1.5 M ⊕ , and a bulk density of 3.8 +1.2 −1.1 g cm −3 . The RV data show only one additional signal with a semi-amplitude of ∼ 2.3 m s −1 and a period of 35 d associated with the stellar rotation, as suggested by our analyses of the photometry and spectral line indicators.

While the occurrence rate of planets around early M dwarfs (3500 K < T eff < 4000 K) has been investigated in detail with Kepler and K2 samples (see e.g., Dressing & Charbonneau 2013 , 2015 Montet et al. 2015) , the number of currently known planets transiting low-mass stars is still much smaller with respect to those discovered around solar-type stars. While none of these surveys were optimized for M dwarfs, we expect more statistically significant results from the TESS mission for these stars. Figure 9 shows the confirmed transiting planets around M dwarfs as a function of the orbital period and the effective temperature of the host star. However, very few of these systems have precise determinations of the planetary masses (i.e densities), eccentricities and orbital architectures that would be required to link the statistical properties of this population with planet formation and evolution models in the low stellar mass regime.

There are several validated transiting systems similar to TOI-776 in terms of planetary architecture. Kepler-225, Kepler-236 and Kepler-231 are two-planet transiting systems composed of super-Earth and/or mini-Neptune sized companions with similar periods and semi-major axes, all validated by Rowe et al. (2014) . However, these systems are on average 5 mag fainter than TOI-776 and the planets do not have a mass determination nor precise stellar parameters. Similarly, K2-240 (Díez Alonso et al. 2018) has two transiting super-Earths with periods of 6 and 20.5 d, although they do not have mass determination and orbit an active star that is 2 mag fainter with a clear photometric rotational period of 10.8 d.

If compared to systems with mass determination, TOI-776 shows some similarities with Kepler-26 (Steffen et al. 2012 and 17.2 d, respectively, and bulk densities compatible with those of sub-Neptunes determined from transit timing variations. However, the system has two more planets without mass determination, an inner Earth-sized planet and an outer mini-Neptune sized planet. Kepler-138 is a very interesting system of three small planets, whose densities were estimated through photodynamical modeling (Almenara et al. 2018) . The most similar to the TOI-776 planets in terms of orbital period, Kepler-138 b (10.3 d) and c (13.8 d), are very different in composition, the former being a Mars analogue and the latter a prototypical rocky super-Earth. The third, outermost planet seems to have retained a substantial volatile-rich envelope. Finally, K2-3, the brightest of all three systems, has three small transiting planets and only the two inner ones (with periods of 10 and 24.6 d) have a mass determination using HARPS-N, HARPS, HIRES and PFS RVs (Almenara et al. 2015; Damasso et al. 2018; Kosiarek et al. 2019) , only an upper limit is measured for the third (with a period of 44.5 d). The planets have a similar composition, compatible to that of water-worlds or water-poor planets with gaseous envelopes, however the poor bulk density estimations of planets c and d impede further conclusions. The right panel of Fig. 9 shows all of the aforementioned systems, color-coded by bulk density and with the J-band magnitude of their host stars indicated.

Therefore, we conclude that, although multi-planetary systems of super-Earths and/or sub-Neptunes are common around early-type M dwarfs, only TOI-776 has all of its planets well characterized, with a bulk density measurement above 3σ, precise stellar parameters and a host star bright enough for atmospheric follow-up observations with current and planned facilities.

We investigated possible TTVs through a 3-body simulation, using the Python Tool for Transit Variations (PyTTV; Korth 2020). We simulated the estimated TTVs and RVs using the stellar and planetary parameters reported in Table 2 , 4 and 5 and found an expected TTV signal with a period of ∼ 150 d and a maximum amplitude of ∼ 2 min for the inner planet. Thus, the time span of the photometric observations, their cadence and signal-to-noise would prevent a detection of TTVs with the currently available data.

Additionally, we carried out a set of dynamical simulations to study the long-term stability of the system. We used the parameters in Table 4 and 5 and randomly drew 1000 samples from the posterior distributions as initial parameters for the dynamical simulations. We integrated each parameter set for 10 6 orbits of the inner planet, using the tool REBOUND (Rein & Liu 2012) with the standard IAS15 integrator (Rein & Spiegel 2015) . We also explored the stability using the MEGNO criteria as implemented in REBOUND. In the cases of close encounters between the bodies or one body ejection, the system would be flagged as unstable for the specific set of parameters. We found that the systems is dynamically stable over the entire integration time and for the whole parameter posterior space.

Moreover, since we could not constrain the eccentricity of the planets with the RV data, we tested whether an upper limit could be estimated from this analysis. For both planets, we randomly extracted values of orbital eccentricities from uniform distributions between 0 and 0.3, since RV analyses ruled out higher values. Unfortunately, our results do not provide better constraints than the RV data. Figure 10 shows the location of the TOI-776 system in a massradius diagram. Both planets occupy a scarcely populated region, characterized by a lack of planets around M dwarfs and with precise bulk density measurements. A comparison with the theoretical models by Zeng et al. (2016) , reported in the left panel of Fig. 10 , shows that TOI-776 b and c are consistent with mixtures of silicates and water in a 50-50 proportion. We adopted the three-layer models from Zeng & Sasselov (2013) and Zeng et al. (2016) to infer the interior structure of the planets. However, given the mass and radius input, the solution of the model is degenerate. As a consequence, the same mass-radius pair can lead to a broad range of combinations of iron, silicate and water mass fractions. On the other hand, when we applied the latest models by Zeng et al. (2019) , assuming a 1 mbar surface pressure level and an equilibrium temperature of 500 K (from Table 5), we found that an Earth-like rocky core with a 0.1% and a 0.3% molecular hydrogen atmosphere is consistent with the bulk densities of TOI-776 b and c, respectively. Nonetheless, it is clear that both of the planets in the system have an internal composition ranging from water worlds to rocky planets that have retained a significant atmosphere.

For a better understanding of the nature of the two exoplanets, we performed a more detailed modeling of their interior compositions, using their masses, radii and surface temperatures. Our model considers a canonical four-layer structure consisting of a two-component iron and silicate core, a layer of H 2 O and a H/He envelope. We assume that the core is Earth-like in composition (a third of iron, two-thirds of silicates by mass), meaning the core, water and H/He envelope mass fractions (x core , x H 2 O , x H/He ) are free parameters which sum to unity. The model solves the planetary structure equations of mass continuity and hydrostatic equilibrium assuming spherical symmetry. Further detail regarding the internal structure model can be found in Madhusudhan et al. (2020) and Nixon & Madhusudhan (in prep) .

The equation of state (EOS) prescriptions for the iron and silicate layers are adopted from Seager et al. (2007) , who used a Vinet EOS of the phase of Fe (Vinet et al. 1989; Anderson et al. 2001 ) and a Birch-Murnaghan EOS of MgSiO 3 perovskite (Birch 1952; Karki et al. 2000) . Thermal effects in these layers are ignored, since they have a small effect on the planetary radius (Howe et al. 2014) . However, thermal effects in the outer envelope can alter the mass-radius relation significantly (Thomas & Madhusudhan 2016) . 8 a purely terrestrial (iron plus rock) composition. Therefore, the planets must possess an envelope with some amount of H 2 O and/or H/He, in order to explain their masses and radii. The right panel of Fig. 10 shows limiting cases for each planet in which the envelope composition is either purely H 2 O or purely H/He. The mass and radius of TOI-776 b can be explained to within 1σ with a pure H 2 O envelope of 8-51% by mass or a pure H/He envelope with mass fraction 5.3 × 10 −5 -4.6 × 10 −3 . Best-fit solutions for pure envelopes are found at x H 2 O = 0.206 and x H/He = 6.6 × 10 −4 . TOI-776 c might have larger envelopes; within 1σ, it is consistent with a pure H 2 O layer of 18-73% or a pure H/He envelope with mass fraction 5.4×10 −4 -1.1×10 −2 . The best-fit pure-envelope solutions for TOI-776 c are x H 2 O = 0.52 and x H/He = 3.5 × 10 −3 . Each of the best-fit models, shown in the right panel of Fig. 10 , have a radiative-convective boundary at P rc = 10 bar. It is also possible that the planets in this system have both H 2 O and H/He components, as well as an iron/rock core. For the three components, we explored the full range of plausible values (x core , x H 2 O and x H/He ) that could explain the interior compositions of each planet. We considered two different temperature profiles for each planet, with P rc = 1 and 100 bar. Figure 11 shows the 1σ mass fractions of water and H/He compatible with the masses and radii of TOI-776 b and c. We imposed upper limits on the total H 2 O and H/He mass fractions for each planet: for TOI-776 b, x H 2 O ≤ 51% and x H/He ≤ 0.46%, while for TOI-776 c, x H 2 O ≤ 73% and x H/He ≤ 1.1%. These correspond to cases with pure H 2 O or H/He envelopes as previously discussed. Figure 11 shows as well a significant overlap between the bestfitting shaded regions for the two planets, meaning that the planets could also share the same composition.

The masses and radii of TOI-776 b and c allow for a wide range of possible solutions, from water worlds with steam atmo-spheres to mostly rocky planets with hydrogen-rich envelopes, however they are inconsistent with bare rocks without atmospheres. Our models assume a surface pressure of 0.1 bar, meaning a water-world solution for either planet yields a steam atmosphere. On the other hand, a higher surface pressure could result in liquid H 2 O at the surface. A rocky planet with an outgassed secondary atmosphere which includes carbon compounds is unlikely: Elkins-Tanton & Seager (2008) placed an upper limit on the mass fraction for this type of atmosphere at 5%. The lower mass limits in the case of pure H 2 O envelopes are 8% and 18% for TOI-776 b and c respectively. On the other hand, in a carbonrich atmosphere, the dominant species, CO 2 , has a higher mean molecular weight than H 2 O, leading to a lower atmospheric scale height. All things considered, we can infer that a 5% carbon-rich atmosphere is less than what would be needed to explain the planet radii. However, determining whether the two planets have H 2 O-or H/He-rich atmospheres is impossible with the present data. Atmospheric observations of the planets would be required in order to break this degeneracy.

The occurrence rate distribution of close-in planets exhibits a paucity of planets between 1.7−2.0 R ⊕ (Fulton et al. 2017; Hardegree-Ullman et al. 2020 ) around FGK stars (T eff > 4700 K) and between 1.4 − 1.7 R ⊕ (Hirano et al. 2018; Cloutier & Menou 2020) around mid-K to mid-M dwarfs (T eff < 4700 K). This feature is pointed out as the result of the transition from small rocky planets to larger non-rocky planets with volatile-rich envelopes (Weiss & Marcy 2014; Dressing & Charbonneau 2015) . Recent studies showed that the location of the radius gap depends on the orbital period or, alternatively, on the planet's insolation (Van Eylen et al. 2018; Martinez et al. 2019; Cloutier & Menou 2020) . Additionally, the width and center of the radius gap also depends on whether the host star is single or part of a multiple star system (Teske et al. 2018) .

According to the above discussion, if we consider the radius axis in Fig. 10 , TOI-776 b and c belong, within the uncertainties, to the radius gap in the case of FGK stars. On the other hand, they are well above the radius gap if we account for mid-K to mid-M dwarfs. Similarly, when looking at the distribution of transiting planets in a radius-insulation diagram (Fig. 12 , left panel), the TOI-776 planets lie above the radius valley -the two-dimensional view of the radius gap -that separates rocky super-Earths from volatile-rich sub-Neptunes around FGK stars.

The right panel of Fig. 12 shows the period-radius diagram of all known exoplanets which orbit M dwarfs with precise bulk density measurements. The dashed line marks the empirical location of the radius valley for FGK stars, following Van Eylen et al. (2018) , while the solid line indicates the location for mid-K to mid-M dwarfs as in Cloutier & Menou (2020) . The change in slope as a function of stellar type is the result of a change in the dominant mechanism responsible for sculpting the radius valley. For instance, the thermally driven mass loss, caused by photo-evaporation or core-powered mechanisms, becomes less efficient toward low-mass stars. The measured slope for mid-K to mid-M dwarfs suggests that gas-poor formation (Lee et al. 2014; Lee & Chiang 2016; Lopez & Rice 2018 ) might be the main process from which small planets form. However, thousands of small planets around low-mass stars with precise radii are needed in order to robustly state if the radius valley is the result of the erosion or the gas-poor formation scenarios . Although enriching the sample of exoplanet systems orbiting M dwarfs is nowadays possible thanks to TESS and fu- (2020), consistent with gas-poor formation scenarios. Together with TOI-776 b, the other two systems that are within both lines are K2-146 (solid, Hamann et al. 2019; Lam et al. 2020, translucent) and TOI-1235 (Bluhm et al. 2020; ture space-based missions such as PLATO an alternative is to obtain precise bulk density measurements of exoplanets lying in the region of discrepancy between models. TOI-776 b joins TOI-1235 b (Bluhm et al. 2020; ) and K2-146 b (Lam et al. 2020; Hamann et al. 2019) inside the period-radius region where thermally driven mass loss models disagree with the predictions from gas-poor formation. However, K2-146 b belongs to this region if we refer to the parameters reported in Hamann et al. (2019) , because the radius estimated by Lam et al. (2020) (see translucent points in Fig. 12 ) is more than 2σ higher, causing the planet to be placed outside the radius valley. Our previous analyses show that both TOI-776 b and c are likely to have retained a significant atmosphere, with slightly different envelope mass fractions. This result, given their period and radius, would be consistent with the planets having formed under gas-poor conditions.

On the other hand, the system's composition may be reconciled with thermally driven mass loss because the inner, most irradiated planet, has a smaller envelope mass fraction compared to its outer companion. Unlike other known systems whose planets straddle both sides of the radius gap (e.g., Dumusque et al. 2014; Niraula et al. 2017; Nowak et al. 2020) , TOI-776 is an interesting case where photo-evaporation could have stopped or become inefficient early in the planet's history. But, it is possible that the planets are currently undergoing mass loss under the core-powered mechanism, which erodes sub-Neptune planets into rocky super-Earths in Gyr timescales, contrary to the few Myr timescale when photo-evaporation is effective (Sanz-Forcada et al. 2011) . As reported in Table 2 , the age of TOI-776 is between 2 and 10 Gyr. However, the current data precision and limited number of known planets in this specific regime hamper any further investigation in favor of one or the other mechanism of formation. New studies on the dependence of the radius valley with other stellar parameters such as the age or metallicity, together with a larger sample of well-characterized planets in or near the radius valley, will help discerning between them in a demographic sense Berger et al. 2020; Gupta & Schlichting 2020) .

However, for the first time, we can compare between planets which belong to this region of the parameter space where formation models make opposing predictions. TOI-1235 b has a rocky composition with a 90% confidence upper limit in the envelope mass fraction of 0.5%, thus incompatible with a gaspoor formation scenario. We reach the opposite conclusion for TOI-776 b and c, whose bulk densities imply the presence of a volatile envelope making them compatible with the gas poor formation mechanism, and discarding photo-evaporation. Therefore, although other stellar parameters might need to be taken into account, we can tentatively predict that the stellar mass below which thermally-driven mass loss is no longer the main formation pathway for sculpting the radius valley is probably between 0.63 and 0.54 M , which correspond to the host stellar masses of TOI-1235 and TOI-776, respectively. More planets in this interesting region of the parameter space with precise bulk density measurements are key to reveal the mechanisms responsible of the radius valley emergence around low-mass stars with respect to solar-like stars.

6.5. Atmospheric characterization and habitability 6.5.1. Transmission Spectroscopy Metric (TSM)

We used the proposed metric by Kempton et al. (2018) acterization studies. Figure 13 shows the transmission spectroscopy metric (TSM) for all exoplanets in the Exoplanet Encyclopedia 7 with a radius less than 3 R ⊕ . We used the scale factors listed in Table 1 from Kempton et al. (2018) as opposed to the suggested value for temperate planets, 0.167, to compute the TSM values in Fig. 13 . The estimated TSM of TOI-776 b and c are 55.4 and 48.8 respectively, which places them among the top priority targets for atmospheric follow-ups of small planets around nearby stars. This is not surprising, because TOI-776 is one of the brightest M dwarfs with known transiting planets. However, most of the planets shown in Fig. 13 are well below the radius gap which makes the TOI-776 system a valuable target for atmospheric characterization in order to trace the formation and evolution of multi-planetary systems orbiting low-mass stars and break the degeneracy of internal composition models.

In order to quantitatively assess the possibility of TOI-776 b and c's atmospheric characterization with the James Webb Space Telescope (JWST), we investigated a suite of atmospheric scenarios and calculated their JWST synthetic spectra using the photo-chemical model ChemKM (Molaverdikhani et al. 2019a) and petitRADTRANS (Mollière et al. 2019) . We based the temperature structure of these planets on modern Earth's tempera-7 www.exoplanet.eu ture structure, and we increased the surface temperature for it to be consistent with the equilibrium temperature of TOI-776 b and c (Kawashima & Rugheimer 2019). We followed a similar approach as in Luque et al. (2019) : we estimated TOI-776's (T eff = 3709 K) flux in the range between X-rays and optical wavelengths, using as reference GJ 832 geometric mean spectra (T eff = 3816 K). The stellar data were obtained from the MUS-CLES database . To set up the models, we used Hébrard et al. (2012) chemical network with 135 species and 788 reactions, and an updated version of Hébrard et al.

(2012)'s UV absorption cross sections and branching yields. Figure 14 shows the synthetic transmission spectra of TOI-776 b and c assuming different metallicities, carbon-to-oxygen ratios and haze opacities. For our fiducial model (top left panel of Fig. 14) , we assume solar abundances. Such spectra are predominantly consisting of water and methane features, as expected for this type of planets (Molaverdikhani et al. 2019b ). The significance of these features are on the order of 100 ppm, well above conservative JWST expected noise floor (20 ppm for NIRISS and 50 ppm for MIRI, Greene et al. 2016) . We calculated the NIRISS-SOSS, NIRSpec-G395M, and MIRI-LRS uncertainties with PandExo (Batalha et al. 2017) , assuming two transits and binned for R = 50, supporting the previous statement. In this scenario, the contribution from haze opacity partially obscures molecular features below 2 µm, but it is almost ineffective at longer wavelengths (see left and right upper panels of Fig. 14) . We note, however, that the radiative feedback of haze particles Middle: Enhanced carbon-to-oxygen ratio by a factor of two. Bottom: Enhanced metallicity by a factor of 100. The left column represents spectra without haze opacity and the right column with haze opacity. might significantly affect the temperature structure and the composition of atmosphere (Molaverdikhani et al. 2020 ). We did not take this effect into account in this work in order to keep the temperature profiles consistent with the Earth's profile.

Smaller planets are expected to have enhanced metallicities (e.g. Wakeford et al. 2017 ). Therefore, we investigated two deviations from our solar abundance fiducial model: 1) an enhanced carbon-to-oxygen ratio (C/O) two-times the solar value, and 2) an enhanced metallicity of hundred times higher than solar. C/O enhancement alone does not affect the composition and spectral features substantially, as seen in Fig. 14 middle panels. On the other hand, one might expect a higher metallicity to result in more pronounced spectral features, due to higher species abundances. However, the bottom panels of Fig. 14 discard this possibility. On the contrary, an enhanced metallicity causes a higher mean molecular weight, which in turn shrinks the spectral significance (bottom-left panel of Fig. 14) , and, simultaneously, it results in a higher haze production, which also obscures the spectra significantly (bottom-right panel of Fig. 14) . Therefore, a flat transmission spectrum may indicate a hazy atmosphere with a high metallicity (Kreidberg et al. 2014 ) as opposed to a nonexisting atmosphere (Kreidberg et al. 2019) . Complementary observations, such as ground-based high-resolution spectroscopy or spectroscopy of the reflected light, are required to reveal the true nature of these flat spectra.

We present the discovery and characterization of the two-planet system transiting the early-type M dwarf TOI-776. Both planets were detected by the TESS mission, confirmed from groundbased transit follow-up observations and have their dynami-cal masses determined with precise RV measurements using HARPS. Our findings are summarized below:

-TOI-776 is a bright (V = 11.54 mag,J = 8.48 mag), M1 V star in the Centaurus constellation at 27.2 pc. Fifteen years of ground-based photometric monitoring by ASAS-SN, ASAS, NSVS, Catalina, and SuperWASP help us to measure a rotational period between 30 to 40 d, typical of inactive earlytype M dwarfs. -A joint fit of all the available transit photometry from TESS, MEarth, and LCOGT and the precise RVs from HARPS reveals that the TOI-776 system consists of two transiting planets, namely TOI-776 b, which has a period of 8.25 d, a radius of 1.83 ± 0.11 R ⊕ , a mass of 4.66 ± 0.97 M ⊕ , a bulk density of 4.2 +1.3 −1.0 g cm −3 , and an equilibrium temperature of 513±12 K; and TOI-776 c, which has a period of 15.66 d, a radius of 2.06 ± 0.13 R ⊕ , a mass of 6.1 ± 1.5 M ⊕ , a bulk density of 3.8 +1.2 −1.1 g cm −3 , and an equilibrium temperature of 415±10 K. The RV data show one additional signal, with a period of 35 d, associated with the star's rotation, in agreement from our analyses of the photometry and spectral line indicators.

-The bulk densities of TOI-776 b and c allow for a wide range of possible interior compositions, from water worlds to rocky planets with H/He-rich atmospheres, but they are too low for either planet to have a purely terrestrial (iron plus rock) composition. Thus, an atmosphere is expected for both planets. -From its location in a period-radius diagram, TOI-776 b lies in the transition region where formation and evolution models make different predictions for planetary systems orbiting M dwarfs. For the TOI-776 system, gas-poor formation may explain the composition of both planets although it is possible that the planets are still undergoing thermally driven mass loss under the core-powered scenario.

Article number, page 17 of 24 A&A proofs: manuscript no. main -The TOI-776 system is an excellent target for the JWST. It is the only known multi-planetary system with planets inside and near the radius valley for which all planets: 1) have a more than 3σ bulk density determination, and 2) are extremely suitable for atmospheric characterization. Thanks to the brightness of its host star, it is a remarkable laboratory to break the degeneracy in planetary interior models and to test formation and evolution theories of small planets around low-mass stars.

Acknowledgements. This work was supported by the KESPRINT 8 collaboration, an international consortium devoted to the characterization and research of exoplanets discovered with space-based missions. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. We acknowledge the use of TESS Alert data, which is currently in a beta test phase, from pipelines at the TESS Science Office and at the 

IEEE Transactions on Pattern Analysis and Machine Intelligence

celerite: Scalable 1D Gaussian Processes in C++, Python, and

Article number

A multi-planetary system

The Radial Velocity Method for the Detection of Exoplanets

Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun

Stellar Systems, and the Sun

Proc. SPIE

Tapir: A web interface for transit/eclipse observability, Astrophysics Source Code Library

IAU Symp. 61: New Problems in Astrometry

isochrones: Stellar model grid package

AAS/Division for Extreme Solar Systems Abstracts

Contemporary Physics

Exoplanet Transits and Occultations

Proceedings of the National Academy of Science

8 www.kesprint.science the NSF through grant AST-1824644. J.K., Sz.Cs., M.E., A.P.H., K.W.F.L., S.G. acknowledge support by DFG grants PA525/ 18-1, PA525/ 19-1, PA525/ 20-1, HA3279/ 12-1 and RA714/ 14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets. M.F, I.G. and C.M.P. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19 and 174/18 

Appendix A: Joint fit priors.

Article number, page 21 of 24 A&A proofs: manuscript no. main Table A .1. Priors used for the models presented in Sect. 5 using juliet. The prior labels of N, U, and J represent normal, uniform, and Jeffrey's distributions. The parameterization for (p, b) using (r 1 , r 2 ) (Espinoza 2018 ) and the linear (q 1 ) and quadratic (q 1 , q 2 ) limb darkening parameterization (Kipping 2013) Parametrization for e and ω.