TOI-663: A newly discovered multi-planet system with three transiting mini-Neptunes orbiting an early M star

We present the detection of three exoplanets orbiting the early M dwarf TOI-663 (TIC 54962195; V = 13 . 7 mag, J = 10 . 4 mag, R ⋆ = 0 . 512 ± 0 . 015 R ⊙ , M ⋆ = 0 . 514 ± 0 . 012 M ⊙ , d = 64 pc). TOI-663 b, c, and d, with respective radii of 2 . 27 ± 0 . 10 R ⊕ , 2 . 26 ± 0 . 10 R ⊕ , and 1 . 92 ± 0 . 13 R ⊕ and masses of 4 . 45 ± 0 . 65 M ⊕ , 3 . 65 ± 0 . 97 M ⊕ , and <5.2 M ⊕ at 99%, are located just above the radius valley that separates rocky and volatile-rich exoplanets. The planet candidates are identified in two TESS sectors and are validated with ground-based photometric follow-up, precise radial-velocity measurements, and high-resolution imaging. We used the software package juliet to jointly model the photometric and radial-velocity datasets, with Gaussian processes applied to correct for systematics. The three planets discovered in the TOI-663 system are low-mass mini-Neptunes with radii significantly larger than those of rocky analogs, implying that volatiles, such as water, must predominate. In addition to this internal structure analysis, we also performed a dynamical analysis that confirmed the stability of the system. The three exoplanets in the TOI-663 system, similarly to other sub-Neptunes orbiting M dwarfs, have been found to have lower densities than planets of similar sizes orbiting stars of different spectral types.


Introduction
Small exoplanets discovered around other stars could be similar or very different to those of our Solar System in terms of composition, formation history, and atmospheric properties.They are also difficult to study due to the small transit and radial-velocity (RV) signals they produce.But the improved sensitivity of recent surveys, such as that conducted with the Transiting Exoplanet Survey Satellite (TESS), have enabled the detection of smaller exoplanets.TESS was launched in 2018 to conduct an all-sky survey and discover transiting planets around the nearest and brightest stars (Ricker et al. 2015).
When looking for small and cool exoplanets, low-mass stars have distinct advantages and are particularly interesting.The transit depth and RV semi-amplitude caused by small exoplanets orbiting M dwarfs are much greater than those caused by similar planets orbiting larger stars of earlier spectral types, making such planets easier to detect.More importantly, such systems are excellent candidates for atmospheric characterization via transmission or thermal emission spectroscopy (e.g., Kempton et al. 2018;Batalha et al. 2018).To determine which planets are predominantly rocky or volatile-rich and ⋆ Nasa Sagan Fellow.
whether they could have a significant atmosphere, precise radius measurements from transit surveys must be combined with mass measurements from dynamical observations, such as precise RV measurements.Exoplanets with precise mass, radius, and bulk density measurements orbiting M dwarfs are thus important objects for better understanding the so-called radius valley separating the population of super-Earths and mini-Neptunes (e.g., Fulton et al. 2017;Mayo et al. 2018;Cloutier & Menou 2020;Hardegree-Ullman et al. 2020).
The radius valley around M dwarfs is located at slightly smaller radii than around Sun-like stars (Fulton et al. 2017;Cloutier & Menou 2020) and marks the transition between rocky super-Earths and mini-Neptunes that host H/He envelopes (Owen & Wu 2017).A variety of physical mechanisms have been proposed to explain the emergence of the radius valley, including photo-evaporation, core-powered mass loss (Jin & Mordasini 2018;Gupta & Schlichting 2020;Wyatt et al. 2020), or a direct outcome of the planet formation process wherein super-Earths form at later times in either a gas-poor or completely gasdepleted environment (Lee & Connors 2021;Lopez & Rice 2018).Each of these mechanisms predicts a unique dependence of the radius valley on the orbital period such that measuring the radius valley's slope provides direct insight into the origin of the radius valley.Cloutier & Menou (2020) discovered that the radius valley separating rocky and gaseous planets exhibits a negative slope with insolation around low-mass stars, as opposed to the positive slope found around solar-type stars.This supports models of direct terrestrial planet formation in a gas-poor environment.The different slopes of models of gas-depleted formation and thermally driven mass loss naturally carve out a subset of the period-radius parameter space, within which the two models make opposing predictions about the bulk composition of planets.Characterizing the bulk composition of transiting planets allows us to determine whether the M dwarf radius valley originated directly from planet formation or as a result of post-formation atmospheric escape.
Multiple transiting planet systems are particularly interesting: they provide a unique opportunity for direct comparative planetology (Millholland et al. 2017) because the exoplanets formed within the same protoplanetary disk and evolved around the same host star.These systems also allow studies of planetplanet interactions (Barros et al. 2015) and planetary formation and migration (Delisle 2017).The discovery and accurate characterization of a system with multiple transiting planets thus represents a significant addition to the current knowledge about planetary systems formation.Furthermore, transiting multi-planet systems may contain additional long-period planets not yet detected.Deep searches to find additional exoplanets have already been conducted for a few stars, and both photometry and RV follow-up of transiting systems have successfully detected additional transiting exoplanets (e.g., the Trappist-1 system; Gillon et al. 2017;Luger et al. 2017 andthe LHS 1140 system;Lillo-Box et al. 2020).
Here we present the discovery of three new exoplanets orbiting around the same M1-type star.This discovery was made using TESS, multiple ground-based photometry facilities, including ExTrA and the Las Cumbres Observatory (LCO) network, and precise RV measurements from the High Accuracy Radial velocity Planet Searcher (HARPS), MAROON-X, the InfraRed Doppler (IRD), and the Echelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (ESPRESSO).Section 2 presents an analysis of the stellar properties of TOI-663.Section 3 describes the observations and data used in this study, including photometry from space and from the ground, RV measurements, and high-resolution imaging.Section 4 presents the global analysis of the available data, which is used to constrain the planetary properties.Section 5 presents the discussion and conclusion.

Stellar parameters
TOI-663 is an M-type star located at a distance of 64.23 ± 0.18 pc (Gaia Collaboration 2018, Lindegren et al. 2018, and Bailer-Jones et al. 2018).The astrometry, photometry, and stellar parameters are reported in Table 1.
TOI-663 mass and radius are derived from the stellar parallax and K s -band magnitude, used to compute the absolute K s -band magnitude M K s , and the empirically derived M dwarf massluminosity and radius-luminosity relations from Mann et al. (2019) and Mann et al. (2015), respectively.The analysis by Mann includes stars with close-to-solar metallicities.We used SpecMatch-Emp1 (Yee et al. 2017) and HARPS high-resolution spectra of TOI-663 to estimate the stellar metallicity, for which we found a value around 0.07 ± 0.12 dex (but the five best matches are spread over ∼0.5 dex so the metallicity remains We determine the stellar radius of TOI-663 independently by modeling the spectral energy distribution (SED) with stellar atmosphere and evolution models.To construct the SED, we used the magnitudes from Gaia Data Release 2 (DR2; Evans et al. 2018;Maíz Apellániz & Weiler 2018), the 2-Micron All-Sky Survey (2MASS; Skrutskie et al. 2006;Cutri et al. 2003), and Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010;Cutri & et al. 2013).We modeled these magnitudes using the procedure described in Díaz et al. (2014).We used informative priors for the effective temperature (T eff ), and metallicity ([Fe/H]) from the analysis of the HARPS co-added spectra with SpecMatch-Emp (Yee et al. 2017), and for the distance from Gaia DR2 (Gaia Collaboration 2018;Bailer-Jones et al. 2018).Non-informative priors were used for the rest of parameters.We used the PHOENIX/BT-Settl stellar atmosphere models (Allard et al. 2012).The data with the maximum a posteriori stellar atmosphere model is shown in Fig. 1.We obtained a stellar radius of 0.504 ± 0.010 R ⊙ , consistent with the estimation from Mann et al. (2015) presented in Table 1.We decided to use the stellar radius estimated from the empirical relation for the rest of the paper.
We tried to determine the stellar rotation period of TOI-663.We studied activity indices from the HARPS spectra to see if there was some periodicity in the S index of the sodium lines or in the S index of the Hα line.After removing any long-term variation, we did not find a significant peak in the periodograms.The flux around the calcium lines is very low in the HARPS A19, page 2 of 18 spectra but using the ESPRESSO spectra, we were able to measure the value of log R ′ HK = −4.783± 0.063, derived from the calcium doublet of the ESPRESSO combined spectra.This indicates that the star has likely a rotation period estimated to P rot = 30 ± 3.1 days (Astudillo-Defru et al. 2017a).According to this rotation period and the stellar radius, we expect the vsin(i) to be less than a km s −1 , which is consistent with the values of full width at half maximum (FWHM) of the cross-correlation function (CCF) from the ESPRESSO spectra.No sign of periodicity were found in the All Sky Automated Survey (ASAS; Pojmanski 2002) and the Zwicky Transient Facility (ZTF; Bellm et al. 2019) photometry as well as in the HARPS dataset (see Fig. 2).Unfortunately, this difficulty to find the rotation period of the star prevented us from determining the age of the system using gyrochronology.NASA Ames processed the TESS photometric data and the resulting Presearch Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014) ) light curve of TOI-663 was corrected for dilution in the TESS aperture by known contaminating sources.It is critical to ensure that no visually close-by targets that could affect the depth of the transit are present in the 21 ′′ TESS pixel and to check for contaminating eclipsing binaries.Figure 3 shows a plot of the target pixel file and the aperture mask that is used for the simple aperture photometry (SAP).We can see that one star overlaps the TESS aperture but its faintness results in minimal dilution of the TESS light curve.The SPOC conducted a transit search of the sector 9 light curve on 2019 April 25 and the combined light curve for sectors 9 and 35 on 2021 May 27 with an adaptive, noise-compensating matched filter (Jenkins 2002;Jenkins et al. 2010Jenkins et al. , 2020)).The sector 9 search produced threshold crossing events with orbital periods of 2.60d and 4.70d; the combined search produced a third threshold crossing event with period of 7.10d.An initial limb-darkened transit model was fitted (Li et al. 2019) and a suite of diagnostic tests were conducted to help make or break the planetary nature of each signal (Twicken et al. 2018).The transit signatures passed all diagnostic tests presented in the SPOC Data Validation reports, and all TIC objects other than the target star were excluded as the transit source in each case.The TESS Science Office reviewed the vetting information and issued alerts for TOI-663.01 and TOI-663.02 on 2019 May 07, and TOI-663.03 on 2021 July 08 (Guerrero et al. 2021).

Observations
Canto Martins et al. (2020) conducted a search for rotation period in TESS TOIs using fast Fourier transform, Lomb-Scargle, and wavelet techniques, accompanied by a rigorous visual inspection, and found that TOI-663 does not show any sign of chromospheric activity throughout the observed time ranges.

Ground-based photometry
The TESS pixel scale is ∼21 ′′ pixel −1 and photometric apertures typically extend out to roughly 1 arc-minute, generally causing multiple stars to blend in the TESS aperture.To definitively exclude the presence of another star causing the signal A19, page 3 of 18 Cointepas, M., et al.: A&A, 685, A19 (2024) 2.

ExTrA photometry
ExTrA (Bonfils et al. 2015) is a near-infrared (0.85 to 1.55 µm) multi-object spectrograph fed by three 60-cm telescopes located at La Silla Observatory in Chile.Seven full transits (and three partial transits) of TOI-663.01, six full transits of TOI-663.02, and two full transits of TOI-663.03 were observed using the 2 https://tess.mit.edu/followup ExTrA telescopes (see Table 2).We used 8 ′′ aperture fibers and the low-resolution mode (R ∼ 20) of the spectrograph with an exposure time of 60 seconds.Five fiber positioners are used at the focal plane of each telescope to select light from the target and four comparison stars chosen with 2MASS J-magnitudes (Skrutskie et al. 2006) and Gaia effective temperatures (Gaia Collaboration 2018) similar to the target.The resulting ExTrA data were analyzed using a custom data reduction software to produce a synthetic photometry in a 0.85-1.55micron bandpass, described in more detail in Cointepas et al. (2021).

LCOGT
We observed nine and four full predicted transit windows of TOI-663.01 and TOI-663.02,respectively, using the A19, page 4 of 18 Cointepas, M., et al.: A&A, 685, A19 (2024) Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network nodes.We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen 2013), to schedule our transit observations.The 1 m telescopes are equipped with 4096 × 4096 SINISTRO cameras having an image scale of 0. ′′ 389 per pixel, resulting in a 26 ′ × 26 ′ field of view.The images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018).Differential photometric data were extracted using AstroImageJ (Collins et al. 2017).We detected transit-like signals in all 13 light curves using photometric apertures that exclude flux from all known neighbors of TOI-663, confirming the event on-target relative to all known Gaia DR3 and TICv8 stars.

MuSCAT
We observed one transit of TOI-663.02 on 2020 March 24 UT with the multiband imager MuSCAT (Narita et al. 2015) on the National Astronomical Observatory of Japan (NAOJ) 188 cm telescope at the Okayama Astro-Complex in Okayama, Japan.MuSCAT has three optical channels each equipped with a 1k × 1k CCD camera with a pixel scale of 0. ′′ 361 pixel −1 , enabling simultaneous g-, r-, and z s -band imaging.The observation was performed with exposure times of 30, 15, and 30 s for the g, r, and z s bands, respectively.We used optical diffusers to obtain slightly defocused images.We performed image calibration and aperture photometry to extract relative light curves using a custom pipeline described in Fukui et al. (2011).

MuSCAT2
We observed one transit of TOI-663.01 on 2020 January 17 UT with the multi-band imager MuSCAT2 (Narita et al. 2019) on the 1.52m TCS telescope at the Teide Observatory in the Canary Islands, Spain.MuSCAT2 has four optical channels each equipped with a 1k × 1k CCD camera with a pixel scale of 0. ′′ 44 pixel −1 , enabling simultaneous g-, r-, i-, and z s -band imaging.The observation was performed with exposure times of 60, 30, 30, and 15 s for the g, r, i, and z s bands, respectively.We performed image calibration and aperture photometry to extract relative light curves using a dedicated MuSCAT2 pipeline described in Parviainen et al. (2020).

GMU
TOI-663.01 was observed on the night of 2021 February 07 in the R band with 90 second exposures with the George Mason University (GMU) 0.8m Ritchey-Chretien telescope using the automated software control from Reefe et al. (2022).Intermittent clouds were present throughout the night, and a detection of the transit was not obtained.Data were reduced with a custom Python code and plate-solved with astrometry.net.Aperture photometry, light curve extraction, and detrending were performed with AstroImageJ (Collins et al. 2017).The results are available at ExoFOP at the NASA Exoplanet Archive.

OAAlbanya
TOI-663.02 was observed on the night of 2020 March 19 UT in the Ic band using an exposure of 120 s and with a good focus.The observation was performed from the 0.4 m Schmidt-Cassegrain telescope from Observatori Astronòmic Albanyà, at an altitude of 230 m in the northeast of Spain.The night was mostly transparent during the window of observation, with some minor clouds at the out of transit phase of the observation.Data reduction was performed with aperture differential photometry using the AstroimageJ suite.Results showed a tentative signal coherent with an on-time transit, with high RMS related to the predicted depth, but consistent with an event nevertheless.The dataset was detrended using the airmass values and the movement of the Y-axis of the image sequence to reduce systematics.

PEST
We observed TOI-663.(Demory et al. 2020).The observations were acquired with a deep-depleted Andor Ikon CCD in the z' filter and consisted of 396 data points, which cover a full transit of TOI-663.01.Each observation had an exposure time of 2.0 s spanning 5 h of continuous observations.The data were reduced and the aperture photometry was done, using the custom pipeline PRINCE done on site, as detailed in Demory et al. (2020).For the final light curve, the aperture of 14 pixels was chosen.We then analyzed the light curve with AstroImageJ (Collins et al. 2017) to correct for airmass and total counts in the aperture thus reducing systematics.ESPRESSO observations can be carried out with simultaneous calibration using a Fabry-Pérot etalon on a second calibration fiber.However, given the relatively faint magnitude of TOI-663 (V = 13.67)we did not use the Fabry-Pérot etalon for simultaneous calibration, and the calibration fiber was pointed at the sky for background calibration that is used as a sky subtraction calibration in the data processing.The ESPRESSO RVs were processed using the Data & Analysis Center for Exoplanets (DACE), a facility dedicated to extra-solar planets data visualization, exchange, and analysis hosted by the University of Geneva 5 .The RV and stellar activity indicators provided in DACE are computed by cross-correlating the calibrated spectra with stellar templates for the specific spectral type, using the latest version (3.0.0) of the publicly available ESPRESSO Data Reduction Software6 .The ten ESPRESSO RVs of TOI-663 have a median uncertainty of 1.58 m s −1 .

MAROON-X
We obtained 14 spectra of TOI-663 with MAROON-X (Seifahrt et al. 2018) at Gemini North in two campaigns from UT 2021 February 22 to 2021 March 04 and UT 2021 April 19 to 2021 April 30 under Prog-ID GN-2021A-Q-230.MAROON-X is a fairly new Extreme Precision Radial Velocity spectrograph operating in the red-optical wavelength regime, delivering R ≃ 85 000 from 500-920 nm in two arms, dubbed blue (500-670 nm) and red (650-920 nm).Spectra were taken with a fixed exposure time of 20 min and showed an average peak signal-to-noise ratio (S/N) of 100 in the red arm.The data were reduced by the instrument team using a custom python3 data reduction pipeline to produce optimum extracted and wavelength calibrated one-dimensional spectra.The RV analysis was performed using serval (Zechmeister et al. 2018), a template matching RV retrieval code in a custom python3 implementation.On average, the RV uncertainty per data point was 2.3 m s −1 for the blue and 1.7 m s −1 for the red arm.MAROON-X uses a stabilized Fabry-Perot etalon for wavelength and drift calibration (Seifahrt et al. 2022); it can deliver 30 cm s −1 on-sky RV precision over short timescales (Trifonov et al. 2021) but suffers from inter-run RV offsets, which in the case of our observations had a prior of 1.5 ± 0.5 m s −1 for the blue and 1.6 ± 0.5 m s −1 for the red channel.

IRD
We observed TOI-663 with the IRD instrument (Tamura et al. 2012;Kotani et al. 2018) on the Subaru 8.2 m telescope.IRD simultaneously covers the Y, J, and H bands with the spectral resolving power of λ/∆λ ≈ 70 000.A total of six near-infrared spectra were secured between UT 2021 February 02 and June 08, each with an integration time of 1200−1500 s.One frame was interrupted by bad weather, which we discarded in the subsequent analysis.We reduced the raw frames with the standard reduction procedure, and measured the relative RVs using the custom RV-analysis pipeline for IRD (Hirano et al. 2020).Except for the frame affected by bad weather, the extracted one-dimensional spectra had the typical S/Ns of 51−73 per pixel at 1000 nm, and derived RVs had the statistical errors of 3.5−4.8m s −1 .

High-resolution imaging
As part of our standard process for validating transiting exoplanets to assess the possible contamination of bound or unbound companions on the derived planetary radii (Ciardi et al. 2015), we observed the TOI 663 with near-infrared adaptive optics (AO) imaging at Keck and with optical speckle imaging at Gemini.Gaia DR3 is also used to provide additional constraints on the presence of undetected stellar companions as well as wide companions.

Near-infrared adaptive optics imaging
The Keck Observatory observations were made with the NIRC2 instrument on Keck-II behind the natural guide star AO system (Wizinowich et al. 2000) on 2019 June 10 UT in the standard three-point dither pattern that is used with NIRC2 to avoid the left-lower quadrant of the detector, which is typically noisier than the other three quadrants.The dither pattern step size was 3 ′′ and was repeated twice, with each dither offset from the previous dither by 0.5 ′′ .NIRC2 was used in the narrow-angle mode with a full field of view of ∼10 ′′ and a pixel scale of approximately 0.0099442 ′′ per pixel.The Keck observations were made in the narrowband Br-γ filter (λ o = 2.1686; ∆λ = 0.0326 µm) and in the Jcont filter (λ o = 1.2132; ∆λ = 0.0198 µm), both with an integration time of 10 seconds for a total of 90 seconds on target per filter.
The science frames were flat-fielded and sky-subtracted.The flat fields were generated from a median average of dark subtracted flats taken on-sky.The flats were normalized such that the median value of the flats is unity.The sky frames were generated from the median average of the dithered science frames; each science image was then sky-subtracted and flat-fielded.The reduced science frames were combined into a single combined image using an intra-pixel interpolation that conserves flux, shifts the individual dithered frames by the appropriate fractional pixels, and median-coadds the frames.The final resolutions of the combined dithers were determined from the FWHM of the point spread functions: 0.052 ′′ at 2 µm and 0.064 ′′ at 1.2 µm.
The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 20 • at separations of integer multiples of the central source's FWHM (Furlan et al. 2017).The brightness of each injected source was scaled until standard aperture photometry detected it with 5σ significance.The resulting brightness of the injected sources relative to TOI-663 set the contrast limits at that injection location.The final 5σ limit at each separation was determined from the average of all of the determined limits at that separation and the uncertainty on the limit was set by the rms dispersion of the azimuthal slices at a given radial distance.The Keck data have a sensitivity close-in of ∆.2 mag at 0.06 ′′ (Br-γ) and δ.4 mag at 0.07 ′′ (Jcont), and deeper sensitivity at wider separations (δ.7 mag at ≳0.5 ′′ ); the final sensitivity curves A19, page 6 of 18 for the Keck Br-γ imaging is shown in (Fig. 4).No close-in (≲1 ′′ ) stellar companions were detected by Keck in agreement with the optical speckle imaging.

Optical speckle imaging
TOI-663 was observed on 2019 June 10 and 2020 January 12 UT using the 'Alopeke and Zorro speckle instruments on Gemini South and North, respectively.Zorro and 'Alopeke (Scott et al. 2021) both provide simultaneous speckle imaging in two bands (562 and 832 nm) with output data products including a reconstructed image and robust contrast limits on companion detections (Howell et al. 2011).Both datasets gave consistent results and Fig. 5 shows the Gemini South Zorro resultant 5sigma contrast curves and the reconstructed 832-nm speckle image.We find that TOI-663 is a single star with no companion brighter than about 5-9 magnitudes, respectively, from 0.1 to 1.2 arcsec.At the distance of TOI-663 (64 pc), these angular limits correspond to spatial limits of 6.4-77 AU.

Gaia assessment
In addition to the high-resolution imaging, we used Gaia to identify any wide stellar companions that may be bound members of the system.Typically, these stars are already in the TESS Input Catalog and their flux dilution to the transit has already been accounted for in the transit fits and associated derived parameters.There are no additional widely separated companions identified by Gaia that have the same distance and proper motion as TOI-663 (see also Mugrauer & Michel 2020, 2021).
Additionally, the Gaia DR3 astrometry provides additional information on the possibility of inner companions that may have gone undetected by either Gaia or the high-resolution imaging.The Gaia renormalized unit weight error (RUWE) is a metric, similar to a reduced chi-square, where values that are ≲1.4 indicate that the Gaia astrometric solution is consistent with the star being single, whereas RUWE values ≳1.4 may indicate an astrometric excess noise, possibly caused by the presence of an unseen companion (e.g., Ziegler et al. 2020).TOI-663 has a Gaia DR3 RUWE value of 1.6 -while this value is slightly higher than the nominal value of 1.4, the 1.4 RUWE boundary is more of a guideline as opposed to a hard cutoff.Additionally, TOI-663 does not appear in the Gaia DR3 binary table; as such, the Gaia astrometric fit is still consistent with a single star model.However, this RUWE value could be linked to the long-term RV trend seen in the HARPS dataset, which could be associated with a bound companion.

Analysis and modeling
We used the software package juliet (Espinoza et al. 2019) to model the photometric and RV data.The algorithm is built on many publicly available tools for the modeling of transits (batman; Kreidberg 2015), RVs (radvel; Fulton et al. 2018), and Gaussian processes (GPs; george; Ambikasaran et al. 2015;celerite, Foreman-Mackey et al. 2017).In order to compare different models, juliet efficiently computes the Bayesian evidence (ln Z) using dynesty (Speagle 2020), a python package to estimate Bayesian posteriors and evidence using nested sampling methods.Instead of starting with an initial parameter vector centered on a likelihood maximum discovered through optimization techniques, nested sampling algorithms sample directly from the given priors.Throughout our various analyses, we ensured that we had enough live points given the number of free parameters to avoid missing peaks in the parameter space.Here we decided to conduct two independent analyses of our data.We started the analysis using only TESS photometry, and the resulting planet parameters were used as priors in a subsequent RV analysis.As the result of these two analyses, we were able to constrain the priors of some parameters and use them for a joint analysis of all our data.

TESS photometry
First, using juliet, we modeled the TESS PDCSAP light curve where our planet candidates was initially detected.The transit model fits the stellar density ρ ⋆ along with the planetary and jitter parameters.Except for the stellar density, which we calculated in Sect.2, we used ExoFOP priors for the orbital parameters.We adopted a few parametrization modifications when dealing with the transit photometry.We assigned a quadratic limb-darkening law for TESS, as shown to be appropriate for space-based missions (Espinoza & Jordán 2015), which was then parameterized with the uniform sampling scheme (q 1 , q 2 ) introduced by Kipping (2013).Additionally, rather than fitting directly for the planet-to-star radius ratio (p = R p /R ⋆ ) and the impact parameter of the orbit (b = a/R ⋆ cos i), juliet used the parameterization introduced in Espinoza (2018) and fit for the parameters r 1 and A19, page 7 of 18  We assume circular orbits.The TESS PDC-SAP light curves are detrended from contamination, so we set the TESS dilution factor to one.In addition, we added in quadrature a jitter term σ TESS to the TESS photometric uncertainties, which may be underestimated due to additional systematics in the space-based photometry.Figure 6 depicts the best-fit transit model for the two TESS sectors.In the same figure, we also showed the phase-folded light curve for each planet.

Radial velocity analysis
For the RV modeling using juliet, we tested a three-planet and four-planet model including or not a linear long-term trend and GPs and we compared the evidence to select the one with the highest probability.We used as prior for each candidate the period and time-of-transit center estimated from the TESS photometry analysis and assume circular orbits.After fitting all of our datasets with a three-planet model including a trend seen in the HARPS dataset, we saw that the model did not fit properly our set of RV measurements for the different instruments.As a model with only three planets and a trend does not seem to match our data, we tested two different models: adding a GP or considering the possibility of having a fourth non-transiting exoplanet in the system.We looked for another periodic signal in the RV datasets.We found a potential peak at longer period around 25.8 days in the HARPS dataset compatible at the ∼1.4sigma level with the estimated rotational period derived from the log R ′ HK value (see Sect. 2).We fitted a four-planet model to our RV datasets (considering this new planet is not transiting and with a period range between 8 and 30 days) and found that the evidence of this model (ln Z = 293.33 ± 0.66) was slightly lower than the evidence for the model with 3 exoplanets, a long-term trend and a GP (ln Z = 294.60 ± 0.67), although the model comparison does not favor one model over the other (Kass & Raftery 1995).For the joint fit of photometry and RVs we adopted the model with three exoplanets, a long-term trend and a GP (Matern kernel with non-informative loguniform priors), as the signal of a possible fourth exoplanet was not convincing and the GP allows a better description of the residuals of the three-planet model.The long-term trend was estimated at 0.0229 +0.0091 −0.0085 m day −1 in the joint fit presented in the next section.

Joint fit of all data
We performed a joint analysis using juliet of the TESS photometry complemented by some ground-based photometry from ExTrA, LCOGT, and MuSCAT and velocity data from HARPS, ESPRESSO, MAROON-X, and IRD to obtain the most precise parameters of the TOI-663 system.Each individual groundbased light curve presented in Sect. 3 was analyzed for the transit time variation (TTV) analysis (see Sect. 4.5).From this, we decided to only select a few of these transits for the global analysis to improve on the radius measurement.To choose which transits to use, we kept a first selection of transits for which the posterior of the time-of-transit center was more constrained after the fit (a nonuniform posterior on t 0 ), and then we only selected the nights with the smallest systematics visible in the data.We selected all transits available for TOI-663.03 as it was the least constrained by the TESS photometry.
After a verification that the depth of the transits were consistent in the different wavelength ranges, we decided to fit only one set of r 1 and r 2 parameters for each planet.To detrend the photometry data, we used time-dependent GPs.We chose the approximate Matern kernel proposed by Foreman-Mackey et al. ( 2017) because there are no obvious quasi-periodic oscillations in the light curves.We do not have a sufficient precision on the RV data to constrain the eccentricity of the different planets in this system so we decided to fix them at 0 to decrease the number of free parameters.Multi-planet systems tend to have low eccentricities but this decision is explained in more detail in the analysis in Sect.4.4.We seek the posterior distribution of 127 free parameters in total, and due to the large number of parameters, we perform this fit using dynesty.We constrained our priors using previous juliet results before the final run for the joint modeling to optimize the computation time due to the large parameter space.Because we used rather broad priors for the planetary parameters in the separate fits of TESS and HARPS data, nested sampling is an efficient method for exploring the parameter space, and because most of the planetary parameters are specific to a given data type, they would not change significantly in a joint fit.
Table 3 shows the posteriors of the planetary parameters derived using the stellar parameters presented in Sect.2, and Table C.1 shows the posteriors of the jitter parameters.Figures 7-9 show the results of the photometry of our joint fit to the data for the three exoplanets.In Figs. 10 and

Stability of the system
From our global fit, we find that the three planets currently have period ratios P c /P b =1.807 and P d /P c =1.513.We investigated the short-term dynamics in this three-planet system.Our A19, page 8 of 18 Cointepas, M., et al.: A&A, 685, A19 (2024)  T eq (K) ( †)  674 ± 30 553 ± 24 482 ± 21 Notes.The posterior estimate corresponds to the median value.Error bars denote the 68% posterior credible intervals. (* ) We sample from a normal distribution for the stellar mass, stellar radius, and stellar temperature that is based on the results from Sect. 2. ( †) Equilibrium temperature was calculated assuming 0.3 Bond albedo and the semimajor axis distance considering an uniform temperature across the entire planet's surface.The map has a resolution of 121 × 121 pixels, which represent as many system configurations that were numerically integrated to derive a NAFF chaos estimation as defined in Eq. ( 1).We report the results in the form of a color code, where the bluer regions of the map refer to systems with smaller NAFF, and hence weaker chaos.The nearby principal MMR are indicated, and the vertical line marks the orbital period of planet c as estimated from our global fit.
objective is both to get a first hint on the range of orbital eccentricities compatible with the system stability, and to look for potential mean-motion resonances (MMRs).As such, we built a two-dimensional chaos map.It consists in a 121 × 121 grid of simulations in the parameter space exploring the orbital eccentricity of the middle planet e c on one axis and the orbital period of that same planet P c on the other axis (Fig. 12).All the other dynamical parameters were initially fixed to their estimation resulting from the global fit, hence with the hypothesis that e b = 0 = e d and ω i = 0 for i = b,c,d.Each subsequent box in the grid defines a unique set of system orbital parameters that we used as initial conditions for numerical integrations of the equations of motion.We used the IAS15 adaptive timestep integrator to perform the numerical integrations (Rein & Spiegel 2015).IAS15 is part of the python package REBOUND7 (Rein & Liu 2012).We also included a correction for general relativity effects following Anderson et al. (1975) and implemented in REBOUNDx (Tamayo et al. 2020).Each set of initial conditions was integrated over 10 kyr, during which the mean longitudes of the planets were regularly recorded in output.After the numerical integrations, we used the mean longitude time series to derive the numerical analysis of fundamental frequencies (NAFF) fast chaos indicator (Laskar 1990(Laskar , 1993)).This quantity informs about the strength of chaos in a planetary system via the variations of the orbital mean-motions, according to the following formula: In this equation, n 1 and n 2 depict the orbital mean-motion computed over the first and second half of the integration, respectively.n 0 is the initial mean motion.The mean-motion difference between the two integration halves is computed on each planet j, and the maximum value defines the NAFF of the planetary system.We refer to Stalport et al. (2022) for a more detailed synthesis of this process, and for a discussion of the link between NAFF and the orbital stability.In Fig. 12, we report the results as a color code, where the bluer regions indicate smaller NAFF, and hence, weaker chaos.Orbital eccentricities e c larger than 0.1 are strongly disfavored, given the resulting strong chaos.While this exact threshold should be taken with caution since it results from a partial exploration of the parameter space, it provides us with insights into the range of eccentricity allowed by orbital stability.Only small values are not destructive for the system.This justifies the assumption of null eccentricities to favor the convergence of our global fit.Additionally, two vertical features are detected.They overlap with the 9:5 MMR between planets b and c, and the 3:2 MMR between planets c and d.A vertical line is drawn at the orbital period of planet c as derived from our global fit, which shows that the system does not fall in any of those features.We do not find the planets in any low-order MMR.
In the previous section, the possibility of a fourth planet in the TOI-663 system is mentioned.This additional candidate would have a significantly longer orbital period (more than three times the period of planet d), hence likely not affecting the orbital stability of the inner system.

TTV analysis
Using juliet, we calculated transit timing variations (TTVs) for all photometric datasets.Rather than fitting a period P and a time of transit center t 0 , juliet searches for individual transit times.We combined multiple transits happening at the same time in order to reduce the uncertainties on the estimation of the time of transit center.The results of the analysis showing the difference between the observed transit times and the calculated linear ephemeris from TESS transits and for the different planets are presented in Fig. 13.The data show no significant variation at

Discussion and conclusion
In this paper we present the validation and characterization of three mini-Neptunes transiting the M dwarf TOI-663.The exoplanets were detected by the TESS mission, then confirmed and characterized from ground-based transit follow-up observations and precise RV measurements.

Mass-radius diagram and internal structure
TOI-663 b, c, and d, with respective radii of 2.27 ± 0.10 R ⊕ , 2.26 ± 0.10 R ⊕ , and 1.92 ± 0.13 R ⊕ and masses of 4.45 ± 0.65 M ⊕ , and 3.65 ± 0.97 M ⊕ , 1.9 ± 1.3 M ⊕ , lie at the edge of the population of mini-Neptunes.However, the determination of the mass of planet d is at the ∼1.7 sigma level, which is not sufficient to precisely characterize this exoplanet.We can only provide an upper limit for planet d of <5.2 M ⊕ at 99%.We plan to observe TOI-663 in the coming semesters using a combination of HARPS and the Near Infra Red Planet Searcher (NIRPS; Bouchy et al. 2017) to characterize the long-term trend found in the HARPS dataset, improve the mass measurements, and better characterize possible additional signals in this system.The three detected planets of the TOI-663 system are lowmass mini-Neptunes.Their radii are considerably larger than rocky analogs (Fig. 14), and hence their radii must be dominated by volatiles, for example water.
To estimate the range of compositions for TOI-663 b, TOI-663 c, and TOI-663 d, we employed a layered interior model consisting of up to four layers: a H-He atmosphere, a water layer, a silicate mantle, and an iron core (see Dorn et al. 2017).We solved the standard structure equations for each layer and used a nested sampling algorithm to construct structure models that reproduce the measured masses and radii of the planets, to quantify the degeneracy between various interior parameters, and to produce posterior probability distributions (Buchner et al. 2014) assuming a solar Mg/Si ratio (Lodders 2020) and an age of 5 Gyr.For the water layer, we used the AQUA equation of state (Haldemann et al. 2020).In the atmosphere layer, we adopted the semi-gray atmosphere model of Guillot (2010).This model takes the stellar irradiation and the intrinsic luminosity of the planet into account; the latter was estimated using the scaling relation from Mordasini (2020).
In a first scenario, we allowed the water mass fraction to be up to 50% underneath a H-He layer.The results are summarized in the upper half of Table 4. Within these constraints, we find that the three TOI-663 planets have a similar composition of roughly 30% core, 40% mantle, and 30% water layer mass fraction, with a small H-He atmosphere of log(M atm /M p ) = −4 to −5.
In order to estimate the maximum water mass fraction in the planets, we also consider internal structures with no H-He atmospheres (M atm = 0) and without limiting the total water mass fraction.Additionally, we followed Dorn & Lichtenberg (2021) in assuming an isothermal temperature profile until 0.1 bar, followed by an adiabatic profile.The results for this case are summarized in the lower half of Table 4.We find that the water mass fraction is about 20% higher than in the constrained case (i.e., approximately 50%).Accordingly, the core and mantle mass fractions are reduced by ∼10% (i.e., M core /M p ∼ 20% and M mantle /M p ∼ 30%).
Overall, the TOI-663 planets are expected to be water-rich and may have small H-He atmospheres.It should be noted, however, that while disfavored by the nested sampling algorithm, an Earth-like composition with a small H-He atmosphere (<1%) can also fit the observed mass and radius of the TOI-663 planetary system.Moreover, pollution of the H-He envelope could lead to a contraction of the atmosphere and hence a higher atmospheric mass fraction (Lozovsky et al. 2018).In addition, we have estimated here the amount of surface water and neglected any water dissolved deep in the interior, which could result in higher water mass fractions (Dorn & Lichtenberg 2021).Also, the host star elemental mass fractions can differ from solar  2022) proposed a population of planets around M dwarfs that are water-rich worlds.On one hand, the planets of TOI-663 may indeed be part of water-rich worlds, according to our analysis.On the other hand, the planets do not follow the density trend of the proposed water-rich population as obtained from Zeng et al. (2019), which, however, neglects density variations due to differing equilibrium temperatures or ages.Furthermore, Rogers et al. (2023) show that the water-rich worlds can be equally well described by compositions that include H-He envelopes.In general, additional constraints are needed to further reduce the existing degeneracy between hydrogen-and water-dominated envelopes, for example, spectral data from JWST.
From a formation point of view, significant water budgets on the TOI-663 planets are feasible.High water budgets on the TOI-663 planets may be inherited from a formation beyond the snow line (Mordasini et al. 2009) as well as the production of endogenic water (Kimura & Ikoma 2020).How much endogenic water can be produced by redox reactions between a primordial H-He envelope in contact with an oxidized mantle for such small sub-Neptunes as the TOI-663 planets is still poorly known, but recent studies indicate that this is an important aspect to consider (Misener et al. 2023).
Despite evidence that M dwarfs tend to form rockier planets, the mass-radius diagram of known exoplanets with precise masses (<25%) and radii (<8%), shown in Fig. 14, shows that sub-Neptunes around M dwarfs appear to be less dense than those around other stars.The three exoplanets of the TOI-663 system more or less conform to this observation.

Position around the radius valley
Using data from the Kepler and K2 missions, Cloutier & Menou (2020) discovered that the radius valley separating rocky and gaseous planets has a negative slope with insolation around low-mass stars, as opposed to the positive slope found around solar-type stars (Martinez et al. 2019).This supports models of the direct formation of gas-poor terrestrial planets whose atmospheres were not removed through photo-evaporation (Owen & Wu 2017;Jin & Mordasini 2018;Lopez & Rice 2018), core-powered mass loss (Ginzburg et al. 2018;Gupta & Schlichting 2019, 2020), or impact erosion (Schlichting et al. 2015;Wyatt et al. 2020).In contrast to Cloutier & Menou (2020), Van Eylen et al. (2021) show the opposite trend of the slope of the radius valley with orbital period, which is similar to the positive slope found around solar-type stars.To resolve the discrepancy, they state that a larger sample of well-studied planets orbiting M dwarf stars is required.Unfortunately, the three TOI-663 exoplanets are located above the radius valley and do not help us constrain either of the two models.

Similarity of multi-planet systems
M dwarfs tend to host multi-planet systems, which can provide much stronger tests for planetary formation theories than singleplanet systems (Fang & Margot 2012;Steffen & Hwang 2015;Ballard & Johnson 2016).Weiss et al. (2018) have demonstrated a correlation between the size of planets in a multi-planet system and a regular spacing between them, described as the "peas in the pod" theory.This concept has been employed to restrict planet formation models.According to Otegi et al. (2022), who examined multi-planet systems observed by K2 and TESS for which both the radii and the masses were measured precisely, M-dwarf multi-planet systems exhibit more similarity in terms of density than in mass, indicating a possible connection in their formation mechanism.The densities of the three exoplanets of the TOI-663 system are consistent with each other within 1 sigma and therefore compatible with the Otegi et al. (2022) observation.

Potential for atmospheric characterization
Regarding atmospheric characterization, we calculated the transmission spectroscopic metric (TSM) in the J magnitude using the equation from Kempton et al. (2018) and obtained the following values: TSM = 76 +16 −12 for TOI-663 b, TSM = 75 +29 −17 for TOI-663 c, and TSM = 78 +120 −35 for TOI-663 d.This value takes into account the equilibrium temperature of the exoplanet computed at its semimajor axis distance.Unfortunately, they all fall below the criteria of Kempton et al. (2018) for being the best candidates for transmission spectroscopy.

CointepasFig. 1 .Fig. 2 .
Fig. 1.SED of TOI-663.The solid line is the maximum a posteriori PHOENIX/BT-Settl interpolated synthetic spectrum, red circles are the absolute photometric observations, and gray open circles are the result of integrating the synthetic spectrum in the observed bandpasses.

Fig. 3 .
Fig. 3. TESS target pixel file image of TOI-663 in Sector 9 (created with tpfplotter; Aller et al. 2020).The electron counts are colorcoded.The red bordered pixels are used in the SPOC SAP.The size of the red circles indicates the Gaia DR2 magnitudes of all nearby stars.

Fig. 4 .Fig. 5 .
Fig. 4. Keck near-infrared AO imaging and sensitivity curve for TOI-663 taken in the Br-γ filter.The image reaches a contrast of ∼7.7 magnitudes fainter than the host star within 0. ′′ 5. Inset: image of the central portion of the data.
Fig. 6.Modeling of the TESS data.Two left panels: TESS photometry time series from Sectors 9 and 35 along with the maximum a posteriori (solid black line).Three right panels: phase-folded light curve to the period of the planet (black points correspond to one-sixth of the transit duration for each planet).The three panels are associated with the b, c, and d planets, respectively.The dark area to the left on the far right plot is due to the transit of the other exoplanets that happened close to the transit of TOI-663 d.

r 2
to guarantee full exploration of physically plausible values in the (p,b) plane.
11, we present the RV component of the joint fit and the phase-folded data.The corner plots of the RV semi-amplitude for the three planets are presented in Fig. B.1.

Fig. 7 .
Fig. 7. Phase-folded light curves for TOI-663 b for the TESS (17 transits), ExTrA (three transits with two telescopes each), and LCOGT (four transits) photometry.The black line is the maximum a posteriori model, the black points are data binned to one-sixth of the transit duration.

Fig. 8 .
Fig. 8. Phase-folded light curves for TOI-663 c for the TESS (eight transits), ExTrA (two transits with two telescopes each), MuSCAT (one transit with three filters), and LCOGT (two transits) photometry.The black line is the maximum a posteriori model, the black points are data binned to one-sixth of the transit duration.

Fig. 9 .
Fig. 9. Phase-folded light curves for TOI-663 d for the TESS photometry (seven transits) and ExTrA photometry (two transits, one with one telescope and the other with two telescopes).The black line is the maximum a posteriori model, the black points are data binned to one-sixth of the transit duration.The dark area to the left on the far right plot is due to the transit of the other exoplanets that happened close to the transit of TOI-663 d.

Fig. 10 .Fig. 11 .
Fig. 10.RV measurements from HARPS, ESPRESSO, MAROON-X, and IRD for TOI-663.The data points with errors are shown in colors; the black line and gray band correspond respectively to the median and 1σ of 1000 randomly chosen posterior samples from the joint fit.

Fig. 12 .
Fig.12.Chaos map around our best-fit solution.The map has a resolution of 121 × 121 pixels, which represent as many system configurations that were numerically integrated to derive a NAFF chaos estimation as defined in Eq. (1).We report the results in the form of a color code, where the bluer regions of the map refer to systems with smaller NAFF, and hence weaker chaos.The nearby principal MMR are indicated, and the vertical line marks the orbital period of planet c as estimated from our global fit.

Fig. 14 .
Fig.14.Mass-radius diagram of the detected exoplanets(Otegi et al. 2020) and the TOI-663 planetary system.The green, gray, and dark blue lines show the mass-radius relation for an Earth-like, a pure-iron composition, and a 50% water content with and Earth-like interior, respectively.The blue region shows the M-R relation of pure water for the equilibrium temperature range of the planets in the TOI-663 system.
Mann et al. (2015Mann et al. ( , 2019) )ar radius and obtained a more precise effective temperature than the one computed from theMann et al. (2015Mann et al. ( , 2019) )empirical relations.

Table 2 .
Summary of ground-based follow-up observations.
01 and TOI-663.02 on 2020 April 28 UT and 2020 April 26 UT, respectively, with the Perth Exoplanet Survey Telescope (PEST) in the Rc band.PEST is located near Perth, Australia.At the time, the 0.3 m telescope was equipped with a 1530 × 1020 SBIG ST-8XME camera with an image scale of 1. ′′ 2 pixel −1 resulting in a 31 ′ × 21 ′ field of view.A custom pipeline based on C-Munipack 3 was used to calibrate the images and extract the differential photometry.
3.2.8.SAINT-EX Time series observations of the transit of TOI-663.01 were acquired with the SAINT-EX Observatory (Search And char-acterIsatioN of Transiting EXoplanets) on the night of 2020 February 05 UT.SAINT-EX is a 1-m F/8 Ritchey-Chrétien telescope that resides within the Observatorio Astronómico Nacional located at the Sierra de San Pedro Mártir in Baja California, Mexico (31.04342N,115.45476W) at an altitude of 2780 m (Pepe et al. 2021pectra of TOI-663 with ESPRESSO(Pepe et al. 2021) at the 8.2 m ESO Very Large Telescope (VLT) array, at the Paranal Observatory in Chile.The observations were obtained from 2022 April 05 to 2022 April 18 for an ESPRESSO program designed to test if the density of sub-Neptunes orbiting M-dwarfs is different than sub-Neptunes orbiting more massive stars (Program ID: 109.2391.001;PI: Grieves).We observed TOI-663 in HR2x1 mode (1 UT, R ∼ 140 000) over a spectral range from ∼380 to ∼780 nm.ESPRESSO is contained in a temperature and pressure controlled vacuum vessel to minimize the night spectral drift.

Table 3 .
Posterior stellar and planetary parameters obtained from our joint photometric and RV juliet analysis for TOI-663.

Table 4 .
Nested sampling results for a maximum water mass fraction of 50% (top) and assuming no H-He atmosphere and no limit on the water mass fraction (bottom).
Table A.1.RV measurements for the different spectrographs used in this analysis.